Jean Motsch
email: jean.motsch@st-cyr.terre.defense.gouv.fr

Olivier D Éforges
email: olivier.deforges@insa-rennes.fr

Marie Babel
email: marie.babel@insa-rennes.fr

EMBEDDING MULTILEVEL IMAGE ENCRYPTION IN THE LAR CODEC

A still image codec should not only to be good from compression point of view, it might also provide services. The LAR interleaved S+P codec allows loss-less to lossy coding, presenting both resolution and distortion scalability. The proposed method brings image encryption at no cost in the LAR framework. It exploits the embedded quadtree decomposition of the LAR to provide multilevel protection to the whole bit-stream. Theoretical aspects are considered and experimental results show the effectiveness of the process.

INTRODUCTION

Evolution in still image codecs shows the need for new embedded services. The JPEG2000 standard [START_REF] Taubman | Jpeg2000: Image Compression Fundamentals, Standards, and Practice[END_REF] introduces some features that are among the ones requested for a good codec:

• error resilience for transmission in noisy environments, such as wireless and the Internet,

• region of interest coding,

• scalability, both in resolution and distortion,

• and, of course, good distortion-bit-rate performances. Data protection, through encryption for example, is also a need. Good codecs should provide not only the best distortion-bit-rates performance, but also include the aforementioned services.

This paper deals with embedded image protection. Whereas most methods relies on coding whole or part of the image, using block or stream ciphering, we propose a method based on using the quadtree decomposition bit-stream as a way to protect the content of the image. The main idea is to transmit the data without the quadtree decomposition, using the quadtree as the key to decrypt the image. Figure 1 gives a rough idea of the control that can be done on the quadtree bit-stream. As the quadtree bit-stream is multilevel, the authorization system manages different levels of quality for the same original bit-stream. As an example, low resolution of the picture is always provided to take the decision of buying the picture while highest resolution is furnished after payment. As for most data encryption method, one way to confirm the usefulness of this method is to show that a brute force attack on the bit-stream is too time consuming.

This paper is organized as follows. Section 2 presents shortly the context of the LAR framework and its relationship to quadtree. Section 3 emphasizes on some theoretical results on quadtree partitioning and their coding cost. Section 4 shows some experimental results and comments, and section 5 does conclude the topic.

LAR FRAMEWORK

LAR codec

The LAR (Locally Adaptive Resolution), based on a variable blocksize decomposition, leads to an efficient lossy image compression technique [START_REF] Deforges | Region of interest coding for low bit-rate image transmission[END_REF].

The LAR compression method -initially designed for lossy greyscale image coding -is a two layer codec : a spatial coder and a complementary spectral one. The spatial coder provides a low bitrate compressed image whereas the spectral coder codes the texture. Self extraction of regions of interest is also possible [START_REF] Deforges | Region of interest coding for low bit-rate image transmission[END_REF]. Extensions also enables middle and high quality coding and region-level based chromatic image coding [START_REF] Deforges | Supervised segmentation at low bit rates for region representation and color image compression[END_REF]. Furthermore, subjective quality image evaluation show that the low resolution LAR picture is better than the same bit-rate JPEG2000 coded picture [START_REF] Deforges | Region of interest coding for low bit-rate image transmission[END_REF].

A minimal-redundancy pyramidal decomposition, the LAR-APP introduces a multiresolution framework [START_REF] Babel | Lossless and lossy minimal redundancy pyramidal decomposition for scalable image compression technique[END_REF]. Finally, an interleaved S+P pyramidal decomposition with a refined prediction model shows state-of-the-art results combined to resolution and SNR scalability [START_REF] Babel | Interleaved s+p pyramidal decomposition with refined prediction model[END_REF].

The basic principle of the spatial coder is that the local resolution, i.e. the pixel size, depends on the local activity. This leads to build a resolution variable image using a quadtree structure. The pixel's size (e.g. from 32 × 32 to 2 × 2) is computed thanks to a local morphological gradient. A single parameter, the threshold on the contrast, drives the resolution process. Figure 2 presents an example of quadtree decomposition. Through this decomposition, the pixel size provides implicitly the block's nature: small blocks are located on contours whereas large ones describes smooth areas. An efficient post-processing [START_REF] Muresan | Optimal recovery approach to image interpolation[END_REF] removes perceptible blocks artifacts.

LAR Interleaved S+P bit-stream

The LAR interleaved S+P codec is a two passes codec:

1. first pass constructs a low bit-rate, low resolution and high distortion picture. It plays the role of the spatial coder and it can be stopped at every level of the pyramid.

2. second pass, equivalent to the spectral coder, provides the texture at every level of the pyramid.

Figure 3 shows the pyramidal decomposition on a picture. This two passes process allows to refine the picture and provides both resolution and distortion scalability.

In order to reflect the pyramidal decomposition, the bit-stream is organized in sub-streams, one for each level of the pyramid. There are three different type of sub-streams: the quadtree's refinement, the refinement on the low resolution image (first pass) and the local texture values. Figure 4 gives a basic idea of the LAR interleaved S+P bit-stream. bit-stream

¡ ¡ ¡ ¡ ¡ ¡ ¡ ¢ ¡¢ ¢ ¡¢ ¢ ¡¢ ¢ ¡¢ ¢ ¡¢ ¢ ¡¢ £ ¤£ ¤£ ¤£ £ ¤£ ¤£ ¤£ £ ¤£ ¤£ ¤£ £ ¤£ ¤£ ¤£ ¥ ¤¥ ¤¥ ¤¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¤¥ ¤¥ ¤¥ ¥ ¤¥ ¤¥ ¤¥ ¦ ¡¦ ¡¦ ¦ ¡¦ ¡¦ ¦ ¡¦ ¡¦ ¦ ¡¦ ¡¦ ¦ ¡¦ ¡¦ ¦ ¡¦ ¡¦ ¦ ¡¦ ¡¦ § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § § ¡ § ¤¤¤¤¤¤ ¤¤¤¤¤ ¤¤¤¤¤ ¤¤¤¤© ¤© ¤© ¤© ¤© ¤© © ¤© ¤© ¤© ¤© ¤© © ¤© ¤© ¤© ¤© ¤© © ¤© ¤© ¤© ¤© ¤© ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¤ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ¡ ! ¡! ¡! ¡! ¡! ! ¡! ¡! ¡! ¡! ! ¡! ¡! ¡! ¡!

Error on the quadtree

The LAR interleaved S+P shows to be reactive on errors made on its embedded quadtree decomposition. Figure 5 illustrates two important properties of the codec for data encryption. Firstly, without a perfect knowledge of the quadtree, the reconstructed picture is far from the original one. Secondly, guessing the top level quadtree decomposition doesn't help obtaining a good quality picture. We consider an image I having M lines and N columns with S = M N pixels. Every pixel is coded using q bits and therefore take Q = 2 q different values. QP [L T ..L B ] (I) is a quadtree partition of the image I using blocks size ranging from LT × LT for the top level to LB × LB at the bottom. ˛QP [L T ..L B ] (I) ˛is the number of different quadtree partitions possible for image I.

Let lT = log2LT , lB = log2LB and l = lT -lB +1 the number of levels in the quadtree partition. We note Ω [L T ..L B ] the number of different quadtree possible on a LT × LT block. ω [L T ..L B ] = log2Ω [L T ..L B ] is the number of bits needed to code this quadtree partition.

For sake of illustration, we will consider a 256 × 256 8 bits picture using a QP [32..2] partition with 5 levels.

The number of quadtree partition Ω [L T ..L B ] is related to the following recursively defined function Φ:

Φ(0) = 1 (1) Φ(n) = 1 + Φ 4 (n -1) (2)
In fact, Ω [L T ..L B ] = Φ(l -1). Table 1 gives the first values of Φ(n). Notice that log 2 Φ(n) is the number of bits required to code the value of Φ(n). For n > 0, Φ(n) can be roughly approximated by Φa(n) = 2 (4 n-1 ) . This establishes the quadtree coding cost The number of quadtree partition for the whole image I is simply computed considering that every LT × LT block contains a QP [L T ..L B ] partition. ˛QP [L T ..L B ] (I) ˛is then given by:

C(QP [L T ..L B ] ) to about 4 l-2 bits per square LT × LT block. n Φ(n) log 2 Φ(n) Φa(n) 0 1 0 1 1
˛QP [L T ..L B ] (I) ˛= Φ(l -1) M N L 2 T (3)
Usually, coding cost are expressed in bit per pixel (bpp). This cost is given by:

log 2 Φ(l -1) L 2 T ( 4
)
Following with our example, ˛QP [32..2] (I) ˛≈ 2 4186 . It also means that the maximum coding cost will be 4186 bits or 0.064 bpp.

Quadtree coding

Coding of quadtree partitions is computed using a progressive method. A one is transmitted for every block that has to be splitted in 4, while a zero is transmitted otherwise. And this is done recursively. Figure 6 shows an example of quadtree. The quadtree partitions are QP [8. .2] . The bit-stream to code the first level is 0111. The second, and last, level, is coded with 011010010010. The bit-stream is 16 bits long, or 0.0625 bpp. In fact, there are Φ2 4 = 83521 possible quadtree partitions that need at most 16.35 bits to be coded, or 0.0634 bpp. This is the upper bound of the coding cost.

Remarks

Section 3.1 demonstrates the huge number of quadtree partitions that are possible in coding an image. However, two important remarks have to be presented.

Firstly, the result shown in equation ( 3) relies on the hypothesis that blocks are independent in the image. This assumption implies that the coding cost of the quadtree given by equation ( 4) is a maximum coding cost. If there is some correlation among blocks, the entropy of the generated bit-stream will be less than the coding cost. An entropic coder will then be able to reduce the size of the bitstream canceling the remaining redundancy. Experimental results of section 4 show that the entropy after entropic coding is still important.

Secondly, the space of all images is still order of magnitudes bigger than the space of partitions. In the example we used, every pixel is coded with 8 bits but the maximum quadtree coding cost is 0.064 bit. That explains that a quadtree cannot be used as a unique identifier of an image. Indeed, numerous images share the same quadtree.

EXPERIMENTATION

Experiments are conducted using pictures lena, airplane, baboon, goldhill, man, pepper and woman, all of size 512 × 512 and 8 bits coded. Quadtree decomposition are QP [64..2] i.e. 6 levels in the pyramid. In that case, the maximum quadtree entropy is 0.06387 bpp.

Figure 7 shows the influence of the threshold on the entropy of the quadtree bit-stream. The threshold drives the process of quadtree decomposition. That is, as the threshold increases, the number of regions decreases, with bigger blocks. Usually, the threshold has a fixed value of 30 (on a 256 greyscale). Notice that the quadtree entropy remains above 0.03 bpp, about half its maximum value.

If the quadtree entropy is 0.03 bpp, for a 512 × 512 picture, the quadtree entropy is about 7864 bits for the whole image. That means, the number of possible quadtrees is about 2 7864 ≈ 10 2367 . Even if decoding a bit-stream takes 1 µs, then the brute force approach will be largely too long.

Conversely, if the protection of the picture may stay for 100 years, that means, 3.156 × 10 9 s, the entropy of the quadtree must be only more than 51.5 bits. As the quadtree is progressively transmitted in the bit-stream, the entropy of each level of the quadtree coding is computed and table 2 shows the results obtained with threshold value equal to 30. The entropy of first level is null when there is no 64 × 64 blocks in the quadtree. The second level has an entropy about 50 that might be sufficient to provide a reasonable data protection. The entropy increases as the block's size decreases, making the reconstruction without the partition more and more difficult. Notice also that for image airplane, the 17 bits of entropy are not sufficient. In all case, the top level of the pyramid might be decoded, furnishing only a very low resolution picture with high distortion.

CONCLUSION AND PERSPECTIVES

This paper presented an original image encryption method with following properties: embedded in the original bit-stream at no cost, allowing multilevel access authorization combined with state-of-theart still picture codec. Multilevel quadtree decomposition provides a way to select the quality of the picture decoded. Theoretical and experimental results shows the effectiveness of this tool on an information theory point of view.

One main goal of this work is the implementation of an archiving system for high resolution art pictures of the Louvre's museum. This digital library will provide a selective access with different quality of images.

Figure 1 :

 1 Figure 1: Quadtree bit-stream control

Figure 2 :

 2 Figure 2: Example of LAR quadtree decomposition with threshold equals to 30

Figure 3 :

 3 Figure 3: Pyramidal decomposition in LAR interleaved S+P

Figure 4 :

 4 Figure 4: LAR interleaved S+P bit-stream description

Figure 5 : 3 . 1 .

 531 Figure 5: LAR interleaved S+P reconstructed pictures

Figure 6 :

 6 Figure 6: QP [8..2] on a 8 × 32 image

Figure 8 :

 8 Figure 8: Quadtree entropy vs. number of blocks

Table 1 :

 1 Φ(n) and Φa(n) values for small n For our example, Ω [32..2] = 4.86 × 10 19 and we need up to 65.40 bits to code the quadtree partition on an 32 × 32 block.

		2	1	2
	2	17	4.09	2 4
	3	83522	16.35	2 16
	4 4.86 × 10 19	65.40	2 64

  Figure8shows results on investigating on the quadtree entropy as a function of the number of blocks. This figure is shaped like the h2(p) = -p log 2 (p) -(1 -p) log 2 (1 -p). It exhibits the fact that the quadtree entropy is maximum approximately when the number of blocks is half the maximum number of possible blocks.

	quadtree entropy [bpp]	Evolution of the quadtree entropy vs. threshold
	0.07
	0.06
	0.05
	0.04
	0.03
	0.02
	0.01
	0.00							threshold
	0	50	100	150		200	250	300
	Figure 7: Quadtree entropy vs. threshold
	entropy [bpp]		entropy vs. number of blocks
	0.07
	0.06
	0.05
	0.04
	0.03
	0.02
	0.01
	0.00							#blocks
	0	10000	20000	30000	40000	50000	60000	70000

Table 2 :

 2 Entropy of quadtree levels given in [bits] for various 512 × 512 images

	image	64 × 64 32 × 32 16 × 16 8 × 8 4 × 4
	lena	0	61	656	2916	8886
	airplane	17	175	601	1922	7243
	baboon	0	0	21	1647 10356
	goldhill	0	65	218	2659 12218
	man	0	35	466	2824 11014
	pepper	0	41	667	3148	8245
	woman	0	86	560	2810	9483


P Leclaire

M J Swift

K V Horoshenkov

Determining the specific area of porous acoustic materials from water extraction data

The internal pore surface area per unit mass of solid (specific area) is an important parameter for the physical study of sound wave propagation in porous media. It is shown that this parameter can be evaluated with a reasonable accuracy from the integration of the water extraction curve (volume of water extracted versus applied pressure). This method is nonacoustic, modelless and well adapted to acoustic materials. It is a complementary method to the standard BET (Brunauer, Emmet, Teller) technique, generally used for powders in chemistry and the pharmaceutical industry. Experiments are made on loose glass beads of a few hundred microns diameter and results are compared to theoretical values. Results on heterogeneous materials used in civil engineering and in noise control are also given. The parameters necessary for the description of modem acoustic materials are estimated from the specific area measured by this method.

I. INTRODUCTION

The propagation of sound waves in porous media is a complex phenomenon involving elastic and kinetic energy exchanges between the solid and the fluid and visco-thermal interactions. This topic is studied in the area of geophysics, petroleum engineering and underwater acoustics and in the fields of noise control and civil engineering, where the saturating fluid is air. A current problem which is receiving attention is the description of the high frequency behavior without the need to make assumptions about the pore shape and using appropriate parameters, i.e., measurable parameters with a physicai meaning. Johnson et a/. 1 and Champoux and Allard 2 -4 have expressed the high frequency behavior in terms of overall coefficients integrating the complexity of the pore shape; these are the viscous length A and the thermal length A 1 • Yamamoto and Turgut 5 and Attenborough 6 have shown that the size of the pore can affect the results and have included a distribution in pore size in the acoustic models. The parameters defined in all these models can be measured under certain conditions and for certain materials. In the case of polyurethane foams, simple experiments using ultrasound in air or in different gases have been proposed to measure the tortuosity 7 and the viscous and thermal lengths 8 • 9 but these methods are limited to high porosity materials producing a fairly low sound attenuation. In addition, scattering of ultrasonic waves 10 due to the small wavelengths involved can be a limiting factor.

In this article, we show that the parameters involved in the acoustic description of porous media can be evaluated from the measurement of specific area. The specific area, defined as the internal pore surface area per unit mass of solid, is obtained from a nonacoustic method applicable to many materials without making any assumption on the pore shape. In particular, the experiments lead to a direct determination of A 1

• This parameter is useful for the study, at the pore scale, of modem acoustic materials from the physical standpoint and is of a major importance for the description of • thermal exchanges associated to the sound wave propagation in air-filled porous structures. [START_REF] Lafarge | [END_REF] This method using extraction of water has never been used in acoustics for measuring the specific area.

II. CAPILLARITY AND SPECIFIC AREA

A. Surface tension and capillary law

Water rising in a narrow tube by the capillary effect is a very well known phenomenon (see Lowell and Shields 11 for example). It is described by the concepts of surface tension and contact angle, and a capillary law. The surface tension is used to describe the mechanical tension of the interface separating two fluids while the capillary law gives the condition of mechanical equilibrium in a capillary system. For an interface defined by its two principal radii of curvature r 1 and r 2 , the capillary law is

fl.P = (]' ( ~ + ~), , ri r2 (1)
where fl.P is the pressure difference between the two fluids and (]' the surface tension in N/m. The curvature of the interface depends on the pore geometry and the contact angle 0.

For a narrow circular tube of radius r the interface is a portion of the surface of a sphere and r 1 = r 2 = r/cos 0. The pressure difference fl.P is then given by 2 (]'cos 0 fl.P=---

r (2)
The capillary law for slit-like pores is obtained when either r 1 or r 2 tends to infinity in Eq. ( 1):

a cos() dP=--, a
where a is the semiwidth of the slits.

B. Specific area from the water extraction curve

(3) Rootare and Prenzlow [START_REF] Rootare | [END_REF] have proposed a method for measuring the pore surface area of various powders from mercury intrusion data and compared their results to those provided by the standard BET method (Brunauer, Emmet and Teller 13 ). The method of Rootare and Prenzlow is based on the calculation of the work dW required to raise a liquid in a capillary. This work is equal the work required to force the liquid out [START_REF] Lowell | Powder Surface Area and Porosity[END_REF] and the method can be adapted and applied in the case of water extraction. The contact angle for wetting liquids being between oo and 90°, cos () is greater than zero and

dW=awaCOS () dA, ( 4
)
where a wa is the surface tension of the water/air interface and dA is the surface area covered by the liquid along the pore wall during the process. This work is also equal to the product P dV of the volume dV of displaced liquid and the increase of pressure P required to move the liquid. Thus the following relation is obtained:

awa COS ()dA =P dV. ( 5
)
Integrating Eq. ( 5) from zero to the total pore volume VTot yields the total pore surface area:

S= PdV. 1 fvtor awa COS() o (6)
Equations ( 4)- (6) show that no particular assumption on the pore shape is made for the calculation of the specific area.

Ill. EXPERIMENTS

Evans [START_REF] Evans | Measuring the pore size distribution of a model granular material using a water suction method[END_REF] has accomplished extensive experimental work on the measurement of pore size distribution for the description of the acoustic properties of porous media. 6 This article concerns the specific area and not the pore size distribution but the same experimental principle is used. An improved version of Evans' setup giving more accurate and reliable results has been developed. The experimental principle (Fig. 1) consists of applying a pressure drop to a sample initially saturated with water and measuring the volume of water extracted by capillary action. The pressure drop is obtained by lowering the burette from the zero position, taken at the top of the sample. This configuration insures a constancy of the head level and therefore, of the applied pressure variations. The excess water drained out of the sample flows from the top of the burette into a measuring cylinder.

Preliminary experimental work consists of washing each sample and evaluating the volume, the porosity and the dry density. During the experiment, all the necessary information is obtained from the measurement of the head level and of the cumulative volume of water extracted. The remaining water in the sample after an experiment was less than 15% for the worst case (Coustone samples).

IV. RESULTS

The results for three grades of loose glass beads and three acoustic samples (Coustone) made of flint particles consolidated with a binder are presented. Basic information on each sample is given in Table I. For the Coustone samples, N corresponds to the nominal concentration of binder {3% in mass) while t and 2 correspond respectively to half and twice the nominal concentration. The extrar;:tion curve, i.e., the volume of water removed normalized by the total volume extracted from the sample, is plotted in Fig. 2 as a function of the applied pressure for the three grades of glass beads and for Coustone N (similar results as for the grades N have been obtained for the grades t and 2). Applying Eq. ( 6), the total surface area is deduced for each sample from the area of the left half space defined by the curve V(P) and the axes (P,O), (P,l) and (O,V). The constants chosen for the calculation of the integral are a wa = 0.07275 N/m and ()= oo. Since the pressures involved are small when using water, it is assumed that the constants do not to vary during the experiment. The applied pressure is obtained from the head level difference H with respect to the zero position: P=pg H, (7) where p is the fluid density and g the gravitational constant.

•••••••••••••;;~~•~
The specific area S s is then obtained by dividing the total surface area S by the mass of the sample. Theoretical values can be calculated in the case of glass beads. For a collection of identical spherical particles of radius r,, of density Ps CPs= 2480 kg/m 3 for glass) and without surface roughness, the specific area is given by 3 SsTh=--.

rs Ps (8) The results are summarized in Table II. The thermal length appearing in this table is explained and discussed in Sec. VI.

V. DISCUSSION OF RESULTS

The comparison of specific area measured by water extraction with theoretical values is possible only for glass beads. A reasonable agreement of 6.5% and 9.9% has been found between the measured values and the calculated ones for the two smallest grades (Table II). The larger discrepancy of 17.5% obtained for the largest beads is likely due to the gravitational force acting on the fluid. This effect adds to the capillary force and increases the actual force applied to the water. An order of magnitude of the ratio q of the gravitational to capillary forces between the top and the bottom of the sample can be calculated by 11Pg pger q = !:J.p = 2 (J' wa COS (}' (9) where e is the thickness of the sample. For a sample of 1.5 em thickness made of 1.64 mm diameter beads, the average pore size is about 400 JLm and q is of the order of 20% while it is less for smaller beads. According to Eq. ( 9), the method can be optimized by reducing the thickness of the sample or by increasing the surface tension while keeping the contact angle close to 0°. For the grades 136-010 and 136-021, the experiments were initially performed in samples of 3 em thickness. The results obtained for the specific area were respectively 2.42 and 1.94 m 2 /kg and the deviation from the theoretical values respectively 11.0% and 31%. This confirms the influence of gravity and shows that the precision can be improved for coarse materials by reducing the sample thickness.

The limitations of the method are the dependence of the contact angle on the material, the effect of gravity and the effect of the remaining water trapped in the cracks at the end of the experiment. In addition, no chemical reaction must occur between water and solid. This method is limited to materials with small pores where the capillary law holds and where capillary forces are much greater than gravitational forces. The effect of gravity can be accounted for through a correction factor that should be developed in order to enlarge the applicability range of the method to coarser materials.

Van Brakel et a/. [START_REF] Van Brakel | [END_REF] have argued that measurements of surface areas by mercury porosimetry will never be reliable because of the structure hysteresis observed in the curves V(P) depending on whether the mercury is injected or extracted. Indeed, the assumption made on the reversibility of the work during an intrusion or an extraction process may break down because of the interconnectivity of the pores and because of the existence of ink bottles and dead ends in the material. In addition, results can be affected by the sample deformations due to the high pressures required to inject mercury in powders. However, the assumption of reversibility can be checked and if it breaks down, an average value and limits can be worked out and compared to BET measurements for finely dispersed media. We did not investigate the structure hysteresis at the moment because this would have required a more complicated experimental setup. In our experiments, the same assumptions are made as by Rootare and Prenslow but the effect of sample deformation is thought to be of less importance because the pressures involved are less when using water to fill bigger pores. Considering the reasonably good agreement between the measured and the calculated specific areas, the method can be considered as reliable and more adapted to the coarse materials encountered in acoustics than the BET used for very finely dispersed materials with higher specific areas.

VI. SPECIFIC AREA AND ACOUSTICS OF POROUS MEDIA

The specific area is widely used in chemistry and in other fields of physics. It is an important parameter for the physical study of wave propagation in porous media.

For some new acoustic materials such as those fabricated from loose solid particles consolidated with a binder or such as porous foams impregnated with a chemicals, 16 the specific area can be used to evaluate the average thickness of the layer h, provided it is thin with respect to the pore size. Assuming a constant film spreading over the solid surface, a simple calculation yields the following expression that can be used in a first approximation:

hS 5 =Y 1 , PI (10)
where y 1 and p 1 are respectively the mass content and the density of the added matter. In the case of Coustone, the density of binder being 1130 kg/m 3 and the binder content being respectively 0.03, 0.015 and 0.06 for the grades N, t and 2, the application of Eq. ( 10) gives the values 29, II and 41 JLm, respectively, for the average thickness of binder. In fact, the layer thickness is not constant for high concentrations of binder and only a fraction spreads over the solid surface while the rest accumulates between the grains. 17 For the same reason, this simple evaluation would not hold for materials with a strong surface roughness or with a fractal porosity. Nevertheless, Eq. ( 10) can be considered as a good starting point for a more refined physical analysis.

Important contributions to the physical study of elastic wave propagation in porous media have been brought out by Johnson et al. 1 and by Champoux and Allard 2 who have introduced the concepts of viscous characteristic length A and of thermal characteristic length A 1

• These two parameters are necessary for the acoustic description in the full frequency range without making assumptions on the pore shape.

The viscous characteristic length A is given by the volume-to-surface ratio weighted by the fluid velocity field in the pore: 1 (II) At this time, the experimental determination of A is possible only under particular conditions or for certain materials. The intrinsic definition (II) does not allow any direct measurement for every materials, but an estimation via the thermal length is possible knowing that A is generally between 1.5 and 3 times less than A 1

• A direct evaluation of A would require the additional knowledge of the velocity field in the pores.

For a material of porosity </J, the thermal length is given by2

fdV </J Al=2 fdS =2-S -. S Ps (12)
Lemarinier et al. 18 have shown the connection between the specific area, the thermal length and the dynamic compressibility of air in high porosity plastic foams. They have evaluated A 1 from measurements of specific area by the BET method. The BET method is an original method based on the adsorption of the molecules of a gas by the particles of a powder. For a given quantity of powder (generally a few grams), the precision depends on the gas used and on the grain size. The quantity of gas adsorbed being proportional to the surface, the problem consists in determining the number of moles contained in a single layer and the total amount adsorbed on the solid surface. The model generally used to evaluate this amount is the one released by Brunauer, Emmet and Teller. 13 This method is widely used for the characterization of powders in chemistry and pharmaceutical industry where the order of magnitude of the measured specific areas is greater than a few dozens of m 2 /g. The method proposed in this article is an alternative and complementary method to the BET technique that should apply to many acoustic materials whose specific areas, of a few m 2 /kg, are much lower than those of powders. The overall thermal lengths calculated for our samples from the experimental data are displayed in Table II.

On the basis of the theory of Johnson et a/., 1 Lafarge et a/. [START_REF] Lafarge | [END_REF] have extended the model of Champoux and Allard 2 and proposed a more refined analysis of thermal exchanges in air-filled porous structures introducing the concept of a dynamic thermal permeability k 1 ( w) as a function of the angular frequency w. The low frequency limit of k 1 ( w) is a geometrical constant called the thermal permeability kh .

This new parameter has been identified as the inverse of the trapping constant r which is a parameter defined in problems of diffusion-controlled reactions in disordered media. 19 The trapping constant has been measured recently by an acoustical method.Z 0 Lafarge et a!. have shown that the knowledge of the trapping constant (or the thermal permeability kh) and of the specific area (or the thermal length A 1 ) is enough to predict, in the whole frequency range, the dynamic compressibility of air and therefore the acoustical properties of air-filled porous media.

VII. CONCLUSION

A method for measuring the specific area of acoustic materials has been proposed and tested for several grades of glass beads and materials used in civil engineering. This method is modelless and is based on the integration of the work required to extract water from a porous sample by capillary action. A good accuracy (6.5%-9.9%) has been obtained for the smallest particle size (0.68-1.1 mm) while the method fails for higher grades (17.5% for 1.64 mm beads) because of the effect of gravity acting on water. Nevertheless, the precision can be improved for coarse materials and the method can be considered as reliable and well adapted to acoustic materials.

The surface area is an important parameter for the study of porous media from the physical standpoint simply because the fluid/solid interface represents the frontier where the physicochemical reactions take place and where the energy exchanges occur. It is needed for the microscopic characterization of modem acoustic materials and for the description of the acoustic properties of air-filled materials in the whole frequency range.

/ 1 -

 1 FIG. I.Experimental principle for measuring the surlace area of porous materials.

FIG. 2 .

 2 FIG.2. Normalized volume of water extracted of the porous sample as a function of the applied pressure.

TABLE I .

 I Samples characteristics.

	0.9
		Particle size		Dry density Volume Thickness
	Grade	(mm)	Porosity	(kg/m 3 )	(ml')	(em)
				362)X2480 200	=1.5
	Coustone N Coustone t	=I =I	0.337 0.365	1427 1385	182 182	1.2 1.2
	Coustone 2	=I	0.315	1439	182	1.2

Beads 136-014 0.68:!:0.12 0.385 (1-0.385)X2480 300 =2.3 Beads 136-010 1.11 :!:0.15 0.365 (1-0.365)X2480 200 =1.5 Beads 136-021 1.64:!:0.15 0.362 (!-0.

TABLE II .

 II Results.

		Specific area Theoretical value Difference Thermal length
	Grade	(m 2 /kg)	(m 2 /kg)	(%)	(~J-m)
	Beads 136-014	3.79	3.56	6.5	133
	Beads 136-010	1.96	2.18	9.9	236
	Beads 136-021	1.74	1.48	17.5	263
	Coustone N	0.92			511
	Coustone ±	1.2			431
	Coustone 2	1.3			340

ACKNOWLEDGMENTS

The authors thoroughly thank D. C. Hothersall, J. F. Allard, A. Cummings and A. A. Javadi for their comments on this article and for their enlightening remarks. They are also grateful to E. Chadwick of Sound Absorption LTD., Bacup, Lancashire, UK, for providing acoustic materials and for stimulating discussions. This work has been made under the British Engineering and Physical Sciences Research Council (EPSRC) Grant No. C36.


Daniel Balouek-Thomert

Arya K Bhattacharya

Eddy Caron

Karunakar Gadireddy

Laurent Lefevre

Parallel Differential Evolution approach for Cloud workflow placements under simultaneous optimization of multiple objectives

Keywords: Parallel computing, workflow placement, cloud computing, energy, makespan, middleware, Differential Evolution, multi-objective evolutionary algorithm, Pareto front convergence

The recent rapid expansion of Cloud computing facilities triggers an attendant challenge to facility providers and users for methods for optimal placement of workflows on distributed resources, under the oftencontradictory impulses of minimizing makespan, energy consumption, and other metrics. Evolutionary Optimization techniques that from theoretical principles are guaranteed to provide globally optimum solutions, are among the most powerful tools to achieve such optimal placements. Multi-Objective Evolutionary algorithms by design work upon contradictory objectives, gradually evolving across generations towards a converged Pareto front representing optimal decision variables -in this case the mapping of tasks to resources on clusters. However the computation time taken by such algorithms for convergence makes them prohibitive for real time placements because of the adverse impact on makespan. This work describes parallelization, on the same cluster, of a Multi-Objective Differential Evolution method (NSDE-2) for optimization of workflow placement, and the attendant speedups that bring the implicit accuracy of the method into the realm of practical utility. Experimental validation is performed on a real-life testbed using diverse Cloud traces. The solutions under different scheduling policies demonstrate significant reduction in energy consumption with some improvement in makespan.

I. INTRODUCTION

Many of the IT and Analytics services that organizations utilize currently depend on large computing infrastructures that are hosted either locally or at remote data centers [START_REF] Foster | Computational grids[END_REF]. A popular business model for renting out resources of a data center is provided by Cloud Computing, which enables customers to allocate computing, storage and network capacity over the Internet and pay by the hour of use. In recent years, concerns about energy consumption are increasingly becoming common as Clouds often consume a large amount of electricity to power and cool computing resources in their datacentres [START_REF] Dongarra | The International Exascale Software Project roadmap[END_REF]. This situation is partially caused by an overprovisioning to ensure service delivery at peak hours, leading to underutilized resources at other times [START_REF] Orgerie | A survey on techniques for improving the energy efficiency of large-scale distributed systems[END_REF]. Efficient allocation of tasks to resources can improve consolidation on a minimum number of nodes, while transiting remaining unused nodes to low-power modes or shutdown [START_REF] Beloglazov | Managing overloaded hosts for dynamic onsolidation of virtual machines in cloud data centers under quality of service constraints[END_REF], [START_REF] Orgerie | When Clouds become Green: the Green Open Cloud Architecture[END_REF]. The implementation of the task-to-resource allocation policy consists in picking in real time at Cloud provider's end the best combination of resources, in order to fit the customer's needs at lowest cost, risk and energy consumption. Server allocation policies usually involve two actors: the Cloud provider who defines placement policies according to available resources while the customer submits sets of tasks to be executed. Such policies must benefit both the provider and customer in terms of the above metrics.

In previous work [START_REF] Balouek-Thomert | Energyaware server provisioning by introducing middlewarelevel dynamic green scheduling[END_REF], the present authors proposed methods for provisioning resources and distributing requests with the objective of meeting performance requirements while reducing energy consumption. GreenPerf, a hybrid metric, was introduced as a ratio of performance and power consumption for energy efficiency. The proposed solution considered willingness to perform energy savings by balancing user's and provider's preferences when scheduling the requests over the physical nodes. However, considering the contradictory nature of these objectives, GreenPerf could not fully explore the large domain space of possible solutions. The search and computation of these solutions is a NP-Hard problem, which can be formulated as an optimization problem with multiple contrary objectives: minimizing both energy consumption and completion time. In this work, we have used Non-Dominated Sorting Differential Evolution to obtain the best Pareto front with a spectrum of solutions representing minimum energy at one end of the front and minimum makespan (completion time) at the other.

The field of multi-objective optimization, particularly techniques using evolutionary algorithms, has advanced significantly since the first attempt [START_REF] Schaffer | Multiple objective optimization with vector evaluated genetic algorithms[END_REF] using Genetic Algorithms, and is widely used today in numerous applications. Among the most noteworthy developments rank the SPEA2 algorithm by Zitzler et al [START_REF] Zitzler | SPEA-2: Improved performance of the Strength Pareto Evolutionary Algorithm[END_REF] and NSGA-II algorithm by Deb et al [START_REF] Deb | A fast and elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF]. Some aspects of the latter approach have been incorporated in the development of the Non-Dominated Sorting Differential Evolution II technique (NSDE-II), used in the current paper.

Evolutionary

Algorithms that work concurrently on a population of candidate solutions are naturally amenable to parallelization and consequent speedup, because a significant percent of the computations operate on individual candidates independent of the others. There are two broad paradigms for parallelization, the "master-slave" model [START_REF] Lampinen | Differential evolution: New naturally parallel approach for engineering design optimization[END_REF] and the "island" model [START_REF] Tasoulis | Parallel differential evolution[END_REF]. Talbi et al [START_REF] Talbi | Parallel Approaches for Multi-Objective Optimi-zation[END_REF] provides a comparative analysis of various approaches towards implementation of parallelism for Multi-Objective Evolutionary Optimization. This work focuses on workflow applications that consist of multiple components (tasks) related by precedence constraints that usually follow the data flow between them. Although this is the most common situation, precedence constraints may exist for other reasons, and be arbitrarily defined by the user. We intend to integrate NSDE-2 as a Multi-Objective Optimization engine within a large scale infrastructure. NSDE-2 would be accessible as a remote service that accepts a workflow as an input and computes a set of placement solutions that minimizes energy consumption and makespan as an output. This output is to be placed and executed on the infrastructure using the DIET Middleware. In this framework the time spent on NSDE-2 optimization contributes to the makespan, hence this work addresses speedup of NSDE-2 through parallelization. The current version uses only energy and makespan as the objectives for concurrent minimization, the intention is to gradually integrate more independent objectives into the optimization process.

This paper introduces several contributions: (i) an evolutionary approach to workflow placement (ii) a choice of solutions to the user based on his priorities, ranging from best-energy to best-makespan, and intermediates, (iii) an experimental protocol using a real life testbed and (iv) parallel launching of evolutionary optimization process on the same Cloud infrastructure targeted for workflow placement.

The remainder of this paper is structured as follows.

Section II presents Differential Evolution, the DIET toolkit and a short summary of related works from the literature. In Section III, we introduce the Non-Dominated Sorting Differential Evolution II algorithm. In Section IV we present the problem formulation. In Section V, we propose a generic and customizable infrastructure for workload placement on a large scale infrastructure. In Section VI we evaluate the quality of workflow placement using our approach. Finally we draw Conclusions and discuss further lines of development.

II. BACKGROUND AND RELATED WORK

In this section we provide an overview of evolutionary optimization. We then introduce the DIET middleware and the features used in this paper. Finally, we present a short summary of related work on workload placement using multiobjective evolutionary optimization.

A. Differential Evolution

The developments in Multi-Objective Evolutionary Algorithms referred in Section I have been along the track of Genetic Algorithm (GA) [START_REF] Holland | Adaptation in natural and artificial systems: An introductory analysis with applications to biology, control, and artificial intelligence[END_REF], the baseline Evolutionary optimization approach, applied to the direct multi-objective paradigm. At the basic algorithm level, Differential Evolution (DE) was formulated as an alternate approach to GA by Storn and Price [START_REF] Storn | Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces[END_REF]. The present authors have applied both GA and DE in a few complex industrial processes [START_REF] Bhattacharya | Optimization of Continuous Casting Mould Oscillation Parameters in Steel Manufacturing Process using Genetic Algorithms[END_REF][START_REF] Bhattacharya | Estimation of operating Blast Furnace reactor invisible interior surface using Differential Evolution[END_REF][START_REF] Bhattacharya | Optimization of oscillation parameters in continuous casting process of steel manufacturing: Genetic Algorithms versus Differential Evolution[END_REF]; the latter work also provides a comparison in computational efficiency for that industrial process between GA and DE demonstrating that DE comes out favourably. Due to these developments the authors decided to use their version of DE as the baseline algorithm for the current multi-objective problem.

Evolutionary Algorithms have also been successfully parallelized on Cloud frameworks. Lee et al [START_REF] Lee | Designing a parallel evolutionary algorithm for inferring gene networks on the Cloud computing environ-ment[END_REF] implemented a parallel GA-PSO method for inferring gene networks in a Cloud computing environment using the Hadoop MapReduce programming model. Tang et al [START_REF] Deng | A parallel version of Differential Evolution based on Resilient Distributed Datasets model[END_REF] parallelized the DE algorithm using a resilient distributed datasets model, and compared consequent performance improvements relative to MapReduce on a wide range of benchmark problems.

The above examples represent parallelization of single-objective evolutionary algorithms on Cloud clusters to solve specific optimization problems, and not scheduling of actual workflows based on multiple objectives.

B. The DIET Middleware

DIET [START_REF] Caron | DIET: A scalable toolbox to build network enabled servers on the grid[END_REF] is an open-source middleware that enables a scalable execution of applications. Tasks are scheduled on distributed resources using a hierarchy of agents, as shown in Figure 1. DIET comprises several elements, including:

• Client application that uses the DIET infrastructure for remote problem solving.

• Server Daemon (SeD), which acts as a service provider exposing functionality through a standardized computational service interface. A single SeD can offer any number of computational services.

• Agents, deployed alone or in a hierarchy, facilitate service location and invocation interactions between clients and SEDs. Collectively, a hierarchy of agents provides high-level and scalable services such as scheduling and data management. The head of a hierarchy is termed as Master Agent (MA) whereas the others are Local Agents (LA). Applications are given a degree of control over the scheduling subsystem using plug-in schedulers (available in each agent) that use information gathered from resources via estimation functions (filled by each SeD). When a SeD receives a user request, by default it uses a pre-defined function to populate an estimation vector with system related information. A developer can create his own performance estimation function and include it into a SeD so that when the SeD receives a user request, the custom function is called to populate an estimation vector. These estimation vectors are used by agents to locate and invoke services required to execute a user application. Typically, a client request is made to a MA, which in turn broadcasts it to its agent hierarchy.

Another feature used in this work is DIET workload management capabilities. The DIET engine can handle workflow by assigning tasks to SeDs using one DIET service call. This assignment is made internally and dynamically by the MA, which receives requests from clients containing the description of a workflow. In this context, the MA determines how to schedule the workflow according to:

• Precedence constraints between tasks

• Scheduling policies/current plug-in schedulers

• Service performance properties

• Available resources on the infrastructure. This work uses the design of a new DIET plugin scheduler to express information about servers' performance and power consumption, which is then taken into account when servers are provisioned to applications. Estimation vectors are used to determine the suitability of different SEDs while considering energy efficiency for executing the workflow and performance when executing the optimization engine service.

C. Related Work

Several approaches using multi-objective optimization to manage workload placement are present in the literature [START_REF] Talukder | Multiobjective differential evolution for scheduling workflow applications on global grids[END_REF], [START_REF] Tsai | Optimized task scheduling and resource allocation on cloud computing environment using improved differential evolution algorithm[END_REF]. Objectives refer to load balancing [START_REF] Abdulmohson | Energy efficient load balancing of virtual machines in cloud environments[END_REF], load prediction or platform reconfiguration [START_REF] Legillon | A multi-objective evolutionary algorithm for cloud platform reconfiguration[END_REF], among others. A Pliant logic approach is used in [START_REF] Benyi | Energy-aware VM scheduling in IaaS clouds using pliant logic[END_REF] to improve energy efficiency in simulation based experiments. The authors conclude with the need to find trade-offs between energy consumption and execution time for optimization. Although most of the above works deal with workflow scheduling on Clouds using Multi-Objective Evolutionary algorithms, they have not explored the parallelism potential of the Cloud infrastructure in the scheduling process itself.

One of the first developments in that direction is seen in [START_REF] Mezmaz | A parallel bi-objective hybrid metaheuristic for energy-aware scheduling for cloud computing systems[END_REF], where a Genetic Algorithm is used for optimization and Dynamic Voltage Scaling to minimize energy consumption. A comprehensive review of the state of the field is presented in [START_REF] Tsai | Metaheuristic Scheduling for Cloud: A Survey[END_REF]; work on parallelism of Differential Evolution algorithms in this context is yet to be reported.

Moreover, existing work [START_REF] Kim | Power aware scheduling of bag-oftasks applications with deadline constraints on DVS-enabled clusters[END_REF] commonly assume that nodes from a homogeneous cluster are identical in power consumption and performance, which is not always true in practice. Causes of variation include external environmental factors, such as temperature and node location in a rack, aging of components due to use and leakage power that varies over time [START_REF] Diouri | Your cluster is not power homogeneous: Take care when designing green schedulers![END_REF]. We conclude that scheduling decisions based on performance and energy consumption values of the machines should be evaluated and dynamically adjusted using live monitoring.

From a resource management perspective, Grids and Clouds use meta schedulers to schedule jobs across multiple sites and local resource managers that control computational resources at a site level. Users commonly submit batch jobs to request resources over a period [START_REF] Capit | A batch scheduler with high level components[END_REF]. Cloud aggregators such as RightScale provide application-specific Cloud management and load balancing. At an application level, distributed OS such as [START_REF] Wentzlaff | An operating system for multicore and clouds: Mechanisms and implementation[END_REF] offer programming models that allow OS services to scale to match demand. Most of these systems, however, neither take energy efficiency into account nor offer means for users to specify how they want to schedule their applications while exploring tradeoffs between energy efficiency and performance [START_REF] Tu | A power-aware cloud architecture with smart metering[END_REF].

III. NON-DOMINATED SORTING DIFFERENTIAL EVOLUTION II (NSDE-II)

Differential Evolution (DE) belongs to the broad class of evolutionary optimization techniques that developed as distinctive variants of classical Genetic Algorithms (GA). DE was selected as the evolutionary method of choice on the basis of the authors' prior studies on the relative efficiency and merits of this against GA, as reported in [START_REF] Bhattacharya | Optimization of oscillation parameters in continuous casting process of steel manufacturing: Genetic Algorithms versus Differential Evolution[END_REF]. This section presents the concept of differential evolution for a single objective and the key aspects to adapt it to Multi-Objective Differential Evolution.

A. Baseline Differential Evolution

Formally, if the dimensionality of the solution space is denoted as D and the number of candidate solutions is N, then the elements of the i th vector of the solution Xi,G at generation G may be denoted as ,

....., ) i G i G i G i G D i G X x x x x = for all i N …  (1) 1, , 2, , 3, , , , , ( , , ...
The DE process fundamentally generates new solutions from the current candidate set by adding the weighted difference between two randomly selected candidate solution vectors to a third to generate a "mutant" vector, and then creating a crossover between an existing vector and the mutant, that is called the "trial" vector. The latter is allowed to replace the existing vector only if it is found to be more "fit" -the complexity of this "fitness determination" exercise depending entirely upon the nature of the problem under consideration.

If Vi,G represents the mutant vector, then according to the baseline DE process called DE/rand/1 [START_REF] Storn | Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces[END_REF] ,

1,

i G r G r G r G V X F X X = + × - … (2) 2, 3, ( )
where r1, r2 and r3 are random integers less than N, different from each other and from 'i', and F usually lies between 0.5 and 1. There are many variations of this baseline process where two instead of one difference terms are sometimes considered, the best solution in a population is taken into account, etc.; descriptions of alternative schemes may be seen in [START_REF] Caron | Deployment of a hierarchical middleware[END_REF], among others.

Crossover is performed between the 'mutant' vector Vi,G and the target vector Xi,G to generate a 'trial' vector Zi,G according to , , j , , , , if rand (0,1)

otherwise j i G j i G j i G v Cr z x ≤   =    … (3)
where zj,i,G is the element j of the trial vector Zi,G, randj (0, 1) denotes a random number between 0 & 1 applied to the element j, Cr is the crossover threshold usually set between 0.4 and 1. At the final selection step the choice for candidate 'i' in the next generation is made between Zi,G and Xi,G on the basis of higher fitness by direct one-to-one comparison.

The present work generates the mutant vector according to the alternate scheme (proposed in [START_REF] Storn | Differential evolution: a simple and efficient heuristic for global optimization over continuous spaces[END_REF] and also used by current authors in [START_REF] Bhattacharya | Optimization of Continuous Casting Mould Oscillation Parameters in Steel Manufacturing Process using Genetic Algorithms[END_REF][START_REF] Bhattacharya | Estimation of operating Blast Furnace reactor invisible interior surface using Differential Evolution[END_REF][START_REF] Bhattacharya | Optimization of oscillation parameters in continuous casting process of steel manufacturing: Genetic Algorithms versus Differential Evolution[END_REF] where it is found to work better than other DE variants) ,

i G r G best G r G r G r G X X R X X F X X = + × - + × - ( 1, , 1, 2, 3, ( ) ( )
) 4
where R is set at 0.5 and F varies randomly between -2 and +2 across generations (and are same for all 'i' within a generation). The crossover probability Cr in eq. ( 3) is set at 0.7.

B. Multi-Objective Differential Evolution NSDE-II

Compared to single-objective DE discussed in Sec. 3A, the mechanisms of selection to a new generation in multi-objective DE are radically different. The basis for this difference lies in the fact that one cannot uniquely order the candidate solutions based on their "fitness" when the number of axes for ordering is more than one. Thus when comparing two candidates Zi,G and Xi,G, the former may be better for the first objective and the reverse for the second. This problem on transiting from single-to multipleobjectives is equally relevant for any evolutionary algorithm and not for DE alone. This work has adopted the basic multipleobjective handling techniques of NSGA-II [START_REF] Deb | A fast and elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF] while replacing the baseline GA operations to those of the DE variant outlined in (3) and ( 4) for generation of a trial vector. Hence this is named as NSDE-II. It may be noted that in DE "elitism" and hence selection of the better performing solutions executes a more integral and critical role in the evolutionary process compared to typical GA where elitism is not mandatory.

In a problem with K objectives FFk, k 1…K, a  candidate solution vector Xp is said to dominate another solution Xq if

{ } ( ) ( ), 1,..., k p k q FF X FF X k K ≥ ∀ ∈ … (5)
and for at least one k, ( ) ( )

k p k q

FF X FF X >

; where p, q {1, …, N}, N is population size; and in  turn Xq is said to be dominated by Xp. Now it is apparent that for a population of candidate solutions and with multiple objectives, there will be either one of three types of relations between any pair of candidate solutions -either one dominates the other according to (5), or one is dominated by the other, or neither dominates or is dominated by the other.

The NSDE-II algorithm exploits these three basic types of relationships to select candidates from the population pool into the next generation, by stratifying the pool into fronts based on descending degree of non-dominance. At a secondary level it also exploits the degree of diversity among solutions, giving preference to candidates with higher diversity, using the Crowding Distance Algorithm [START_REF] Deb | A fast and elitist Multiobjective Genetic Algorithm: NSGA-II[END_REF].

IV. PROBLEM FORMULATION

The aim of this work is to improve the energy efficiency of a set of machines while concurrently reducing completion time of a given set of jobs, through optimized workload placement. A server (computing node) is modeled with three resources: CPU, DISK and NETWORK and runs processes which consume these resources. Each process is to be assigned to a single machine, and cannot be moved from one machine to another.

A. Decision Parameters

Any optimization problem will have design parameters whose best possible values from the viewpoint of the objectives are sought to be attained in the optimization process. The optimization task here is to map a given set of jobs in a certain sequence onto the available resources.

Suppose there are m number of resources and n tasks. Then, for any resource j, j [1,..m], all   possible permutations of subsets of all sizes of the set of tasks of size n, constitute the total solution space. If we call the size of this solution space as Sj , then

0 ! ( , ) ( , ) ( )! n j k n S P n k where P n k n k = = = - ∑ … ( 6
) Since Sj is independent of j, we may write it simply as S. It then follows that the size of the total solution space is m S .

The following information is assumed to be known for each server s, and task i, at any time: The knowledge of these variables enables the scheduler to extract the energy consumption and makespan related to the completion of tasks over the available resources.

fs

B. Objective Functions

We have two objective functions:

1) Minimize Makespan (i.e. time taken for completion of all tasks in the workflow)

If Tis is the completion time of the task i launched on server s, then

i i i i is s s s s nf nbw nbr nnet T f dw dr net = + + + … (7a)
If the task i is launched on a server other than s, then

0 is T = … (7b)
The completion time of the workflow is expressed as

1 1 m n mn is s i T T = = = ∑∑ … (8a)
Eq. (8a) is the total time taken to execute all tasks on all servers. However, by definition, makespan is the time when the last server completes its task in the workflow. Hence, (8a) is modified to yield:

(1,.., ) 1 Makespan n mn is s m i T Max T ∈ = = ∑ … (8b)
2) Minimize Energy consumption (i.e. total energy consumed in a workflow) If Cis is the energy consumption of the task i per unit time running on server s, then energy consumption on server s required for the workflow may be expressed as

1 n s is is i W C T = = ∑ … (9)
and the total energy consumed in the workflow is

1 1 1 m m n mn s is is s s i W W C T = = = = = ∑ ∑∑ … (10)
In most cases, faster machines (low Tis) will have higher energy consumption (high Cis), implying that objectives Tmn and Wmn are contradictory -forming the basis for multiobjective optimization.

V. FRAMEWORK FOR WORKLOAD PLACEMENT

To cope with real conditions such as the increasing scale of modern data centers, as well as the workload dynamics and application characteristics that are specific to the Cloud Computing paradigm, DIET allows users to study large-scale scenarios that involve thousands of nodes, each executing a specific workload that evolves during the computation.

The aim of the current framework is twofold: (i) to relieve researchers of the burden of dealing with deployment, resource selection and workload fluctuations when they evaluate new optimization engines and (ii) to offer the possibility to compare them. To perform placement decisions, users encapsulate their optimisation engine in a program, and express the workflow along with the precedence between tasks. The program typically leverages DIET API that allows end users to create and execute remote services 1 . The Master Agent keeps a description of the physical resources, dynamically updated by the nodes hosting the services. Finally, the workload execution is orchestrated by the DIET workflow engine that internally relies on a customizable scheduler (cf. Section II) to assign the resources during the entire execution. We chose to base our framework on DIET since (i) the latter's relevance in terms of performance and validity has already been demonstrated [START_REF] Das | Differential Evolution: a survey of the state-of-the-art[END_REF] and (ii) because it has been recently extended to integrate energy-efficient decision capabilities [START_REF] Balouek-Thomert | Energyaware server provisioning by introducing middlewarelevel dynamic green scheduling[END_REF].

The workflow execution is performed in 3 phases: (i) service discovery, (ii) computation of mapping solutions and (iii) workload placement. The service discovery phase corresponds to the search of an optimization engine within the infrastructure by a given client. As multiple engines can be instantiated on the platform, the user can submit its workload to different engines and compare the cost of generated solutions. The computation of mapping solutions is performed by at least one server with a platform performance description provided by the Master Agent. This description is either based on historical data (past computations) or userdefined benchmarks. Finally, the workload placement is performed and results are returned to the client based on the platform available metrics and monitoring resolution. Mapping solutions are defined as a collection of JSON objects. Each solution contains the mapping between a SeD and a task and an associated cost in terms of workflow completion time and energy consumption.

Two kinds of experiments have been performed to validate this approach. The objective of the first one is to evaluate the computation phase of the engine (i.e., the step where the optimization engine generates a spectrum of solutions) while the second is a comparison of algorithms to evaluate the 1 http://graal.ens-lyon.fr/diet/UsersManualDIET2.9/ concrete gain of NSDE-2 compared to an online placement of workload.

VI. EVALUATION OF WORKFLOW PLACEMENT

In this section we first briefly describe the evaluation testbed. Then we look at the approach and consequences of parallelization of NSDE-2 algorithm. Finally we evaluate the quality of workflow placement using this multi-objective evolutionary optimization approach compared to extant methods.

A. Evaluation Testbed

Experiments used resources from GRID'5000 [START_REF] Cappello | Grid'5000: a large scale, reconfigurable, controllable and monitorable Grid platform[END_REF], a testbed designed to support experimentdriven research in parallel and distributed systems.

Located in France, GRID'5000 comprises 29 heterogeneous clusters, with 1,100 nodes, 7,400 CPU cores with various generations of technology spanning 10 physical sites interconnected by a dedicated 10 Gbps backbone network.

The power measurement in the studied clusters is performed with an energy-sensing infrastructure composed of external wattmeters produced by the SME Omegawatt. This energysensing infrastructure, also used in previous work [START_REF] Balouek-Thomert | Energyaware server provisioning by introducing middlewarelevel dynamic green scheduling[END_REF], [START_REF] Assuncao | An analysis of power consumption logs from a monitored grid site[END_REF], collects at every second the power consumption in averaged watts of each monitored node [START_REF] Assuncao | The green grid'5000: Instrumenting and using a grid with energy sensors[END_REF]. A node's consumption is determined by averaging past consumption over more than 6,000 measurements, whereas its performance is given by the number of FLOPS achieved when using a single CPU cores to execute benchmarks.

We deploy the DIET middleware on 113 physical nodes as follows: 111 dedicated nodes for SeD's, 1 dedicated node for the Master Agent and 1 dedicated node for the Client. The machines are picked among six different clusters as presented in Table I. Detailed description of trace files and their usage in experiments are provided in a related work [START_REF] Balouek | Minimizing Energy and Makespan concurrently in Cloud Computing workloads using Multi-Objective Differential Evolution[END_REF].

B. Parallelization Investigations

We first investigate parallelization of the NSDE-2 algorithm on a handy 8-core Intel laptop with chipset i7-4710HQ@2.5Ghz. These are offline simulations using data for a set of 500 tasks to be placed on 85 servers, where task and server data have been extracted from the GRID'5000 testbed. A population size of 200 is considered for all simulations as well as online optimization executions.

The master-slave approach is followed in parallelization, using the Open MP Library. In the NSDE program the functions not amenable to parallelization include the selection operations using non-dominated sorting and crowding distance algorithms that take up about 3.3% of runtime, and some I/O operations taking approximately another 1%. It follows from Amdahl's law that the Theoretical Maximum Speedup factor is approximately 22.

Table 2 shows results obtained using the sequential NSDE-2 program, parallelized NSDE-2 program running on a single core, and on 4 cores. The data shows computation times and the average values over all candidates for the two objectives, energy and makespan, for a workflow of 500 tasks on 85 servers. The relevance of the average values in these multi-objective simulations is purely to check if the sequential and 1-core-parallel solutions match exactly, which they are observed to do. This, first and foremost, demonstrates the correctness of parallelization. Second, it shows that the speedup factor on 4 cores is 2.36.

An interesting observation from Table 2 is that the speed of evolution of candidates across generations varies between the parallel solutions and the sequential, reflected in different numerical values of the objectives at the same generation levels. It may be difficult to pinpoint the reasons for this; the evolutionary algorithm being a stochastic process is likely to behave differently when executed concurrently on different numbers of nodes, and these differences are likely to amplify over generations.

Figure 2 shows the parallelization speedup factor when running selected sets of 100 and 1000 tasks on 85 servers, real time on the GRID'5000 testbed on a Stremi node (see Table 1). It may be noted that in the NSDE solution framework, each task is effectively a decision (design) variable, and obtaining optimized solution with 1000 variables is itself a challenging task. In fact, this number has been extended to 5000 decision variables on 85 servers launched in parallel on 24 nodes, though comparative sequential runs could not be obtained due to runtime constraints. Fig. 2 shows that as the size of the workflow increases, the parallelization speedup factor gradually approaches its maximum limit.

B. Quality of NSDE-2 optimized workflow placements Figure 3 plots the evolution of NSDE-2 solutions from an initial level of 100 generations up to 10000 generations, with minimization of energy consumption and makespan as the objectives on the two axes. Each dot represents a candidate solution. At any selected generation, at one end of the solution front we have the best energy solution, and at the other end, the best makespan solution. We can observe that the quality of solutions improves as the number of generations increases, and the "cloud" of solutions gradually transforms into a Pareto front. The computation time increases linearly with the number of generations. We choose to retrieve the solution at 3000 generations, after which the improvement in solution becomes less significant.

It may be noted that if jobs are submitted on the cloud for execution after prior reservation, it can be valuable to drive the NSDE-2 to its full potential to obtain the best optimal solution placement. Next we compare the distribution of tasks among nodes on GRID'5000 attained under three different policies, namely NSDE-2 Best Energy, NSDE-2 Best Performance and FIRST FIT. NSDE-2 Best Energy and NSDE-2 Best Performance correspond, respectively, to the candidate solutions on the same Pareto front with the smallest energy consumption and the smallest makespan. These solutions establish the bounds of the Pareto Front. The FIRST FIT policy selects the first available server in an ordered list according to the GreenPerf metric as a nonweighted average ratio between makespan and energy consumption for the said type of task.

For any of considered cases there exists a proper balance between short and long tasks within the dataset. A server is restricted to the execution of, at most, one task at a given time. Considering that the scheduler does not have specific information on the nodes and does not make assumptions about the hardware, the dynamic information is gathered by computing a sample of each type of task on the servers prior to initiation of the evolutionary optimization process.

Figures 4567show the results of these studies on 100, 200, 500 and 1000 jobs launched on 111 servers. The x-axis presents the different policies used to execute the workflow; the y-left-axis shows the total energy consumption of the solution and the y-right-axis shows the makespan value.

We observe that as the complexity and size of the optimization space increases (cf. Sec IVa where the size varies as m n ), the multi-objective evolutionary algorithm provides better solutions compared to the single-metric based ranking approach of FIRST FIT. This is also reflected in Table 3. It may be seen that energy consumption improves in all NSDE scenarios (best-energy and best-makespan) up to 25%, while makespan improves for the best-makespan solution for the larger cases. Further, NSDE provides a spectrum of intermediate solutions to the user to select based on his preferences between energy and makespan.

When the NSDE-2 solution was run up to 10000 generations, it provided a 30% saving in energy with a 50% reduction in makespan without considering the algorithm convergence time. This can be of value in cases of jobs submitted by prior reservation, as typical of jobs involving large computation times.

CONCLUSIONS AND PERSPECTIVES

This work describes the design, implementation and evaluation of an energyefficient resource management system that builds upon DIET, an open source middleware and NSDE- II, an Evolutionary Multi-Objective Optimization engine. It investigates the nuances of parallel launching of the computation-intensive evolutionary algorithm on the baseline cloud infrastructure where the workflow is targeted for optimal placement. The implementation supports an IaaS Cloud and currently provides placement of workflows, considering non-divisible tasks with precedence constraints. Real-life experiment of this approach on the GRID'5000 testbed demonstrates its effectiveness in a dynamic environment. Results shows that our method can offer providers and decision makers an aid to make decisions when conflicting objectives are present, or when in search for realistic trade-offs for a given problem.

Investigations on parallelization of NSDE-2 show that speedup values approaching the theoretical maximum limit have been obtained on a Cloud cluster. It is observed that speed of evolution of solutions across generations fluctuates with the number of active cores. At the next stage, we will investigate multi-core integration of servers in a cluster, which is expected to yield significant energy savings.

ACKNOWLEDGMENTS

Fig. 1 .

 1 Fig. 1. An example of DIET Hierarchy

Fig. 2 .

 2 Fig. 2. Parallelization speedup factors for different job sizes and cores.

Fig. 3 .

 3 Fig. 3. Convergence towards a Pareto front across generations.

Fig. 4 .

 4 Fig. 4. Energy and Makespan comparison for 100 jobs and 111 servers.

Fig. 5 .

 5 Fig. 5. Energy and Makespan comparison for 200 jobs and 111 servers.

Fig. 6 .

 6 Fig. 6. Energy and Makespan comparison for 500 jobs and

TABLE I .

 I EXPERIMENTAL INFRASTRUCTURE

	Cluster Orion Sagittair	Node s 4 38	CPU 2x6 cores @2.30Ghz 2x1 core @2.40Ghz	Memory 32 GB 2 GB	Role SeD SeD
	e Taurus Stremi Graphite 4 10 38 17 Parasilo 1 Parasilo Parasilo 1	2x6 cores @2.30Ghz 2x12 cores @1.70Ghz 2x6 cores @2.00Ghz 2x6 cores @2.40Ghz 2x6 cores @2.40Ghz 2x6 cores @2.40Ghz	32 GB 48 GB 64 GB 128 GB 128 GB 128 GB	SeD SeD SeD SeD MA Client

TABLE 2 .

 2 SEQUENTIAL AND PARALLEL SIMULATIONS

	Paral leli-zatio n Seque ntial	Time for 3000 gens (mins) 33:52	Average value at 100 generations (Makespan in mins, Energy in kJ) Makesp an Energy 39.36 4491.4	Average value at 3000 generations Makes Energy pan 37.05 3546.6
	1-core	33:49	39.36	4491.4	37.05	3546.6
	4-core	14:22	40.46	4874.8	35.85	4061.3

TABLE 3

 3

	. MAKESPAN AND ENERGY USING NSDE-2 (figs. show % COMPARATIVE IMPROVEMENTS IN gains against FIRST FIT)
	Ca se s	No. of job s	NSDE -2 Comp uta-tion time (mins )	Make span for NSDE -2 Best Make span (%)	Make span for NSDE -2 Best Energ y	Energ y for NSDE -2 Best Make span	Energ y for NSDE -2 Best Energ y
	1	100	3.46	-82	-82	0	0
	2	200	6.0	-23	-59	6.1	17.3
	3	500	13.63	15	-19	12.5	17.1
	4	100 0	26.5	7	-14	17.6	24.3


Jong Uhn Kim

A BOUNDARY THIN OBSTACLE PROBLEM FOR A WAVE EQUATION

In this paper we shall discuss a hyperbolic variational inequality associated with the following initial-boundary value problem:

Here n is a bounded open subset of Rn, n ;:::: 2, with smooth boundary an E (0-1) (0-2) (0-3) (0-4) (0-5) C2. We assume that rl and r2 are disjoint, nonempty open subsets of an such that an = rl u r2 = rl u r2 and ar1 = ar2 E C 2 . The above functions J(x,t),uo(x),u1 (x) and 4>(z) are given, and f., denotes the outward normal derivative on an. We assume that 4> $ 0 and Uo ;:::: 4> on f 2 • We can put the above problem in the form of a variational iUP.quality: find u(x,t) which satisfies (0-2), (0-3), (0-4) and a 2 u < at 2 , w -u > +a( u, w -u) ;:::: < J, w -u > (0-6) for every wE G, for all t E (0, T).

< u,w >=in uwdx, a(u,w) =in V'u • V'wdx, G ={wE H 1 (!!): w = 0 on r1 and w 2: ¢on r2}.

If a solution u of this new problem is sufficiently smooth, one can show that u also satisfies (0-1) and (0-5) so that u is a solution of the original problem (0-1) through (0-5). The corresponding stationary variational inequality: find u E G which satisfies a(u,w-u) 2:< J,w-u >,for all wE G, (0-7)

is called the boundary thin obstacle problem (or the Signorini problem), and its mathematical theory has been well-established: see [START_REF] Friedman | Variational principles and free-boundary problems[END_REF] and references therein. However, to the author's knowledge, nothing has been known about the solution of (0-1) through (0-5). The purpose of this paper is to establish the existence of a (weak) solution of (0-1) through (0-5). This problem is an example of a hyperbolic variational inequality where a constraint is imposed on the unknown function rather than on the time derivative of the unknown function. Hyperbolic variational inequalities have been extensively investigated in the case where a constraint is imposed on the time derivative of the unknown function; see [START_REF] Brezis | Problemes Unilateraux[END_REF] and [START_REF] Lions | Quelques methodes de resolution des problemes aux limites non lineaires[END_REF] among others. On the other hand, when a constraint is imposed on the unknown function, existence theorems seem to be known mainly for special problems of space dimension one; see [START_REF] Amerio | Unilateral problems for the vibrating string equation[END_REF], [START_REF] Bamberger | New results on the vibrating string with an obstacle[END_REF], [START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF], [START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF], [START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints[END_REF], [START_REF] Schatzman | Un probleme hyperbolique du 2eme ordre avec contrainte unilaterale: La corde vibrante avec obstacle ponctuel[END_REF] and [START_REF] Schatzman | The penalty method for the vibrating string with an obstacle[END_REF]. In the case of several space dimensions, [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF] and this paper seem to be the only known works. [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF] resolved the above problem by microlocal techniques when n is a half-space; but the same method does not extend to a general domain, as is explained in [START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF].

In this paper, the basic procedure to prove the existence of a solution consists of a penalty method [START_REF] Lions | Quelques methodes de resolution des problemes aux limites non lineaires[END_REF], a multiplier technique [START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF] and compensated compactness [START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF]. Finally, we note that the questions of uniqueness and energy conservation for solutions of our problem are still open. For s=nonnegative real, H'(l1) is defined by interpolation and H-'(l1)=the dual of HJ(lt). We denote by <, > the duality pairing between H-"(!1) and H"(n), 0 ::; a ::; ~• Then, for g E L 2 (!1), hE H"(n), 0::; a::; t, we have < g, h >= fn ghdx.

As above, the bilinear form a(, ) is defined by a(g, h)= j 0 'lg • 'lhdx and

G ={wE Hf 1 (n): w <:</>on r2},
where </>(x) E C 1 (f!) and </>(x)::; 0 on r 2 .

Next we define the distance function d by d(x)=dist(x,an) and set

v 5 = { x E IT : d( x) < a}.
Let So be a positive number such that 1/So bounds the principal curvatures of an. Since n is bounded and an E C 2 , it is known that d(x) E C 2 (V 00 ); see [START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF]. We also set

Us= {x En: d(x) <a}, S0 = {x En: d(x) = S}.
Then, Ss is a C 2 -surface for 0 ::; S < So. We shall use the following lemmas later.

LEMMA 1.1. There is a function h(x) E C 1 (f!)n such that h(x)=the outward unit normal vector at

x on Ss ifd(x) = S::; 6o/2 and h(x) = 0 ifd(x) <: 6o.

PROOF: For each x E 50,0::; 6 <So, there is a unique point y(x) E an such that dist(x,y(x)) = 6 and the outward normal vector at y(x) on an is also outward normal at X on S5. Furthermore,

y(x) E C 1 (Vs 0 ).
Next we denote by v(y) the outward unit normal vector at yon an. Then, v(y)

is a C 1 function of y. Let us choose a function p(t) E C 1 (R) such that p(t) = 1 fort ::; 6o/2 and p(t) = 0 fort> ~60 . We then set h(x) = v(y(x))p(d(x)
), which achieves our goal.

The next lemma is known and stated without proof.

LEMMA 1.2. Let g E L 1 (n). Then, for 0 < 6* < 6o, { gdx = r'* { gdl7sd6 lu 5 • lo ls 5 (1-1)
holds where d!7 5 is the surface element on Ss.

Next we write a well-known consequence of the Div-Curllemma of compensated compactness.

LEMMA 1.3. Let {uk} be a sequence of functions such that

(1-2) (1-3)
where M is a positive constant independent of k and each Uk satisfies

(1-4)
where {gk} is a bounded sequence in L 2 (n X (0, T)), and, ask--+ oo,

(1-5)

Then, as k --+ oo [START_REF] Amerio | Unilateral problems for the vibrating string equation[END_REF][START_REF] Bamberger | New results on the vibrating string with an obstacle[END_REF][START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF] in 'D'(n x (0, T)).

i;;,i = 1, . .. ,n, and vk(n+t) -t,,wk(n+l) ~. Then, apply Corollary 4.3 of [START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF].

Finally, we state LEMMA 1.4. Let {uk} be a sequence of functions such that ask--+ oo,

(1-7)

(1-8)

where -1 ::; a< (3::; 1. Then, we have

Uk--+ u in C([O,T];Hr(n)) (1-9)
for any r < (3.

PROOF: A special case was proved in [START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF]. The same argument can still be employed. Without loss of generality, we may assume a# -!,(3 # -!,because we can take smaller a or (3 if necessary. possibly after a modification on E. Since the embedding Hr2(fl) C Hrt(fl) is compact for every r 1 < r 2 , we can use the Ascoli theorem to arrive at the conclusion.

Let r = ea + (1-8)(3,0 < e < 1,
2. Approximate Problem.

In this section we shall consider an approximate problem with a penalty parameter k > 0. The problem is: find a function uk(z,t) such that

11k(z,0) = Uok(z),8t11k(:t,0) = UJk(:t) (2-2)
< a .. u.,w > +a(u.,w) =< fk,w > + { (k(uk-4>)--~a,u.)wdu

(2-3) 1r 2
holds for all w E H/-1 (!'!), for almost all t E (0, T). Here du is the surface element on 8fl and g-=max(O,-g).

Our strategy consists of three steps.

(Step 1) Establish the existence of solution of (2-1), (2-2) and (2-3).

(

Step 2) Obtain a priori estimates independent of k > 0.

(Step 3) Pass k -+ oo to arrive at a solution of the original problem.

In this section, the first two steps will be carried out.

Existence of solution

Throughout this subsection, we shall suppress the subscript k.

PROPOSITION 2.1. Let u 0 ,u 1 E Clf(!'l),/ E C 1 ([0,T];L 2 (fl))
and 4> E C 1 (fi) with 4>::; 0 on r2.

Then there is a solution of (2-1), (2-2) and (2-3).

PRooF: We set u = v + tto + tu1. Then, and (2-3) is equivalent to

v(z,O) = 0,81v(z,O) = 0, < a,,v,w > +a(v,w) =< J,w > + < ~(u 0 + tut),w > +k { (v-4>)-wdu 1r 2 _.!_ f (8 1 v)wdu, k lr 2 (2-4) (2-
=:; 118tfiiL2(0) ll8ttVm IIL2(0) + ll~utiiL2(0) ll8ttVm IIL2(0) + ~k3 ll8t{(vm-4>nll~2(r 2 )-~ll8ttvmll~2(r 2 r
We consider (2-8) at t = 0 and substitute 8ttvm(O) for Wj:

We have used the fact that vm(O) = 81vm(O) = 0. Now (2-16) implies (2-23)

We then make use of (2-20) through (2-23) to derive (2-4) and (2-5) through a standard procedure.

Remark 2.2. The purpose of the term t81uk in (2•3) is to obtain approximate solutions which are regular enough to justify manipulations in using a multiplier technique in the next subsection.

The utility of that term is seen in (2• 15).

2 2 A priori estimates independent of 1;.

We still suppress the subscript k and let u(z,t) be a solution constructed i. n Proposition 2.1 for a fixed k > Q. We can substitute 8,u for w in (2•3):

for almost all t E (0, T). We now derive from (2•24)

lluiiHt (O) !> M, for all t E (O,T) rt ll(u-¢niL1(r,) $ M/Vk, for all t E [O,Tj, LT 1 18.ull~2(r,)dt $ Mk, (2•24) (2• 25) (2• 26) (2•27) (2-28)
where M denotes positive constants which are independent of k, &nd depend only on the number

L such tbt lluoiiHJ(O) + l!ut iiL2(0) + 11/IIL2(0x (O,T)) $ L .
Next we shall adapt a multiplier technique to derive some uniform estimates in a neighborhood of 81!. Let us set fl6 = {:t E fl: dist(z,80) > 6} for 0 < 6 !> 6o. Obviously, 8fl5 = S5, which was defined in Section 1. We need to observe tha.t solution u of (2•1), (2-2) and (2•3) has sufficient regularity to justify manipulations in using a multiplier technique. First of all, (2•3) implies tha.t 8uu -6u = f L 00 (0,T; H 2 (f!)). We also note that (2-1) implies possibly aft er a modification on a set of measure zero in t.

We now proceed to obtain the estimates. Using the function h(z) constructed in Lemma 1. 1, we multiply both sides of (Z-29) by (h.• V)u • and Integrate over fls x (0, T ), O < 6 $ 6of2:

{ I (8.,u)((h. V)u)dzdt io in 6 = I (81u(T ))((h • V)u(T))dz-I u,((h • V)ua)dz 1n 6 1o 6 -IT I (o,u)((h . V )8,u)dzdt = lo 1n 6
(applying the divergence theorem to the last integral) ;:>: 1T < j, W-u > dt, [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF] for every w such that w E L 00 (0, T; Hf-1 (fl)), 8 1 w E L""(O, T; L 2 (fl)), w(t) E G, for almost all t E (O,T).

IT I ( -~u)((h • V )u)
We shall justify this definition. Suppose that u is a solution according to the above definition.

By taking w = u ± ~,~ E Cg'(fl x (O,T)), in [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF], we find that (0-1) is satisfied in 7J'(n x (O,T)).

(0-3) and (0-4) are satisfied in the sense u(t) E G, for almost all t E (0, T). Now we assume further regularity on u:

(3-5) [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF] for each 6 > 0, where we set Es = {x E fi: dist(x,8r 1 ) :S: 6}.

Recall that 8r 1 = 8r 2 is an (n-2)-dimensional, C 2 -manifold by assumption. Then, (0-1) holds in a stronger sense:

1T < 8,.u,w > dt -1T < Ll.u,w > dt = 1T < j,w > dt (3-7)
holds for every wE L 2 (fl X (O,T)) such that supp w(x,t) n E, is empty, for almost all t, for some 0 > 0.

Next choose any nonnegative function ~(x) E C.\(r 2 ). Then, we can extend ~(x) so that ~(x) E C 1 (f!) and supp ~ n {E 8 U rl} is empty for some 6 > 0. Choose any nonnegative function 77(t) E CJ((O,T)) and take w = u + ~(x)77(t) in [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF]. By means of integration by parts and the divergence theorem (this procedure is justified by [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF] and (3-6)), we find 1T < 8,.u,~(x)1)(t) > dt-1T < Ll.~,~(x)77(t) > dt + 1T j 2 (~~) ~(x)77(t)dadt ;:>: 1T < J,e(x)7J(t) > dt.

  &;•a"= &tl' V'g = &xl , ... , &xn . When E is an open subset of RN, N ;::: 1, D'(E) denotes the space of distributions over E and CIJ'(E) = {w E cm(E) : supp w is a compact subset of E}. We assume that &i! E C 2 and &r1 = &r 2 E C 2 • For m=nonnegative integer, m ( a ) "' ( a ) "n H (!1) ={wE L 2 (l1): ax, . . . axn wE L 2 (!1),0::; a 1 + ... +an::; m}, H 0 (!1) =the completion of CQ"(l1) in Hm(!t), H-m(n) = the dual of HQ'(lt), Hf 1 (!1) ={wE H 1 (!1): w = 0 on r!}.

  r # -!, !• There is E c (O,T] such that meas(E) = 0 and for all t1,t2 E (O,T] \ E, t1 < t2, llu•(t2)-uk(tt)liw(n) ::; ::; M91iu•(t2)-u•(tt)ll~"(n)lfu•(t2)u•(tt)1111(n)' M9( li8tuk(t)IIH"(n)dt) ::; M9(t2tt) , It (1-10) holds for all k and 0 < 9 < 1 such that 9a + (1-9),8 ~ ~,-!,where M9 denotes positive constants independent oft1,t2 and k. Consequently, each uk belongs toC([O,Tj;Hr(fl)) for every a< r < ,B

M 2

 2 being a positive constant independent of m. Combining (2-13), (2-15) and (2-17), we derive (2-18) (2-19) where M denotes positive constants independent of m. By means of (2-11), (2-18) and (2-19), we can extract a subsequence still denoted by { Vm} such that (v. By virtue of Lemma 1.4, we have (vm-(W _, (v-4>)-in C([O, T]; L 2 (8fl.)).

( 2 •

 2 29) holds in 'l>'(O x (0, T)). By virtue of (2•1) and the assumption that f E C 1 ((0,TJ; L 2 (0)), we infer that (2-30) ln fact, (2•29) holds in L 00 (0, T; L'(fl)). Now (2•30) Implies 8 (2-31) for eacll 6 > 0. Because of the mixed bounda.ry condition on {Jfl, one cannot claim that u E

( 3 - 8 ) 2 :lr 2 lo lr 2 ( 3

 38223 It follows from[START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF] and[START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF] that[START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF][START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF] which is true for all nonnegative functions ~(x) E C.l(r2 ) and 7J(I) E C.l((O,T)). Hence we can conclude that we obtain (3-18) By virtue of (2-25) and (2-26), we can extract a subsequence still denoted by { uk} such that as k ...... oo, Uk ...... u weak+ in L 00 (0, T; Ht 1 (f!)), Since each Uk satisfies (2-29), we infer from (3-15) and (3-19) that Hence, by Lemma 1.4, we find that uk ...... u in C([O,T];H!(n)), so that < 8 1 uk(T),w(T)-uk(T) >-->< 8tu(T),w(T)-u(T) > as k ...... oo, for each w. In the meantime, we notice that ~~ {T f (8tuk)wdudtj lo lr ::;; ~( {T f (8 1 uddudt)!( {T f w 2 dudt)!, lo using (2-28)), :::;; M( fT f w 2 dudt)!j,fk--> 0, lo lr 2 as k--> oo, for each w. Next we shall show that (

IC(8)J $ M, for all 0 < 6 $ 6o/2, where M is a positive constant which is independent of k and depends only on the number L such tha t UuoiiHJ(O) + Uut iiL2(0) + II/ II L2(0x(O.T)) $ L.

We next choose any p(s) E C(R) such that for all s, 0 S lp(s) l $ 1, and set t,!>(x) = p(dist(x, 8!1)).

Then , 1{1( z) i s a continuous function and is constant on each S6, 0 < 6 $ 6 0 • Hence, it follows from (2-35) that (using Lemma 1.2) for each 0 < 6• $ 6o/2, where 0(6") i_ s a number such that 10(6")1 $ M 6", (2•39) where M is the same positive number as in (2•37). (2•38) will be used in the next section.

Convergence of approximate solutions.

We assume that uo E H6(l~), u1 E £ 2 (!1),/ E L 2 (!l x (O,T)),¢ E C 1 (fi) with¢>$ 0 on r 2 , and define a solution of (0• I) t hrough (0.5) .by < 81u(T), w(T)-u(T) > -< "" w(O) -1to > :~ 2: 0, almost everywhere in fz X (0, T). [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF][START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF] It remains to show the second part of (0-5). Let us choose a function {s(x) E C 1 (fi) such that 0 $ {s(x) $ 1, for all x E fi,{0 = 1 on fz \ Es and supp {s n {E 5 ; 2 U ft} is empty. We also choose any 'l(t) E CJ((O,T)) such that -1 $ 'l(t) $ 1 for all t. Then, we set

in (3.4). As above, we use integration by parts, the divergence theorem and (3-7) to arrive at

Here, we have the equality, because -'l(t) can be also chosen in place of 'l(t). Since u-¢ ;:: 0 almost everywhere in fz X (O,T) and b > 0 is arbitrary, (3-10) and [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF][START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Lions | Quelques methodes de resolution des problemes aux limites non lineaires[END_REF] imply

This completes the justification of the above definition of a solution.

We now proceed to present the maln result. fk --+ fin L 2 (fl X (0, T)).

(3-15) Using these, for each k, let Uk be a solution constructed in Proposition 2.1. Then, each Uk satisfies [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF][START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Lions | Quelques methodes de resolution des problemes aux limites non lineaires[END_REF][START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints[END_REF][START_REF] Schatzman | Un probleme hyperbolique du 2eme ordre avec contrainte unilaterale: La corde vibrante avec obstacle ponctuel[END_REF][START_REF] Schatzman | The penalty method for the vibrating string with an obstacle[END_REF](15)(16) for every w(x, t) such that wE L""(O,T; Hf-1 (fl)),& 1 w E L 00 (0,T; L 2 (fl)) and w(t) E G, for almost all t E (0, T). Using integration by parts and

It is convenient to introduce the notation: We then set tP6(x) = P6(dist(x,&fl)). Fix any 0 < o:::; 6of2. According to (3-28), as k -> oo. Furthermore, with the aid of (2-38), we notice that where we can estimate Ot(6) by using the number Min (2-37), and (3-13), [START_REF] Brezis | Problemes Unilateraux[END_REF][START_REF] Dacorogna | Weak continuity and weak lower semicontinuity of nonlinear functionals[END_REF][START_REF] Do | On the dynamic deformation of a bar against an obstacle[END_REF][START_REF] Friedman | Variational principles and free-boundary problems[END_REF][START_REF] Gilbarg | Elliptic partial differential equations of second order[END_REF][START_REF] Kim | A one-dimensional dynamic contact problem in linear viscoelasticity[END_REF][START_REF] Lasiecka | Nonhomogeneous boundary value problems for second order hyperbolic operators[END_REF][START_REF] Lebeau | A wave problem in a half-space with a unilateral constraint at the boundary[END_REF][START_REF] Lions | Quelques methodes de resolution des problemes aux limites non lineaires[END_REF][START_REF] Schatzman | A hyperbolic problem of second order with unilateral constraints[END_REF][START_REF] Schatzman | Un probleme hyperbolique du 2eme ordre avec contrainte unilaterale: La corde vibrante avec obstacle ponctuel[END_REF][START_REF] Schatzman | The penalty method for the vibrating string with an obstacle[END_REF] and (3-15), IOt(6)1:::; Mo, for all k.

(3-32) Meanwhile we also have for each k, since Bk is nonnegative and integrable, and 0:::; tP6(x):::; 1, for all x. Since A-BE L""(O,T; L 1 (fl)), it is apparent that