Quentin Rendu

Emmanuel Mignot
email: emmanuel.mignot@insa-lyon.fr

Nicolas Riviere

Barbara Lamberti-Raverot

Sara Puijalon

Florence Piola

Laboratory investigation of Fallopia x bohemica fruits dispersal by watercourses

Keywords: Buoyancy, Settling velocity, Achene, Fruit hydrodynamics, Hydrochory, Helical motion, Shear flow

Introduction

Biological invasions are known to be an important issue in the loss of biodiversity. The japanese knotweed, also referred as genus Fallopia or Reynoutria, is one of the 100 most invasive species [START_REF]GISD: Species profile : Polygonum cuspidatum Sieb. & Zucc. (=Fallopia japonica (Houtt. Dcne[END_REF]. Recent studies showed that the sexual reproduction of the japanese knotweed highly contributes to its dispersion in Belgium [START_REF] Tiébré | Hybridization and morphogenetic variation in the invasive alien Fallopia (polygonaceae) complex in belgium[END_REF], Czech Republic [START_REF] Mandák | Isoenzyme diversity in reynoutria (polygonaceae) taxa: escape from sterility by hybridization[END_REF] and USA [START_REF] Grimsby | Genetic diversity and clonal vs. sexual reproduction in Fallopia spp. (polygonaceae)[END_REF]. The japanese knotweed mainly settles on riverbanks in riparian zones where the hydrochory (i.e. its spreading by watercourses) is the main process of propagation [START_REF] Skoglund | Seed dispersing agents in two regularly flooded river sites[END_REF] even though Varshney et al. [START_REF] Varshney | Unsteady aerodynamic forces and torques on falling parallelograms in coupled tumbing-helical motion[END_REF] state that multiple-winged seeds (such as Fallopia) exhibit complex aerial movements which assist their long distance propagation. Rouifed et al. [START_REF] Rouifed | Achene buoyancy and germinability of the terrestrial invasive Fallopia x Bohemica in aquatic environment: A new vector of dispersion?[END_REF] confirmed that watercourses are an efficient way of dispersion for the japanese knotweed. Note that while rare events such as floods or zoochory are, even though improbable, very efficient for very long distance dispersal [START_REF] Nathan | Long-distance dispersal of plants[END_REF], the present paper focuses on the deterministic long-distance dispersal of japanese knotweed fruits by watercourses.

Hydrochory is mainly governed by the diaspores density, size and shape [START_REF] Chambert | Sorting of seeds by hydrochory[END_REF] and can be split in two main behaviours, depending on the diaspore density: buoyant diaspores that float on the surface and non-buoyant ones that sink. The buoyancy depends on the anatomical structure of the pericarp (the envelope of the fruit) and more particularly the presence of aeriferous tissue (conducting air) in the exocarp, the external layer of the pericarp [START_REF] Hroudová | Effect of anatomical structure on the buoyancy of achenes of two subspecies of Bolboschoenus maritimus[END_REF][START_REF] Praeger | On the buoyancy of the seeds of some britannic plants[END_REF][START_REF] Ridley | The Dispersal of Plants Throughout the World[END_REF]. Moreover, buoyancy is highly correlated with the shape of the diaspores [START_REF] Yoshikawa | Role of seed settleability and settling velocity in water for plant colonization of river gravel bars[END_REF]. Oppositely, even if the size of the achenes highly differs from one individual to another, it did not appear as a main factor in buoyancy studies by Hroudová et al. [START_REF] Hroudová | Effect of anatomical structure on the buoyancy of achenes of two subspecies of Bolboschoenus maritimus[END_REF] and Coops and Van der Velde [START_REF] Coops | Seed dispersal, germination and seedling growth of six helophyte species in relation to water-level zonation[END_REF]. The dispersal of non-buoyant diaspores can be depicted by the well studied mechanisms governing mineral sediment transport in rivercourse [START_REF] Nakayama | Seed deposition of eragrostis curvula, an invasive alien plant on a river floodplain[END_REF]. The fundamental property characterizing the deposit of a non-buoyant seed is the ratio between its settling velocity and the mean stream velocity [START_REF] Chambert | Sorting of seeds by hydrochory[END_REF][START_REF] Yang | Sediment Transport: Theory and Practice[END_REF]. Regarding buoyant seeds, studies were dedicated to the influence of wind, obstacles and vegetation on their dispersal along the free-surface. The results show that diaspores follow the freesurface stream and that differentiation between species is not significant without wind [START_REF] Chambert | Sorting of seeds by hydrochory[END_REF].

Rouifed et al. [START_REF] Rouifed | Achene buoyancy and germinability of the terrestrial invasive Fallopia x Bohemica in aquatic environment: A new vector of dispersion?[END_REF] show that when detaching from the tree and falling in a watercourse, the dry fruits of the japanese knotweed float. While floating they absorb water, become waterlogged and finally sink. They thus exhibit both buoyant and non-buoyant behaviours one after the other. Both floating and non-floating hydrochory mechanisms thus have to be investigated to fully understand the japanese knotweed dispersal. A key parameter for the japanese knotweed fruit dispersal is the floating time after which the fruits are waterlogged enough to sink (i.e. with a density higher than unity). Furthermore, the floating time should be regarded in respect to the watercourse flow rate. Indeed, if the floating time of a diaspore is higher than the ratio of the watercourse length by its bulk velocity, a non buoyant study is irrelevant [START_REF] Andersson | Plant dispersal in boreal rivers and its relation to the diversity of riparian flora[END_REF]. On the other hand, as clearly stated by Varshney et al. [START_REF] Varshney | Unsteady aerodynamic forces and torques on falling parallelograms in coupled tumbing-helical motion[END_REF]:

The falling dynamics of seeds are critical for understanding their dispersal properties. Similarly the settling characteristics of seeds in water may strongly affect the dispersal characteristics and will thus be investigated herein.

The influence of environmental parameters on the vegetation growth has been widely studied [START_REF] Sculthorpe | The biology of aquatic plants[END_REF], however their influence on the floating time and the settling velocity remains poorly documented. Indeed, the floating time is usually measured in controlled laboratory conditions for a unique temperature, under a defined luminosity and for one quality of water while these parameters vary along time and space in a given watercourse.

The aim of the present work is to better understand the influence of environmental parameters on the time and motion characteristics of the fruit as they reach a watercourse, i.e. as it floats and then as it sinks.

In a first section, we aim at evaluating the influence of the environmental parameters on the fruit's floating time: agitation of water, temperature, water quality and enlightenment. Then, the second section is dedicated to the measurement of settling velocity and the three-dimensional motion of the fruits in water at rest. Finally, the falling motion of fruits in a shear flow is investigated. In each of these three sections, the corresponding materials and method are first presented, followed by the results. The last section is dedicated to a general discussion and the conclusions of this work.

Floatability analysis

Achene sampling

The achenes considered herein were collected by the Laboratoire d' Écologie des Hydrosystèmes Naturels et Anthropisés (LEHNA) in december 2012 from Fallopia x bohemica stands growing in Bourgoin-Jallieu in the Département de l'Isère (France). They were randomly selected from several Fallopia shoots of the same stand to reduce variability and then pooled. The achenes were stored in dry conditions at 4°C until they were used for experiments. Preliminary tests showed that these sampling and storage methods do not affect the germination rates [START_REF] Rouifed | Achene buoyancy and germinability of the terrestrial invasive Fallopia x Bohemica in aquatic environment: A new vector of dispersion?[END_REF].

Experimental setup

The experimental setup for the floatability study is made of two 20x30x20cm pools sketched in Fig. 1. These pools are immersed in a large transparent tank of water. Some holes in the bottom of these pools guarantee the same temperature (measure with a thermocouple) and water quality between both pools and the tank. A long and thin pipe (8mm diameter) is fixed all along the bottom's perimeter of the pools and is linked to an air compressor (Fig. 1a). Small holes in this pipe allow Fig. 1: Sketch of the experimental device for free surface studies -a) top view b) front view air bubbles to escape and rise in the pools towards the free surface, thus agitating the volume of water in each pool as shown on Fig. 1b. The pipes are equipped with a valve allowing to adjust the air discharge and then to vary the agitation intensity. With such a device, air bubbles are periodically shed from each hole and generate a vertical motion towards the free surface as they rise. This agitation limits the concentration of the fruits in specific regions of the free surface where they would accumulate and interact with each other.

This whole set-up is placed either in a large hall with a glass-roof, i.e. receiving natural lighting during day time with a non-controlled temperature or in a closed room with artificial light and controlled ambient temperature. The water used in the experiments is either "clean" (degassed tap water, surface tension of 70 mN/m) or "highly polluted" by adding soap to decrease its surface tension (30 mN/m). Sridhar and Reddy [START_REF] Sridhar | Surface tension of polluted waters and treated wastewater[END_REF] collected surface tension data from water from natural, industrial or domestic waters and reveal that in polluted natural watercourses the surface tension decreases from about 73 mN/m (for clear water at 20˚C) to about 50 mN/m.

Both front and rear side of the pools are made of glass in order to ensure an optimal optical access. A 640*480 pixels camera is placed in front of each pool to visualize the bottom of the pools. A mirror located behind the pools and slightly inclined allows to visualize the free surface with the same camera.

For each experiment, one hundred fruits are deposited on the free surface of each pool at t=0 and afterwards a picture is taken every five minutes. A semiautomated image analysis allows to count the number of floating fruits in each pool along the time.

Additionally, another set up consisting of a simple tank allows to study the floating time of immersed fruits (see Table 1). In this case, the fruits are placed under a net placed one centimeter below the water surface, as shown in Fig. 2. In order to study the behaviour of the fruits under different environmental conditions 16 floatability experiments are run with the following varying parameters: agitation of water or immersed fruits, water temperature and surface tension (see columns Fig. 2: Sketch of the specific experimental device used for studying the immersed fruits 2 to 4 in Table 1). Besides, to evaluate the influence of illumination, experiment 8b was run in the dark.

Results

For each configuration of Table 1, the ratio of floating fruits is plotted against time as in Fig. 3 for experiment 1a. Note that apart from experiment 0 for which the fruits do not sink, the figures corresponding to all experiments exhibit the same trend: as they are deposited on the water surface (t=0) the fruits first float, then after a few hours, the first fruits start to sink, the sinking rate increases, reaches a maximum value and decreases. The few fruits remaining on the free-surface at the end of the experiment germinate on the free surface. Note that in Table 1 the subscripts a and b refer to the first and the second pool respectively. t 50% and t 75% are then defined as the times after which 50% and 75% of the fruits have sunk, respectively. Note that for an agitated flow, as the density of the fruit reaches 1, the fruit does not sink towards the bottom but rather remains in the water column and follows the water motion. In such case, a fruit is considered as sunk when it finally settles at the bottom of the pool.

All experimental conditions and corresponding results are summarized in Table 1. Following observations can be made:

-Experiment 0 reveals that in clean water at rest, the fruits do not sink. The sketch in Fig. 5a and the photograph in Fig. 6a reveal that as the fruit is deposited on the water surface, an air bubble is trapped between the water and the fruit which prevents from direct contact and thus from an efficient water absorption. -Experiments 1,2 and 3 reveal that as clean water is agitated, the fruits end up sinking after an average of 40 hours. Moreover these experiments show a small, yet noticeable influence of the strength of agitation on the fruits' floating time: it decreases as the agitation intensity increases (see Fig. 4a). -Experiments 1a, 2a and 3a, driven at the same agitation intensity but at different temperatures reveal that the floating time strongly decreases as the water temperature increases. It ranges from t 50% = 57h for a constant temperature Fig. 3: Time evolution of the ratio of number of floating fruits N over initial number of fruits N 0 . Raw data (squares) is obtained from 1a experiment (see Table 1), linear interpolation (solid line) is used to extract the characteristic floating times t 50% and t 75% at wich 50% and 75% of the initial number of fruits have sunk, respectively of 20°Cto t 50% = 29h for a temperature of 27°C(the same tendency is obtained comparing experiments 2b and 3b, 4 and 5, see Fig. 4b). -Experiments 4 and 5 are driven without agitation in "polluted" water to analyse the effect of surface tension on the fruits buoyancy. First observation is that the fruits actually sink while this was not the case for "clean" water without agitation. Moreover, the average sinking time t 50% appears to be lower for the "polluted" water at rest than for the agitated "clean" water (see Fig. 4c). -Experiments 6 and 7 are driven with fruits maintained below the water surface, performed respectively with "clean" and "polluted" water, showing that the immersion reduces the floating time.

The following conclusions can be made: by comparing experiments 0, 2a, 4, and 6a, it is concluded that the contact between the fruits and the water (agitation or immersion) strongly influences the fruits' floating time (twice shorter for the immersed fruits). The effect of the agitation is thus to limit or cancel the influence of the air trapped below the fruit in water at rest (see Fig. 6a). Nevertheless even with a strongly agitated flow, water contact remains limited to the bottom half of the fruit surface. Also, as the surface tension decreases, the contact between the fruits and the water increases. As a fruit is deposited on the surface of "polluted" water, the center part is situated below the surface, two wings lie on the water whereas the third one is immersed, facing downwards as it is shown in Fig. 5. The pictures in Fig. 6b confirm this behaviour which could explain the increased ability of the fruits to absorb water in a polluted medium. With the immersed fruits, this contact surface between fruit and water is maximized, increasing the water logging and thus reducing the floating time with regards to all other configurations. Furthermore, the immersed fruits sink earlier in polluted water than in clean water, indicating that pollution has an influence on water absorption mechanisms of the fruit. Finally, in all configurations (except in clean water at rest), the floating time decreases as the temperature increases. Temperature is an important factor inducing seed germination: under wet conditions, increased temperature favors and accelerates induction of seed germinations [START_REF] Probert | The role of temperature in germination ecophysiology. Seeds: the ecology of regeneration in plant communities[END_REF][START_REF] Baskin | Seeds: ecology, biogeography, and, evolution of dormancy and germination[END_REF]. The changes occurring during germination (embryo development, use of reserves, alteration of seed coat..) are probably involved in the change of buoyancy observed for the seeds of Fallopia. This may explain why the temperature, which controls these processes, may greatly influence the changes of buoyancy observed.

It should be additionally noted that the influence of the three parameters are cumulative. To conclude this floatability analysis, the floating time reduces as:

the floating fruit/water contact increases the temperature of the water increases the water is polluted

The experimental range of these parameters seems representative of natural watercourse parameters. "Weak agitation" corresponds to limited secondary currents and turbulent energy (just strong enough to evacuate air bubbles) while "strong agitation" corresponds to the presence of, for instance, helical secondary cells in river bends, to boils, emerging boulders, or soaking branches in rivers; surface tension ranges from pure to polluted water in rivers, temperatures from 20°Cto 28°Ccorrespond to seasonal variations of typical river temperatures. These variations of quantitative parameters lead to floating times in the range of 17 to 60 hours. This floating time can be considered as O(10 h) i.e. of the order of some tens of hours.

3 Free fall in water at rest

Experimental setup

The fruits' settling velocity measurements are performed in a 20x18x70cm transparent tank filled with clean water. A 1920*1080 pixels camera and a plane mirror are used to record the three-dimensional fruits' motion. A light source is located behind the tank and the shadow of the fruit is recorded by the camera. Thanks to the plane mirror, the pathline of the fruits in two orthogonal planes is recorded, which permits to reconstruct the 3D trajectory. This technique, known as 3D Particle Tracking Velocimetry (PTV) is presented in Varshney et al. [START_REF] Varshney | Unsteady aerodynamic forces and torques on falling parallelograms in coupled tumbing-helical motion[END_REF]. Because of the small size of the fruits compared to the field view of the camera, the angular position of the fruits during the fall cannot be obtained. To obtain fruits with a density larger than unity, a large amount of dry fruits is first deposited at the surface of the pools (see Sect. 2.2) filled with polluted water (in order to reduce their floating time). When a fruit sinks, its floating time is recorded, it is removed with care from the pool and deposited in a tank filled with clean water, just below the surface to avoid the influence of surface tension. It then sinks within the clean water tank. Through image analysis of the recorded pictures, the three-dimensional pathline of the fruit against the depth is obtained and the mean settling velocity is computed. This allows to investigate the influence of the floating time of a fruit on its mean settling velocity.

Results

Settling velocity

The mean fruits settling velocity is plotted against the floating time of the fruit in Fig. 7 for 42 fruits. For a given floating time, these results show a great variability of velocity from a fruit to another, mostly due to specific geometry parameters. Averaged velocities (blacks squares) are thus also included in the figure. No clear influence of floating time on the settling velocity can be observed. Hence the influence of the floating time will be discarded in further analyses of sinking fruits. Averaging all data from Fig. 7, the ensemble averaged vertical settling velocity v y finally equals 0.41 cm/s with a standard deviation of 0.25 cm/s where • stands for the ensemble averaging of the mean velocity values. By comparison with literature, non buoyant fruits' settling velocity ranges from 0.9 to 2.7 cm/s [START_REF] Chambert | Sorting of seeds by hydrochory[END_REF]. The reasons for which the present fruits have a lower settling velocity is likely to be related to their lower relative density which may be closer to 1.

3D motion Varshney et al. [START_REF] Varshney | Unsteady aerodynamic forces and torques on falling parallelograms in coupled tumbing-helical motion[END_REF] review the four well-known type of bodies falling motion including "autogyrating, "rolling gyration", "floating" and "undulation" motion. In the present case, when sinking towards the bottom in water at rest, all fruits exhibit a coupled autogyration/rolling gyration. Such a coupled motion has already been observed on tulip seeds falling in air by McCutchen [START_REF] Mccutchen | The spinning rotation of ash and tulip tree samaras[END_REF] and is similar Fig. 4: Influence of (a) agitation, (b) temperature and (c) surface tension on fruits floating time -t 50% (respectively t 75% ) denotes the time after which 50% (respectively 75%) of the initially deposited fruits have sunk to the Fig. 2 presented in [START_REF] Varshney | Unsteady aerodynamic forces and torques on falling parallelograms in coupled tumbing-helical motion[END_REF]. From now on this motion will be referred as "helical" for simplicity. In our experiments, a great variability is observed on the helix radius ([1-15 cm], see Fig. 8) and on the vertical mean velocity v y ([±0.5cm/s]) without any correlation between them. A helical pathline with very large radius is presented in Fig. 8a whereas a helical pathline with small radius is presented in Fig. 8b. The smallest radius recorded is 5 mm and the largest one is 70 mm.

Oppositely, limited variability is observed on the ratio of mean tangential velocity over vertical velocity v ⊥ /v y . The ensemble average value of this ratio is v ⊥ /v y = 0.25 and the corresponding standard deviation is about 0.1.

The following conclusions can be raised for this free-fall section: (i) the settling velocity of the fruits does not depend on its floating time ; (ii) the mean settling velocity is low compared to the literature of non buoyant fruits ; (iii) the fruits exhibit a helical motion with a high variability of vertical velocity and radius ; (iv) the mean tangential velocity of the helical motion is about 25% of the mean vertical velocity. The present section aims at investigating the characteristics of the trajectory and velocity of the fruit as it settles within a steady shear flow, representative of a river flow.

Experimental setup

The experiments are performed in a b = 25cm wide and 8 m long horizontal flume with a steady subcritical flow and a streamwise bulk velocity U = 3cm/s, a water depth h = 22cm, a Froude number F r = 0.006 and Reynolds number Re = 4U bh (b+2h)ν = 9600. The geometrical aspect ratio of the flow section b/h is about equal to 1, which tends to generate secondary currents with a downward velocity in the center section of typical magnitude equal to a few percent of U , according to Nezu and Nakagawa [START_REF] Nezu | Turbulence in open channels[END_REF]. x refers to the streamwise axis and y to the vertical axis. u refers to the velocity of the flow and v to the fruit velocity.

About one hour after it sinks within the polluted water pool used in Sect. 2, each tested fruit is deposited slightly below the surface at the upstream section of the test region (x = 0). The one hour delay ensures a density slightly larger than 1 and thus that the fruit actually sinks in the sheared flow. The two-dimensional pathline of the fruit is recorded by a 1900*1280 pixels camera (50Hz acquisition frequency). Using image analysis, the fruit trajectory is recovered and the corresponding velocity components computed. A total of 24 fruits are successfully deposited and the following data refer to 24 trajectories as the one shown in Fig. 9, similar to Fig. 5 from Zhu et al. [START_REF] Zhu | Settling distance and incipient motion of sandbags in open channel flows[END_REF] with sandbags settling in a shear flow.

The vertical profile of streamwise mean velocity along the centerline of the channel in the middle of the test region of about 1m long is measured with a micro-ADV (Acoustic Doppler Velocimeter, from Nortek) and is plotted in Fig. 10. The vertical profiles of streamwise and vertical turbulent intensity (not shown here) appear to fairly match the trend previously presented by Nezu and Nakagawa [14, Fig. 4.9].

Results

Fig. 10 shows the vertical profile of mean streamwise flow velocity u x and of the ensemble averaged fruits vertical velocity v y and streamwise velocity v x . It was assessed that these ensemble average velocities converged for sample sizes higher than 20. It appears that:

-The mean streamwise velocity of the flow is about uniform in the center part of the water column (equal to U 0 = 3.5cm/s) and decreases when approaching the bed due to the boundary layer and when approaching the free surface due to the velocity-dip phenomenon, i.e. the maximum velocity appears not at the free-surface, but rather just below it in open-channel flows [START_REF] Nezu | Turbulence in open channels[END_REF]. -Along the fruits motion (i.e. from the free surface to the bed) the ensemble averaged fruit vertical velocity v y first increases (for y/h > 0.75) and then becomes about constant over the depth, as in the free fall experiments and similarly to measurements from Elkholy and Chaudhry [START_REF] Elkholy | Drag and added-mass coefficients of large sandbags[END_REF]. These observation indicate a balance between weight and vertical hydrodynamic forces in the center part of the channel (0.15 < y/h < 0.75). The magnitude of the vertical velocity v y is about 0.7cm/s in the center part of the channel, which is higher than the ensemble average mean settling velocity in water at rest (see Sect. 3.2). This difference of velocity is attributed to the variations of fruit density between the two experiments and the downward secondary current in the center section of the flume. Finally, when approaching the bed, the fruit vertical velocity decreases. This deceleration is attributed to lift forces appearing in the strongly sheared flow close to the bed. -The profile of ensemble averaged streamwise velocity of the fruits v x is qualitatively similar to that of the streamwise flow velocity u x . However, discrepancies are observed: (i) the fruit velocity gradient is constant in a large portion of the water column (0.15 < y/h < 0.85) ; (ii) v x exceeds u x in the upper part of the channel (y/h > 0.6) and is smaller in the lower part (y/h < 0.5).

In order to investigate the interaction between the flow and the fruits, and thus explain the ensemble average fruit trajectory, it would be useful to compute the lift and drag forces acting on the fruit, based on the measured pathline. This methodology is proposed for instance by Elkholy and Chaudhry [START_REF] Elkholy | Investigations on the trajectory of large sandbangs in open channel flow[END_REF] on large sandbags. However, estimating exact volume and projected surfaces of the fruits, as well as their density is particularly challenging. The approach selected herein reduces to the analysis of the velocity of several individual fruits represented by gray lines in Fig. 10.

It appears that:

-The vertical fruit velocity shows a small yet noticeable variability from one fruit to another. The fluctuations of vertical velocity during the fall for a given fruit can be either small or moderate depending on the individual. -The streamwise fruit velocity shows a large variability from one fall to another and also moderate to large fluctuations during the fall for each fruit. These oscillations are attributed to helical motions as previously observed in the tank with water at rest: the vertical settling velocity produces a horizontal lift force which results in a horizontal circular motion during the settlement. This circular motion alternatively adds or substracts a streamwise relative velocity to the mean fruit velocity. Hence, the fruit streamwise velocity can alternatively exceed or remain slower than the local flow velocity. It can bee seen in Fig. 10 that the fruit streamwise velocity is oscillating around the ensemble averaged value with a difference of velocity of about 0.15U 0 . The tangential velocity of the helical motion was found in Sect. 3.2 to be about 0.25v y which would lead to at most ±0.25v y fluctuations on the streamwise velocity, represented by a double arrow in Fig. 10. Consequently, the streamwise velocity fluctuations for a given fruit are two or three times higher than those measured during free fall. This indicates a complex three-dimensional asymmetric motion of the falling fruits due to the interaction of the helical motion with the one-dimensional turbulent shear flow.

Finally, the main results regarding the streamwise motion obtained in this section are: (i) the ensemble averaged streamwise velocity of the fruits differs from the streamwise velocity of the flow ; (ii) the three-dimensional motion differs from the helical motion in free fall described in Sect. 3.2.

Discussion and Conclusion

The typical cycle of a Fallopia x bohemica fruit falling on a river surface after it detaches from the tree is to 1) fly in the air, 2) deposit and float on the free surface and 3) sink and 4) eventually deposit. The present work allowed to analyse the behaviours of these achenes as they float and settle in controlled flow conditions. It appeared that the floating time strongly depends on the characteristics of the flow: temperature, agitation and surface tension of the water. The floating time decreases as the temperature increases and as the contact with water increases (through agitation, low surface tension or immersion).

The distances covered by the fruit before it stops and germinates (namely, the mean dispersion length) can be roughly estimated as follows. In the present experiment, the mean floating duration ranges between 14 and 66 hours (see Table 1). In a river with a bulk velocity of 50cm/s, where the surface velocity equals about 62.5cm/s, assuming an average floating time of 40 hours (corresponding to test 1a with strong turbulence and a warm temperature varying between 23 and 28°C), the distance covered by the floating fruit reaches about 90km. However many processes may limit the fruit transport such as a storage in dead zones, in a wake downstream obstacles, in recirculation zones downstream a sudden channel expansion or upstream a sudden width decrease, which would result in shorter travel distances.

Then, assuming that the river is 2m deep, and that the mean vertical settling velocity equals 0.7cm/s (Fig. 10), the settling duration will be of about 5 minutes (neglecting secondary currents and large scale structures) and the fruit will thus cover a streamwise distance of about 150m. This rough estimate reveals that the distance covered within the water column is negligible compared to the floating distance. As a consequence, in order to fairly estimate the dispersal length of a fruit specie, a knowledge of the floating duration as performed in section 2 is essential. Indeed, as depending on the temperature, agitation, and surface tension, the floating time can vary of about ±100%, the mean distance covered by the fruit species varies of the same amount.

Predictive models for the sinking phase would require an estimation of, for instance, drag and lift coefficients of the fruits. Such coefficients are however very dependent on the individual and specie of knotweed, so as on the river flow characteristics. It would be also difficult to account for additional phenomena such as change in river sections, soaking branches, boils, etc. A sinking predictive model was thus not undertaken in the present work, especially as present work clearly shows that the settling time is negligible compared to the floating time. In other words, sinking phase has a negligible influence compared to the floating phase of several tens of hours.

Fig. 5 :

 5 Fig. 5: Schemes of fruit position when falling or deposited on water surface with a high tension surface (clean water, a) and low tension surface (polluted water, b)

Fig. 6 :

 6 Fig. 6: Pictures of buoyant fruits deposited on clean water (a) and on polluted water (b)

Fig. 7 :

 7 Fig. 7: Mean settling velocity (v y : •) obtained for 42 free fall experiments as a function of the floating time. Average ensemble value ( v y ) and standard deviation σ are computed for each hour interval ( and bars) and for the whole set of fruit (dash line)

Fig. 8 :Fig. 9 :

 89 Fig.8: Top and front view of the 3D helical motion of two sinking fruits. For case a, only a short portion of the settling motion is recorded due to a helical radius much larger than the window recorded by the camera.

Fig. 10 :

 10 Fig. 10: Dimensionless streamwise flow velocity (u x : •) plotted along the dimensionless height y/h as well as ensemble average streamwise ( v x : ♦) and absolute vertical ( v y : ) fruit velocities, grey solid lines represent the streamwise and vertical velocity of a selection of ten individuals ; h = 22.0 cm, U 0 = 3.5 cm.s -1

Table 1 :

 1

	Name	Agitation	Twater(˚C) Surface tension	t 50%	t 75%
				(mN/m)
	0	None	23-28*	70	-	-
	1a	Strong	23-28*	70	41h	57h
	1b	Weak	23-28*	70	49h	62h
	2a	Strong	20	70	57h	80h
	2b	Mean	20	70	66h	80h
	3a	Strong	27	70	29h	36h
	3b	Mean	27	70	31h	40h
	4	None	20	30	42h	59h
	5	None	27	30	21h	24h
	6a	Immersed 20	70	28h	38h
	6b	Immersed 27	70	25h	31h
	7a	Immersed 20	30	22h	29h
	7b	Immersed 27	30	15h	17h
	8a	Strong	22-26*	70	51h	65h
	8b (in the dark) Strong	22-26*	70	51h	65h

Conditions and results of the 17 buoyancy experiments. a and b refer to the first and the second pool respectively. * refers to the non controlled temperature experiments for which only min-max values of the water temperature are given.


Héctor Maravillo
email: hmaravillo@hotmail.com

José-Fernando Camacho-Vallejo

Justo Puerto
email: puerto@us.es

Martine Labbé
email: mlabbe@ulb.ac.be

A market regulation bilevel problem: a case study of the Mexican petrochemical industry

Keywords: Market regulation, Bilevel programming, Petrochemical industry

In this paper, a bilevel programming model is proposed to study a problem of market regulation through government intervention. One of the main characteristics of the problem herein analyzed is that the government monopolizes the raw material in one industry, and competes in another industry with private firms for the production of commodities. Under this scheme, the government controls a state-owned firm to balance the market; that is, to minimize the difference between the produced and demanded commodities. On the other hand, a regulatory organism that coordinates private firms aims to maximize the total profit by deciding the amount of raw material bought from the state-owned firm. Two equivalent single-level reformulations are proposed to solve the problem. The first reformulation is based on the strong duality condition of the lower level and results in a continuous nonlinear model. The second reformulation resorts to the complementarity slackness optimality constraints yielding a mixed-integer linear model. Additionally, three heuristic algorithms are designed to obtain good-quality solutions with low computational effort. In this problem, the feasible region of the dual problem associated to the follower is independent from the leader's decision. Therefore, the proposed heuristics exploit this particular characteristic of the bilevel model. Moreover, the third heuristic hybridizes the other two algorithms to enhance its performance. Extensive computational experimentation is carried out to measure the efficiency of the proposed solution methodologies. A case study based on the Mexican petrochemical industry is presented. Additional instances generated from the case study are considered to validate the robustness of the proposed heuristic algorithms. Numerical results indicate that the hybrid algorithm outperforms the other two heuristics. However, all of them demonstrate to be good alternatives for solving the problem. Additionally, optimal solutions of all the instances are obtained by using good quality solutions (given by the hybrid algorithm) as initial solutions when solving the second reformulation via a general purpose solver.

Introduction

Market regulation through government intervention has appeared in different ways, such as, fiscal policies (taxes and subsidies), anti-monopoly legislation, price control, quantity control, and nationalization of firms [START_REF] Wiens | Government firm regulation of a vertically integrated industry[END_REF]. Two of the main motivations for regulating the economy is to seek social equity and to handle market failures [START_REF] Stiglitz | Regulation and failure[END_REF]. A specific manner of government intervention is through the participation of state-owned firms, which are under the government control and compete with the private firms. The goal to regulate the market by participation has been very noticeable in developing countries, mainly in raw materials industries [START_REF] Vernon | Uncertainty in the resource industries: the special role of state-owned enterprises[END_REF]. For instance, a list of countries in which market regulation by participation exists in oil industry is presented in [START_REF] Kobrin | The nationalisation of oil production[END_REF], [START_REF] Minor | The demise of expropriation as an instrument of ldc policy 1980-1992[END_REF], [START_REF] Guriev | Determinants of nationalization in the oil sector: A theory and evidence from panel data[END_REF]. At this point we would like to emphasize that this situation occurred in Mexico in the petrochemical industry from 1958 until 2014.

The first paper that considers the idea of regulate the market by including a state-owned firm is [START_REF] Merrill | Government firms in oligopoly industries: A short-run analysis[END_REF]. In that paper, a short-term analysis that concerns with the entrance of a stateowned firm into a three-firms oligopoly market is done. They conclude that the existence of a state-owned firm may improve the performance of the market, which it is also shown in [START_REF] Blackwell | The state and the macro-industrial economy: Lessons from the aerospace sector[END_REF] in the aerospace market. Following up with the existence of state-owned firms, [START_REF] Harris | Government enterprise: An instrument for the internal regulation of industry[END_REF] considers a situation in which the difference between the production and the production level is made up by the government. The novelty in that paper is to realize that state-owned firm acts as the dominant decision-maker, that is, it has complete information of the market and announces its decision. Hence, each private firm reacts to this decision by establishing its output level so that its marginal cost equals the price. Furthermore, [START_REF] Sertel | Regulation by participation[END_REF] explicitly considers the state-owned firm as a leader in a Stackelberg game. Later, in [START_REF] Cremer | The public firm as an instrument for regulating an oligopolistic market[END_REF], a mixed oligopoly model that helps to compare the differences among the nationalization of a state-owned firm and the entrance of a new state-owned firm is studied. This last paper also analyzes the cost effectiveness of the state-owned firms in the long-term and moreover, the authors consider for the first time a budgetary constraint associated with the state-owned firm to guarantee a minimum profit.

Literature reviews regarding mixed oligopolies can be found in [START_REF] De Fraja | Game theoretic models of mixed oligopoly[END_REF] and [START_REF] De Fraja | Mixed oligopoly: An overview[END_REF]. In those reviews the common characteristics of these models are identified. Particularly, the context of the problem, the cases when the government provides complete information, the goal of the state-owned firm, the technology assumptions, the costs structure, and the hierarchy among the firms (in case of Stackelberg games). Another literature review that deals with some foundations for a theory of mixed oligopoly markets is presented in [START_REF] Nett | Mixed oligopoly with homogeneous goods[END_REF]. It is important to remark that all the above mentioned papers coincide in that there is a lack of an unified and accepted general mixed oligopolies modeling framework because each model pursues its own characteristics and set different basic hypothesis. At this point, it is also important to emphasize that, to the best of our knowledge and after an intensive literature review, the intervention of the government via state-owned firms to regulate an industry conformed by two interrelated markets has not been studied before. However, these characteristics appear in real situations, as it will be described in the case study herein considered.

The problem studied in this paper is inspired by a situation where a state-owned firm controls the monopoly of a market but this is interrelated with a second market in which state-owned and private firms compete. The problem considers two interdependent industries, in which the first one produces the necessary supplies for the second one. The government has the monopoly of the first industry, while in the second industry a state-owned firm competes with remaining private firms in the production of commodities. Private firms aim to maximize their profit. Government aims as a general social welfare to balance the market. The latter is achieved by minimizing the difference between the supply and demand of the final commodities.

Given the privileged position that the government occupies in this market as the leader decision-maker, the government decides the amount of commodities that will be produced at the state-owned firms and the amount of supplies that will be offered to the private firms. By doing this, the government will regulate indirectly the production. One can easily appreciate a hierarchical relationship between both economic agents, in which the government makes a decision and the private firms reacts to it affecting the balance in the market. Furthermore, technical and technological issues in the raw material market limit the production, for example, the scarcity of natural resources and the existing production capacity of the firms.

The hierarchy and inter-relationship among the actors of this situation allow to formulate a bilevel programming model, in which the government will act as the leader and the private firms as the follower. Other bilevel models in which the government intervenes to regulate social aspects are [START_REF] Demirci | Designing an intervention strategy for public-interest goods: The california electric vehicle market case[END_REF] and [START_REF] Whittaker | Spatial targeting of agri-environmental policy using bilevel evolutionary optimization[END_REF]. In the former, the government employs an intervention policy based on subsidies in the automotive market, while in the latter, the government apply taxes via an agro-environmental policy imposed to the agriculturists. Additionally, the general model proposed by us perfectly fits to the petrochemical industry in several developing countries in the last decades and it was the case of Mexico for almost 56 years, where it was apparently run without a theoretically recognizable model.

The main contributions of this paper can be listed as follows: a novel mathematical bilevel programming model to study a mixed-oligopoly market, three heuristic algorithms based on an iterative exploration of vertices in the inducible region of the bilevel problem, an extensive computational study analyzing the performance of our exact reformulations and the proposed heuristic algorithms, and a case-study based on the Mexican petrochemical industry.

The remainder of the paper is organized as follows. Section 2 defines the mathematical programming bilevel model and sets the notation. Two reformulations of the bilevel model that reduce it to a single-level one are presented in section 3. The first one is a continuous 3 non-linear problem, and the second one reduces to a mixed integer linear program. In both cases, the resulting programs are hard to solve for medium to large sizes. Then, section 4 describes three proposed heuristic algorithms to provide good quality solutions for the bilevel model. Section 5 reports the numerical results according to the case study under consideration, and compares the results obtained through the reformulations and the algorithms. Also, the results obtained from additional experimentation with random instances are summarized. Conclusions and recommendations for future research are given in section 6.

A market regulation bilevel problem

Problem's description

The problem herein studied considers an industry conformed by two economic markets, one of them associated to raw material and the other one of the final commodities. The supplies (raw material) of the second economic market are produced in the first one. In this industry, there is a state-owned firm vertically integrated, that is, that produces in both economic markets. The state-owned firm monopolizes the production of the supplies in the first economic market, but in the other one, this firm competes against the private firms. All the firms have a maximum production capacity and the state-owned firm requires to obtain a minimum profit (there is a lower bound over its net profit). The objective of this latter firm is to balance the market by its intervention. To achieve this goal, it minimizes the lack and surplus of the offered commodities with respect to the demand. Hence, this firm determines its production of commodities and the amount of supplies to be offered to the private firms. On the other hand, private firms' goal will be to maximize their benefit.

The decisions taken by the state-owned firm limit the admissible production by the private firms, and the decisions of the private firms affect the achievement of the government.

As it is mentioned before, the government has the monopoly of the supplies production, and it is assumed that there exists a centralized organism that regulates the supplies demand. This organism determines the amount of supplies to each private firm, which is common in economic markets with government intervention. A typical example occurs when the government creates organisms to regulate the competition among firms and establishes particular contracts with each private firm fixing the supplies to be sold to each one. Another case is when the private firms get together and create a centralized mechanism to which the responsibility of distributing supplies among them is delegated. In the latter case, that mechanism seeks for a global welfare of the market; this is the case of common lands cooperative organization in agriculture. Under this scheme, the government will be the leader and the centralized mechanism the follower.

Mathematical formulation

In this section, the mathematical formulation of the problem described in the previous section is formally introduced. Let I be the set of commodities and let J be the set of private firms. Each commodity i ∈ I is sold at a price p i , and it has a demand d i in the market. To produce a commodity i ∈ I, it costs c G i to the state-owned firm and c E ij to a private firm j ∈ J. The government fixes a minimum profit t to be obtained by the stateowned firm. The amount of primary resources that a private firm j ∈ J needs to produce a unit of commodity i ∈ I is denoted by b ij and the amount of supply required to produce a unit of commodity i ∈ I is denoted by a ij . Both types of firms have a limited production capacity. The state-owned firm has a maximum production capacity q B i to obtain supplies of the commodity i ∈ I and a maximum amount q A i to produce for each commodity i ∈ I.

Also, each private firm j ∈ J has a maximum production capacity m j .

In order to present our mathematical programming formulation, we will use the following decision variables. The leader's decision variables are:

x i , which denotes the production of the state-owned firm for each commodity i ∈ I.

z i , which denotes the supply offered by the government to the private firms to produce a commodity in ∈ I.

The follower's decision variables are:

y ij , which denotes the amount produced by private firm j ∈ J of a commodity i ∈ I.

In our model, non-negative auxiliary variables are introduced to express the leader's objective function. Let r i be the shortage of commodity i ∈ I and let s i be the corresponding surplus.

With the above elements, the proposed bilevel programming formulation to model the regulation of the economic market described above results as follows: min

x,z,r,s,y i∈I

(r i + s i ) (1)
s.t. j∈J y ij + x i d i + r i -s i = 1 ∀i ∈ I (2) i∈I (p i -c G i )x i ≥ t (3) 0 ≤ x i ≤ q A i ∀i ∈ I (4) 0 ≤ z i ≤ q B i ∀i ∈ I (5) r i ≥ 0 ∀ ∈ Ii (6) s i ≥ 0 ∀i ∈ I (7) y ∈ argmax j∈J i∈I (p i -c E ij )ȳ ij (8) s.t. j∈J a ij ȳij ≤ z i ∀i ∈ I (9) i∈I b ij ȳij ≤ m j ∀j ∈ J (10) ȳij ≥ 0 ∀i ∈ I, ∀j ∈ J (11)
The leader's problem is defined by ( 1)-( 8), in which (1) represents the leader's objective function that aims to minimize the inefficiency of the market, namely the overall shortage and surplus production with respect to the total demand. Constraint (2) seeks to balance the demand of each commodity with the corresponding production of the state-owned and private firms. In (3), it is ensured that the state-owned firm obtains a minimum required profit from its commodities production. Constraints (4) and ( 5) establish the production capacity associated with the commodities and offered supplies, respectively. The non-negativity of the auxiliary variables is expressed in ( 6) and [START_REF] Merrill | Government firms in oligopoly industries: A short-run analysis[END_REF]. Constraint (8) plays a key role in this model, due to the fact that it requires that the value of the follower's variables is given by the optimal solution of another mathematical programming problem. The follower's problem is defined by equations ( 8)- [START_REF] Cremer | The public firm as an instrument for regulating an oligopolistic market[END_REF], which aims to maximize the total profit of all the private firms. The production of the private firms is limited by the amount of supplies provided by the government, this is enforced by [START_REF] Harris | Government enterprise: An instrument for the internal regulation of industry[END_REF]. Constraint [START_REF] Sertel | Regulation by participation[END_REF] imposes that each private firm's production cannot exceed its maximum capacity. Finally, [START_REF] Cremer | The public firm as an instrument for regulating an oligopolistic market[END_REF] expresses the non-negativity of the follower's variables.

In order to have a well-defined formulation for the proposed bilevel model, it is necessary to make an assumption regarding multiple optimal solutions that may appear in the follower's problem. We assume the optimistic approach defined in [START_REF] Dempe | Bilevel programming: A survey[END_REF]. In other words, among all the optimal follower's solutions, the one that minimizes the leader's objective function is selected. The optimistic approach is suitable for the situation under study since it may express the existence of a certain cooperation degree among the government and the mechanism (committee) that represents the private firms.

Exact solution methodologies

A common methodology that applies, at times, to solve a bilevel programming model is to transform it to a single-level reformulation. In order to achieve this goal, the characteristics of the follower's problem have to be exploited. Once a leader's solution is established, the leader's variables are fixed in the follower's problem. Particularly, z is going to be considered as a parameter in the problem defined by ( 8)- [START_REF] Cremer | The public firm as an instrument for regulating an oligopolistic market[END_REF]. Therefore, the follower's problem results in a linear programming problem that has a corresponding dual, in which α and β correspond to the associated dual variables. The dual problem associated with the follower's problem is as follows:

min α,β i∈I α i z i + j∈J β j m j (12) s.t. a ij α i + b ij β j ≥ (p i -c E ij ) ∀i, ∀j (13)
α i ≥ 0 ∀i (14) β j ≥ 0 ∀j (15)
For both, the primal and dual of the follower's problem, the space of feasible solutions is bounded. Hence, the fundamental duality theorem states that both problems have optimal solutions and their objective function values are equal [START_REF] Bard | Practical Bilevel Optimization, Nonconvex Optimization and Its Applications[END_REF]. Hence, two single-level reformulations that are equal to the bilevel model are presented in this section. The first reformulation ensures that the optimal solution of the reformulated model are in the inducible region of the bilevel problem by using a corollary of the strong duality theorem.

The second reformulation substitutes a non-linear constraint of the first reformulation via the complementarity slackness constraint. The first reformulation results in a continuous non-linear model and the second one is a mixed-integer linear model.

Non-linear reformulation based on the strong duality condition

The first reformulation consists in adding the constraints associated with the dual of the follower's problem, that is, [START_REF] De Fraja | Game theoretic models of mixed oligopoly[END_REF][START_REF] De Fraja | Mixed oligopoly: An overview[END_REF][START_REF] Nett | Mixed oligopoly with homogeneous goods[END_REF][START_REF] Demirci | Designing an intervention strategy for public-interest goods: The california electric vehicle market case[END_REF], and to use the necessary and sufficient optimality conditions to ensure that the follower's decision belongs to the set of rational reactions. Therefore, a constraint that equals the follower's primal and dual objective functions is added. For a predefined leader's fixed value of z the resulting constraint is as follows:

j∈J i∈I (p i -c E ij )y ij = i∈I α i z i + j∈J β j m j (16)
The first proposed reformulation consists in the following non-linear single-level problem: min

x,z,r,s,y,α,β i∈I

(r i + s i ) (Ref.1) s.t. j∈J y ij + x i d i + r i -s i = 1 ∀i ∈ I (17) i∈I (p i -c G i )x i ≥ t (18) 0 ≤ x i ≤ q A i ∀i ∈ I (19) 0 ≤ z i ≤ q B i ∀i ∈ I (20) r i ≥ 0 ∀i ∈ I (21) s i ≥ 0 ∀i ∈ I (22) j∈J a ij y ij ≤ z i ∀i ∈ I ( 23
) i∈I b ij y ij ≤ m j ∀j ∈ J (24)
y ij ≥ 0 ∀i ∈ I, j ∈ J ( 25
)
a ij α i + b ij β j ≥ (p i -c E ij ) ∀i ∈ I, j ∈ J ( 26
)
α i ≥ 0 ∀i ∈ I (27) β j ≥ 0 ∀j ∈ J (28) j∈J i∈I (p i -c E ij )y ij = i∈I α i z i + j∈J β j m j (29)
Note that the linearity of the model is lost in [START_REF] Johnson | La transferencia del manejo de la irrigación en México: una estrategia para lograr la sostenibilidad de los distritos de riego[END_REF]. Further, it can be represented in the following manner since weak duality implies that the other inequality is always satisfied.

- j∈J i∈I (p i -c E ij )y ij + i∈I α i z i + j∈J β j m j ≤ 0 (30)

Mixed-integer linear reformulation based on complementarity slackness

The second reformulation is based on the complementarity slackness constraints that ensures the optimality of the follower's problem. A single-level reformulation is recovered by introducing the following constraints:

α i z i - j∈J a ij y ij = 0 ∀i ∈ I (31) β j m j - i∈I b ij y ij = 0 ∀j ∈ J (32)
y ij a ij α i + b ij β j -(p i -c E ij ) = 0 ∀i ∈ I, j ∈ J (33)
Constraints ( 31), (32), and (33) are non-linear, nevertheless they can be linearized in a straightforward manner by introducing positive big-M constants M 1 i , M 2 i , and M 3 i for all i ∈ I. Also, binary variables γ i , δ j , and ij for all i ∈ I, j ∈ J are included. Constraints (31)-(33) are replaced by the following constraints:

α i ≤ M 1 i γ i ∀i ∈ I (34) z i - j∈J a ij y ij ≤ M 1 i (1 -γ i ) ∀i ∈ I (35)
β j ≤ M 2 i δ j ∀i ∈ I, j ∈ J (36) m j - i∈I b ij y ij ≤ M 2 i (1 -δ j ) ∀i ∈ I, j ∈ J ( 37
)
y ij ≤ M 3 i ij ∀i ∈ I, ∀j ∈ J (38) a ij α i + b ij β j -(p i -c E ij ) ≤ M 3 i (1 -ij ) ∀i ∈ I, j ∈ J (39)
Therefore, this reformulation results in the following single-level linear mixed-integer programming model. min x,z,r,s,y,α,β,γ,δ, i∈I

(r i + s i ) (Ref.2) s.t. j∈J y ij + x i d i + r i -s i = 1 ∀i ∈ I (40) i∈I (p i -c G i )x i ≥ t (41) 0 ≤ x i ≤ q A i ∀i ∈ I (42) 0 ≤ z i ≤ q B i ∀i ∈ I (43) r i ≥ 0 ∀i ∈ I (44) s i ≥ 0 ∀i ∈ I (45) j∈J a ij y ij ≤ z i ∀i ∈ I (46) i∈I b ij y ij ≤ m j ∀j ∈ J ( 47
)
y ij ≥ 0 ∀i, j ∈ J ( 48
)
a ij α i + b ij β j ≥ (p i -c E ij ) ∀i ∈ I, j ∈ J ( 49
)
α i ≥ 0 ∀i ∈ I (50) β j ≥ 0 ∀j ∈ J (51) α i ≤ M 1 i γ i ∀i ∈ I (52) z i - j∈J a ij y ij ≤ M 2 i (1 -γ i ) ∀i ∈ I (53)
β j ≤ M 1 j δ j ∀j ∈ J (54) m j - i∈I b ij y ij ≤ M 2 j (1 -δ j ) ∀j ∈ J ( 55
)
y ij ≤ M 1 ij ij ∀i ∈ I, j ∈ J ( 56
)
a ij α i + b ij β j -(p i -c E ij ) ≤ M 2 ij (1 -ij ) ∀i ∈ I, j ∈ J (57) γ i , δ j , ij ∈ {0, 1} ∀i ∈ I, j ∈ J (58)

Adjusting the value of the big M 's

In the problem under study, parameters a and b are always positives, which means that if more commodity production exists, more supplies are required and more production capacity is used in the private firms. Moreover, if term p i -c E ij is negative, private firms are having losses, which would lead to their exit from the economic market.

Definition. Define the upper bound U B(•) ∈ R of a vector as a real number that is is greater than or equal to all of the components of that vector.

To adjust the value of M 1 i for all i ∈ I, first we seek for an upper bound for α i , such that α belongs to an optimal solution of the follower's dual problem. The worst scenario occurs when

β j = 0 ∀j ∈ J, then α i ≥ p i -c E ij a ij
for all i ∈ I, ∀j ∈ J must hold due to (49). Since [START_REF] De Fraja | Game theoretic models of mixed oligopoly[END_REF] aims to minimize, the optimal solution occurs at equality. To obtain the upper bound for M 2 i for all i ∈ I, constraint (53) is considered. Hence, it is possible to state the following proposition.

Proposition 1.

1. The upper bound for any optimal value of α i for all i ∈ I is as follows:

U B(α i ) = max j∈J { p i -c E ij a ij } 2.
The upper bound for constraint (53) is as follows:

U B z i -j∈J a ij y ij = q B i .
Analogously that for M k i , the upper bounds for M k j i = 1, 2 and k = 1, 2 are computed.

Proposition 2.

1. The upper bound for any optimal value of β j for all j ∈ J is:

U B(β j ) = max i∈I { p i -c E ij b ij }.
2. The upper bound for constraint (55) is:

U B m j -i∈I b ij y ij = m j .
To bound the values of M k ij for all i ∈ I, j ∈ J and k = 1, 2, the following proposition is stated.

Proposition 3.

1.

The upper bounds for y ij for all i ∈ I, j ∈ J are:

U B(y ij ) = m j b ij , i.e. y ij ≤ m j
b ij for all i ∈ I, j ∈ J.

2.

The upper bound of constraint (57) is:

U B(a ij α i + b ij β j -(p i -c E ij )) = a ij U B(α i ) + b ij U B(β j ).
Therefore, the adjusted values for M k i , M k j , and M k ij for all i ∈ I, j ∈ J, k = 1, 2 based on the computed upper bounds are as follows:

M 1 i = max j∈J { p i -c E ij a ij }, ∀i ∈ I (59) M 2 i = q B i , ∀i ∈ I (60) M 1 j = max i∈I { p i -c E ij b ij }, ∀j ∈ J (61) M 2 j = m j , ∀j ∈ J (62) M 1 ij = m j b ij , ∀i ∈ I, j ∈ J (63) M 2 ij = a ij U B(α i ) + b ij U B(β j ), ∀i ∈ I, j ∈ J. (64)

Heuristic solution methodologies

Usually, the reformulations introduced above present computational limitations for largesize instances. Therefore, alternative approaches are required to solve the bilevel problem.

One common approach is to design a heuristic algorithm to obtain good quality feasible solutions with lower computational burden. In this section, three tailor-made heuristic algorithms are proposed to solve the problem under study. The first two heuristics exploit the particular structure of the bilevel problem, while the third one hybridizes the other two. The hybrid heuristic contains key aspects from the previously proposed heuristics to enhance its performance.

Note that it has been already emphasized that the polyhedron that defines the feasible region of the dual problem associated with the follower is independent from the leader's variables. Therefore, if all the vertices of this polyhedron were known, equality ( 16) and the dual constraints of the first reformulation described in section 3.1 could be substituted by a set of constraints guaranteeing that the equality ( 16) would be achieved in one vertex.

Following the latter idea, the dual variables are replaced by parameters that represents those known vertices. This is the approach exploited in the three proposed heuristic algorithms.

Extreme points iterated algorithm (EPIA)

In this algorithm, vertices of the dual polyhedron are iteratively generated, and a mixedinteger programming problem named the Master Problem (MP) is solved for each new vertex. The optimal solution of the MP ensures a feasible solution of the bilevel problem.

The algorithm stops when no improvement in the leader's objective function is obtained.

Let P be the polyhedron associated with the dual constraints of problem ( 13)- [START_REF] Demirci | Designing an intervention strategy for public-interest goods: The california electric vehicle market case[END_REF]. Con- ). Let F be the leader's objective function [START_REF] Wiens | Government firm regulation of a vertically integrated industry[END_REF]. Also, let f and g be the follower primal [START_REF] Blackwell | The state and the macro-industrial economy: Lessons from the aerospace sector[END_REF] and dual [START_REF] De Fraja | Game theoretic models of mixed oligopoly[END_REF] objective functions, respectively. Consider χ = (x, z, r, s) as a vector that groups the leader's variables and Υ = (y 1 , ..., y |J| ) as a vector that groups the follower's variables. Define ϕ as the set of vertices in P that has been already explored. At each iteration K = |ϕ| is updated using k ∈ K as the index of the explored vertices. U B(z i ) represents an upper bound of variable z i , which can be naturally fixed as q B i . Hence, the MP is defined as follows:

sider v k = (α k , β k ) ∈ P, k = 1, ...,
min x,z,r,s,y,λ i∈I (r i + s i ) (65) s.t. j∈J y ij + x i d i + r i -s i = 1 ∀i ∈ I (66) i∈I (p i -c G i )x i ≥ t (67) 0 ≤ x i ≤ q A i ∀i ∈ I (68) 0 ≤ z i ≤ q B i ∀i ∈ I (69) r i ≥ 0 ∀i ∈ I (70) s i ≥ 0 ∀i ∈ I (71) j∈J a ij y ij ≤ z i ∀i ∈ I (72) i∈I b ij y ij ≤ m j ∀j ∈ J (73)
y ij ≥ 0 ∀i ∈ I, j ∈ J (74) i∈I α k i z i + j∈J β k j m j -(1 -λ k )( i∈I α k i U B(z i ) + j∈J β k j m j ) ≤ j∈J i∈I (p i -c E ij )y ij ∀k ∈ K ( 75
) k∈K λ k = 1 (76) λ k ∈ {0, 1} ∀k ∈ K (77)
Constraints ( 75)-(77) guarantee that f (Υ) = g(χ, v k ) holds only for a single vertex, in which λ is an auxiliary binary variable. The obtained solution Υ * belongs to the rational reaction set. Hence, the optimal solution of MP is a feasible solution of the bilevel problem.

The process in which the dual problem ( 12)-( 15) is solved for a fixed vector z obtaining an optimal solution v = (α, β) is represented as D(z) → v. Analogously, solve MP for a set of vertices ϕ is denoted as MP(ϕ) → (χ, Υ). The result is a pair of solutions (χ, Υ). The pseudocode of the EPIA is presented in algorithm 4.1.

Algorithm 4.1 Extreme points iterated algorithm

Step 1. Initialization: ϕ = ∅, ρ = ∞;

Step 2. Find the initial vertices

D(0) → v 1 . D(q B i ) → v 2 . ϕ = ϕ ∪ {v 1 , v 2 }; Step 3. Solve the master problem: MP(ϕ) → (χ k , Υ k ). π = F (χ k , Υ k ); Step 4. Find a vertex k: D(z k ) → v k ;
Step 5. Evaluate leader's objective function.

• If ρ > π then ϕ = ϕ ∪ {v k }, ρ = π. Return to Step 3; • If ρ ≤ π then stops. Output: π, (χ k , Υ k )
Note that optimality of the bilevel problem is not guaranteed by this algorithm. The main reason is that during the process, the algorithm may not improve the leader's objective function. This may occurs when the value z k obtained at iteration k produces a vertex v k that is already in set ϕ, or simply when

F (χ k , Υ k ) ≥ F (χ k-1 , Υ k-1 ).

Follower's gap penalization algorithm (FGPA)

The general idea of FGPA is to find a feasible solution by using vertices of the follower's dual problem, but the dual admits infeasible solutions. The algorithm consists on itera-tively solving a mixed-integer programming problem, named MMP, that is a modification of the MP described in previous section. The MMP permits the existence of a gap between follower's primal objective function and the dual objective function value obtained by evaluating the considered vertices. However, this gap is penalized in the leader's objective function aiming to obtain a good feasible solution.

Since the primal and dual objective values of the follower can be different, the algorithm may explore other dual vertices using solutions that are in the constraint region of the bilevel problem, but not necessarily in the inducible region. Hence, to ensure feasibility, at the end of the iterations another problem named Resulting problem (RP) must be solved. RP is a linear single-level problem defined by constraints ( 1)-( 7) of the leader, constraints ( 9)-( 11) of the follower, and a constraint that equals f (Υ) with the value of the dual objective function g evaluated with the last obtained value of z (ẑ = z k ). Therefore, ψ = g(z k ) is a parameter used in:

j∈J i∈I (p i -c E ij )y ij = ψ (78)
In MMP, the gap is normalized by a constant M and multiplied by a coefficient µ to regulate the impact in the leader's objective function, that is, to balance the supply and demand. The value of M can be bounded by the maximum among all the optimal values of the dual problem obtained at each iteration. The modified master problem (MMP) is follows: min

x,z,r,s,y,λ,ε i∈I

(r i + s i ) + ε M µ (79) s.t. j∈J y ij + x i d i + r i -s i = 1 ∀i ∈ I (80) i∈I (p i -c G i )x i ≥ t (81) 0 ≤ x i ≤ q A i ∀i ∈ I (82) 0 ≤ z i ≤ q B i ∀i ∈ I (83) r i ≥ 0 ∀i ∈ I (84) s i ≥ 0 ∀i ∈ I (85) j∈J a ij y ij ≤ z i ∀i ∈ I (86) i∈I b ij y ij ≤ m j ∀j ∈ J (87)
y ij ≥ 0 ∀i ∈ I, j ∈ J (88) k∈K λ k = 1 (89) λ k ∈ {0, 1} ∀k ∈ K (90) ε k + j∈J i∈I (p i -c E ij )y ij = i∈I α k i z i + j∈J β k j m j ∀k ∈ K (91) ε k ≥ 0 ∀k ∈ K (92) ε ≥ 0 (93) ε ≥ ε k -(1 -λ k )M ∀k ∈ K (94)
Therefore, the RP that is used to obtain a solution in the inducible region is: min

x,r,s,y i∈I

(r i + s i ) (95) s.t. j∈J y ij + x i d i + r i -s i = 1 ∀i ∈ I (96) i∈I (p i -c G i )x i ≥ t (97) 0 ≤ x i ≤ q A i ∀i ∈ I (98) r i ≥ 0 ∀i ∈ I (99) s i ≥ 0 ∀i ∈ I (100) j∈J a ij y ij ≤ ẑi ∀i ∈ I (101) i∈I b ij y ij ≤ m j ∀j ∈ J (102)
y ij ≥ 0 ∀i, j ∈ J (103) j∈J i∈I (p i -c E ij )y ij = ψ (104)
To refer to the MMP, define MMP(ϕ, M ) → (χ, Υ) as in the previous algorithm. The leader's objective function (79) is represented by F * . Also, RP(ϕ, ψ, ẑ) → (χ, Υ) denotes the process of solving RP using the set of vectors ϕ and the parameters ψ and ẑ. The pseudocode of FGPA is shonw in algorithm 4.2.

Hybrid algorithm (HYBA)

As mentioned before, the EPIA cannot ensure to obtain a new vertex in each iteration nor a vertex that improves the leader's objective function. An intuitive idea is to use the FGPA to avoid this issue. The latter algorithm has de advantage of exploring solutions that does not belong to the inducible region, by doing this, unexplored vertices are obtained.

The main idea of the HYBA is to perform the steps of EPIA until stops. When this occurs, a subroutine that solves the MMP is performed to find a new vertex aiming to reach a different feasible solution. If the new exploration improves the leader's objective function,

Algorithm 4.2 Follower's gap penalization algorithm

Step 1.

Initialization: ϕ = ∅, ρ = ∞, M = M ax{1, D(q B i )}; Step 2. Solve the modified master problem: PMM(ϕ, M ) → (χ k , Υ k ). π = F * (χ k , Υ k ); Step 3. Find a vertex k: D(z k ) → v k . M = M ax{M, D(z k )};
Step 4. Evaluate leader's objective funcion.

• If ρ > π then ϕ = ϕ ∪ {v k }, ρ = π. Return to Step 2;

• If ρ ≤ π then go to Step 5.

Step 5. Obtain a feasible solution:

ψ = D(z k ). RP(ϕ, ψ, ẑ) → (χ * , Υ * ). π = F (χ * , Υ * ); Output: π, (χ * , Υ * )
then the FGPA continues until the stopping criterion is reached. The pseudocode of the proposed hybrid algorithm is shown in algorithm 4.3.

Algorithm 4.3 Hybrid algorithm

Step 1. Initialization: ϕ = ∅, ρ = ∞;

Step 2. Find the initial vertices:

D(0) → v 1 . D(q B i ) → v 2 . ϕ = ϕ ∪ {v 1 , v 2 }, M = M ax{D(0), D(q B i )}; Step 3. Solve the master problem: MP(ϕ) → (χ k , Υ k ). π = F (χ k , Υ k ); Step 4. Find a vertex k: D(z k ) → v k . ϕ = ϕ ∪ {v k }, M = M ax{M, D(z k )};
Step 5. Evaluate leader's objective function.

• If ρ > π then Go to Step 6;

• If ρ ≤ π then MMP(ϕ, M ) → (χ k * , Υ k * ). D(z k * ) → v k * . ϕ = ϕ ∪ {v k * }; MP(ϕ) → (χ k , Υ k ). π = F (χ k , Υ k ). D(z k ) → v k , ϕ = ϕ ∪ {v k }, M = M ax{M, D(z k )}.
Step 6. Re-evaluating leader's objective function:

• If ρ > π then ρ = π. Return to Step 3; • If ρ ≤ π then Stops. Output: π, (χ k , Υ k )

Computational experimentation

In this section, the results obtained from an extensive computational experimentation of our solution methodologies are reported. To evaluate both single-level reformulations and

Case study

In Mexico, the government intervention in the petrochemical industry existed by many decades. The main characteristic was that the government monopolized the extraction of the main raw material for this industry, namely petroleum and some other derivatives called basic petrochemicals. Another characteristic of this situation was that state-owned firms competed against private firms in the market of final commodities (secondary petrochemicals).

The institutional framework defined by law the petrochemical industry as the one that performs chemical or physical processes to product compounds from petroleum natural hydrocarbons or from the products derived from refinement operations. Some of these products may serve as raw materials to the industry, and they were classified as part of the basic petrochemicals. The remaining products were included into the secondary petrochemical category [START_REF]Reglamento de la ley reglamentaria del artículo 27 constitucional en el ramo del petróleo, en materia de petroquímica[END_REF].

Specifically, in a law from 1958 [START_REF]Ley reglamentaria del artículo 27 constitucional en el ramo del petróleo[END_REF], it was established that only the government was allowed to exploit the hydrocarbons related with the oil industry, which was concerned with the production, warehousing, distribution, and sales of the petroleum derivatives that can be considered as raw materials (basic petrochemicals) for the industry. However, for the production of secondary petrochemicals products, state-owned and private firms were allowed to be involved in. Therefore, Petróleos Mexicanos (PEMEX) and its subsidiaries was in charge of this industry [START_REF]Acuerdo que dispone que corresponde a petróleos mexicanos u organizaciones o empresas subsidiarias o asociadas a la misma institución, la elaboración de diversos productos de petroquímica básica[END_REF], and from 1992, also their descentralized departments [START_REF]Ley orgánica de petróleos mexicanos y organismos subsidiarios[END_REF].

The classification of basic and secondary petrochemical products varies by law from year to year [START_REF]Reglamento de la ley reglamentaria del artículo 27 constitucional en el ramo del petróleo[END_REF], [START_REF]Acuerdo que dispone que corresponde a petróleos mexicanos u organizaciones o empresas subsidiarias o asociadas a la misma institución, la elaboración de diversos productos de petroquímica básica[END_REF], [START_REF]Reglamento de la ley reglamentaria del artículo 27 constitucional en el ramo del petróleo, en materia de petroquímica[END_REF], [START_REF]Resolución que clasifica los productos petroquímicos que se indican, dentro de la petroquímica básica o secundaria[END_REF], [START_REF]Resolución que clasifica los productos petroquímicos que se indican, dentro de la petroquímica básica o secundaria[END_REF], [START_REF]Resolución que clasifica los productos petroquímicos que se indican, dentro de la petroquímica básica o secundaria[END_REF], [START_REF]Decreto por el que se reforma la ley reglamentaria del artículo 27 constitucional en el ramo del petróleo[END_REF].

To delimit the scope and size of our case study instances based on the above situation, we consider the number of economic units involved in the manufacturing of organic chemical basic products registered by the National Institute of Statistics and Geography in Mexico (INEGI by its acronym in Spanish) during the economical census conducted in 1999, 2004, 2009, and 2014. The biggest number of these units was registered in 2004, which was 159. This value is used as an upper bound of the private firms dedicated to secondary petrochemical industry. For delimiting the amount of commodities considered, we consider the number of commodities classified as basic and secondary petrochemical by the Mexican legislation in different years. This information is summarized in table 1. Obtained by using data from [START_REF]Acuerdo que dispone que corresponde a petróleos mexicanos u organizaciones o empresas subsidiarias o asociadas a la misma institución, la elaboración de diversos productos de petroquímica básica[END_REF], [START_REF]Resolución que clasifica los productos petroquímicos que se indican, dentro de la petroquímica básica o secundaria[END_REF], [START_REF]Resolución que clasifica los productos petroquímicos que se indican, dentro de la petroquímica básica o secundaria[END_REF], [START_REF]Resolución que clasifica los productos petroquímicos que se indican, dentro de la petroquímica básica o secundaria[END_REF] To complete the instances, we extracted data from the statistical report of the energy sector used by the Mexican Secretary of Energy. The data of 47 petrochemical commodities from 1980 to 2002 were obtained. The information consists of the demand, prices, and productive capacity. For each commodity i, the demand d i was generated from a uniform distribution with the following parameters: the average of the demands in that period of time minus the standard deviation, and the average of the demands in that period of time plus the standard deviation. The price p i and the government production capacity q A i were randomly generated between 1 and the maximum price or the maximum capacity for each commodity i, respectively.

The government costs c G i were obtained as the product of the price of each commodity i times a random number between 0.22 and 0.60. The latter range was defined based on the ratio between the production costs and the PEMEX (the Mexican state-owned firm) total income registered from the years 1988-2000, 2011-2013, and 2015-2016. The minimum profit that the state-owned firm must achieved was computed as the 30% of the total benefit that the government may obtain if all the petrochemicals market demand was satisfied by PEMEX.

The state-owned firm has a production capacity q B i regarding raw materials i, which was computed as the product of the maximum production coefficient of the private firms and the maximum capacity of the state-owned firm for producing commodities. Under this assumption, the private firm with best technology could match the production capacity of the state-owned firm for an specific commodity.

Production costs c E ij of private firm j are defined as the product of random coefficients in [0.784, 0.884] multiplied times the price of each commodity i. That range corresponds to the minimum and maximum average values of the ration between the production costs and the income per year of five real private firms. These private firms operate in the Mexican petrochemical sector: Alpek, Kuo (UEN synthetic rubber), Kuo (UEN polystyrene), Mexichem, and Pochteca group.

The production technical coefficients a ij are randomly generated between [0.085, 2.111],

which correspond to the range determined by the average plus/minus one standard deviation of all the petrochemical substances considered in the report presented by CEPAL [? ]. The production levels b ij are generated from the interval [START_REF] Wiens | Government firm regulation of a vertically integrated industry[END_REF]95]. The lower bound is natural since these coefficients must be strictly positives, while the upper bound is equal to the estimated average of the production levels for each substance considered in the petrochemical facilities of PEMEX.

Finally, the production capacity m j of the private firms were randomly generated between 4, 665 and 20, 825, which came from the average of the production capacity of five petrochemical complex of PEMEX.

Based on the above factors and ranges, we generated 30 instances for each one of the following sizes: As it is mentioned above, the instances generated to analyze the case study were taken as a basis for constructing another set of synthetic instances to test the performance of our algorithms on a different data set. This new set of instances has the same sizes but were randomly generated using arbitrary ranges for each parameter. Furthermore, the generation scheme guarantees feasibility of the problem, that is, ensure that the state-owned firm has capacity to produce all the demand.

• |I| =

Numerical results

Our computational experiment consists on solving both subsets of instances described above, that is, the case-study and the synthetic instances. All of them were solved by the two reformulations presented in section 3 and by the three heuristic algorithms proposed in section 4. The leader's objective function value and the required time were registered.

Among the two exact methods, only Ref.2 was able to optimally solve all the tested instances.

Hence, these values were used to compute the optimality gap obtained by the heuristic algorithms. A maximum CPU time of 1000 seconds was set to the solver for solving Ref.1, while Ref.2 and the heuristic algorithms we did not fix a time limit since in all cases the required time was rather small.

We observe from our results, that Ref.1 was not able to solve all the instances within the time limit. We also observe that the case-study instances are harder to be solved than the synthetic ones. For example, for the larger instance sizes (150×200), Ref.1 did not even find feasible solutions for the problem. This may be due to the fact that synthetic instances were well-balanced so that, in most cases, shortages and surplus are both zero. Summarizing, the number of instances solved to optimality by each reformulation is shown in table 2. The following tables 3 and 4 report all the results of our experiment. The results are organized in five column blocks corresponding to the five solution methods that are compared.

In each block, we include two columns with the objective function values (F) and the CPU time (t). In these tables, averages of the registered values for each size of the case-study and synthetic instances are shown. It is worth mentioning that when Ref.1 was not able to solve the instance, that instance was omitted for the computation of the corresponding average value. averages as the size of the instance increases. Note that for most of the synthetic instances, the optimal value is zero, i.e., the supply is perfectly balanced with the demand. Hence, neither shortages nor surplus exist due to the generation mechanism of these instances.

Regarding the heuristic algorithms, it can been observed from table 4 that EPIA and HYBA report values, on average very close, to the optimal ones for all the sizes of the synthetic instances. However, for the case-study instances, HYBA reports the best results among all the three heuristic algorithms (see table 3). Moreover, the performance of EPIA and FGPA significantly depends on the size of the instance. In terms of the average required time, FGPA requires less time for the case-study instances. This finding may be derived from the fact that the algorithm explores points not necessarily in the inducible region, which may help to approach the optimal solution in a faster manner. On the contrary, HYBA is the algorithm that requires the more time, but this may be obvious since it performs all the steps of EPIA and at each iteration it may solve upon four extra mathematical programs. Nevertheless, the registered times are acceptable for a problem of this nature.

To compare the efficiency and the quality of the solutions obtained by the heuristic algorithms, the optimality gap is computed using the values obtained by Ref.2. Given that there are many optimal values equal to zero, the optimality gap (GAP) measures the relative deviation from the optimal value, and it is computed by the following formula:

GAP = (Optimal value) -(Obtained value) (Obtained value) × 100% (105)
Also, an analogous formula is used to compute the relative savings in time (%t). That is, the reduction in computational time when solving the problem using the heuristics rather than Ref. The GAP and %t average values obtained for both type of instances are shown in tables 5 and 6, respectively. HYBA is the algorithm that shows the best quality in the feasible solutions obtained.

For the case study instances, the average optimality gap was lower than 8% for all the sizes (see table 5); while for the synthetic instances, the average optimality gap was zero (see table 6). FGPA presents a very good average gap in comparison with EPIA in the case-study instances, but this behavior is opposite in the synthetic instances. Concerning the savings in the required time, the three heuristic algorithms showed significant savings, which is improved as the size of the instance increases. FGPA is the one that evidence more savings for the case-study instances. It is convenient to mention that for the synthetic instances of size To support the latter findings, the number of optimal solutions obtained by each heuristic algorithm are displayed in table 7. It can be seen from that table that HYBA is the algorithm that was able to obtain more optimal solutions. Also, it an be observed that FGPA obtained more optimal solutions for the case-study instances than EPIA, but for the synthetic instances the behavior was the opposite. These results confirm the suitability for combining both heuristic algorithms to create HYBA.

Conclusions and research directions

This paper studies a market regulation situation, in which the government controls the distribution of raw material and competes against private firms in the production of final commodities. In this market, the government has a privileged position since it determines its production and the amount of raw material offered to private firms and its goal is to balance the market. Although this type of market situation is common in some economies, to the best of our knowledge and after an intensive literature revision, we did not find a similar problem in the literature. The situation that motivated our paper comes from the actual situation in Mexican petrochemical industry with PEMEX from 1958 till 2014. However, this case is not exclusive from the petrochemical industry. Also in the agricultural industry, the Mexican government monopolized the production, imports, and distribution of fertilizers from 1970 to 1986 through the state-owned company named FERTIMEX [START_REF] Janvry | Cooperation on Agriculture, Reformas del sector agrícola y el campesinado en México, Estrategias para mitigar la pobreza rural en América Latina y El Caribe[END_REF]. The government also regulated the water supply for agricultural consumption by fixing prices, defining production goals, and restricting the land transactions [START_REF] Johnson | La transferencia del manejo de la irrigación en México: una estrategia para lograr la sostenibilidad de los distritos de riego[END_REF]. In addition, this type of economic regulation has been applied in other countries. For example, in [START_REF] Vernon | Uncertainty in the resource industries: the special role of state-owned enterprises[END_REF], it is stated that there is a tendency to use government instruments in developing countries, mainly in raw material industries. An indicator of the latter can be seen in the large amount of oil expropriations in the 70's (see [START_REF] Kobrin | The nationalisation of oil production[END_REF], [START_REF] Minor | The demise of expropriation as an instrument of ldc policy 1980-1992[END_REF], [START_REF] Guriev | Determinants of nationalization in the oil sector: A theory and evidence from panel data[END_REF]). Also, nationalizations occurred in Bolivia, Ecuador, Venezuela, and Russia in 2006 (see [START_REF] Guriev | Determinants of nationalization in the oil sector: A theory and evidence from panel data[END_REF]).

Motivated by the lack of models to analyze these type of markets, this paper develops a bilevel programming approach that fits naturally to the problem. To solve the bilevel problem herein proposed, two single-level reformulations are presented. The first reformulation uses the lower level necessary and sufficient optimality conditions to ensure the global optimum of the bilevel problem. As a result of this reformulation, a continuous non-linear problem is obtained. The second reformulation substitutes the above-mentioned conditions by the complementarity slackness conditions generating a MILP. These reformulations are able to solve optimally limited size instances. Therefore, three heuristic algorithms are proposed to find good quality feasible solutions in a reasonable computational time. The first two algorithms (EPIA and FGPA) are tailored for this problem, and exploits the fact that the polyhedron associated to the dual problem of the lower level remains the same for any leader's decision. The third algorithm (HYBA) is a combination of the previous ones.

Extensive computational experimentation is performed to test the exact and heuristic approaches proposed in this paper. Two sets of instances are used: case-study instances, that arise from the Mexican petrochemical industry, and synthetic instances adapted from the previous ones. As a result from the computational experimentation, it can be stated that the reformulation that applies the complementarity slackness constraints is more efficient for the tested instances. On the other hand, the performance of EPIA and FGPA strongly depends on the characteristics of the instance. In addition, the hybridization of both heuristic algorithms (HYBA) yields to the best results in terms of effectiveness and solutions quality.

Moreover, HYBA found the optimal solution for all the synthetic instances, and obtained optimality gaps less than 8% in the case-study ones. It is worth to mention that HYBA found the optimal value in at least 50% for the first four types of instance sizes.

Finally, it can be concluded from our computational experiment, that the most efficient manner to find optimal solutions for the bilevel problem is by obtaining a good quality feasible solution through HYBA and then using it as an initial solution to solve Ref.2. By doing this, a significant computational time reduction is achieved.

Note that this bilevel problem may have multiple optimal solutions, each of them offers different alternatives to the upper level decision-maker (government) for producing and offering raw material to the lower level decision-makers (private firms). Therefore, this characteristic is economically worthy to be explored. Another possible research direction is to neglect, in the model, the assumption there exists an organism that manages the cooperation among private firms. Hence, a new element will emerge in the model since competition among private firms will also exist. As a result of this, a bilevel problem with one leader and multiple followers will appear. Then, in a natural way, a Stackelberg game in which the follower's problem consists of a generalized Nash game fits to analyze this situation.

  |P | as the vertices of P , in which |P | define the number of vertices in that polyhedron with α k = (α k 1 , ..., α k |I| ) and β k = (β k 1 , ..., α k |J|

  the performance of the proposed heuristic algorithms, a set of 360 instances was used. The first subset of 180 instances corresponds to a case based on a real situation occurred in the Mexican petrochemical industry between 1958-2014. The second subset of 180 instances was randomly generated to test the efficiency of our algorithms in more general data sets. All the computational experiments were conducted in a personal computer Dell Inspiron 5558, with the following characteristics: a processor Intel(R) Core i3 with 2.10 GHz, 6.00 GB of RAM under Windows 10 operative system. The mixed-integer linear reformulation based on the complementarity slackness condition (hereafter Ref.2), and the three proposed heuristic algorithms were implemented in C++ using Microsoft Visual Studio 2017. The optimizer used was CPLEX 12.7.1. On the other hand, the non-linear reformulation based on the strong duality condition (hereafter Ref.1) was implemented in AMPL using Baron 18.5.8 as optimizer.

10 , |J| = 10 • |I| = 25 , |J| = 25 • |I| = 25 ,

 1010252525 |J| = 75 • |I| = 50, |J| = 100 • |I| = 75, |J| = 125 • |I| = 150, |J| = 200

Table 1 :

 1 Classification of petrochemical commodities

	Year of the	Basic	Secondary
	classification	petrochemicals	petrochemicals
	1960	16	-
	1986	34	26
	1989	8	13
	1992	20	35

Table 2 :

 2 Number of instances solved to optimality

	Instances	Size	Ref.1 Ref.2
		10 × 10	15	30
		25 × 25	0	30
	Case-study	25 × 75 50 × 100	0 0	30 30
		75 × 125	0	30
		150 × 200	0	30
		10 × 10	30	30
		25 × 25	29	30
	Synthetic	25 × 50 50 × 100	30 29	30 30
		75 × 125	28	30
		150 × 200	13	30

Table 3 :

 3 Summarized results obtained from the case-study instances

	Size	Ref.1		Ref.2	EPIA	FGPA	HYBA
		F	t (s)	F	t (s)	F t (s)	F t (s)	F	t (s)
	10 × 10	0.561 255.56 0.506	0.18 0.774 0.23 0.825 0.16 0.524	0.41
	25 × 25	14.225 500.26 1.000	1.86 1.570 0.30 1.171 0.35 1.020	0.71
	25 × 75	12.435 500.46 0.527	4.36 1.025 0.48 0.565 0.39 0.552	1.02
	50 × 100 211.696 558.80 1.051	53.15 2.041 2.39 1.132 1.14 1.091	5.18
	75 × 125 431.138 501.24 1.546 268.89 2.339 8.89 1.811 3.68 1.621 20.69
	150 × 200	-	-3.279 2754.68 4.129 61.82 5.037 12.06 3.443 112.25

Table 4 :

 4 Summarized results obtained from the synthetic instances

	Size	Ref.1	Ref.2	EPIA	FGPA	HYBA
		F	t (s) F	t (s)	F t (s)	F t (s) F t (s)
	10 × 10	0	0.26 0	0.14 0.007 0.08 0.006 0.10 0 0.09
	25 × 25	0.012 29.01 0	1.43	0 0.11 0.021 0.14 0 0.11
	25 × 75	0 10.94 0	3.65	0 0.17 0.019 0.20 0 0.17
	50 × 100 0.057 73.60 0 35.52	0 0.29 0.170 0.31 0 0.28
	75 × 125	0 118.76 0 90.48	0 0.60 0.086 0.76 0 0.57
	150 × 200 1.657 414.55 0 193.48	0 1.78 0.443 1.74 0 1.57
	We observe from tables 3 and 4 that Ref.1 (recall that Ref.1 is a continuous non-linear
	global optimization problem, while Ref.2 is a MILP) gets values, on average, closed to those
	obtained by Ref.2 for the synthetic instances. However, for the case-study instances the
	situation is different and the average solution values and times move away from the optimal

  2. A negative value indicates that the heuristic consumes more time than the exact

	reformulation. The formula used is presented next:
	%t =	(Required time by Ref.2) -(Required time by a heuristic) (Required time by a heuristic)	× 100%	(106)

Table 5 :

 5 Evaluating the effectiveness of the heuristic algorithms in the case-study instances

	Size	EPIA		FGPA	HYBA
		GAP	%t	GAP	%t GAP	%t
	10 × 10	42.27%	12% 29.60%	15% 7.68%	-48%
	25 × 25	40.15% 746% 11.70%	534% 1.67% 190%
	25 × 75	53.78% 1156% 19.67% 1101% 7.04% 445%
	50 × 100 54.29% 3418% 10.72% 5418% 6.02% 1176%
	75 × 125 33.20% 5269% 12.26% 10696% 5.64% 1499%
	150 × 200 19.29% 9210% 17.20% 30206% 4.92% 3324%

Table 6 :

 6 Evaluating the effectiveness of the heuristic algorithms in the synthetic instances

	Size	EPIA		FGPA		HYBA
		GAP	%t	GAP	%t GAP	%t
	10 × 10	6.67%	79% 10.00%	42%	0%	69%
	25 × 25	0% 1332% 46.67% 1038%	0% 1283%
	25 × 75	0% 2189% 50.00% 1895%	0% 2192%
	50 × 100	0% 13714% 90.00% 12176%	0% 14466%
	75 × 125	0% 16773% 90.00% 14272%	0% 18485%
	150 × 200	0% 11561% 100.00% 11377%	0% 12781%

  10 × 10, HYBA showed a negative saving, this implies that Ref.2 was faster to solve these instances. But, this is expected since HYBA solves at least one linear model and one MILP model, and Ref.2 only solves one MILP model. In spite of that, the advantage of HYBA over the reformulation is evident for medium and large-size instances.

Table 7 :

 7 Number of optimal solutions obtained for each heuristic algorithm for each solution scheme. The final column "Total" sums the time required for the hybrid algorithm and Ref.2 using an initial good solution. As it was expected, the advantages of using this scheme are more notorious for medium and large-size instances. The best example of this finding is given by the instances of size 150 × 200, in which the time

	Instance	Size	EPIA FGPA HYBA
		10 × 10	7	11	20
		25 × 25	0	12	19
	Case-study	25 × 75 50 × 100	5 1	13 13	19 15
		75 × 125	0	4	9
		150 × 200	0	0	1
		10 × 10	28	27	30
		25 × 25	30	16	30
	Synthetic	25 × 75 50 × 100	30 30	15 3	30 30
		75 × 125	30	3	30
		150 × 200	30	0	30
	5.3. Using heuristic solutions to obtain the optimal
	The good performance of the heuristic algorithms lead us to analyze the idea of use the
	near-optimal obtained solutions as an input to Ref.2 seeking to enhance their process. By
	doing this, we expect to significantly reduce the computational time required to optimally

solve an instance. The results obtained from this experimentation are shown in table 8.

The columns with label "Ref.2", "HYBA", and "Ref.2 w/initial sol." display the average required time

Table 8 :

 8 Average required time when Ref.2 uses an initial heuristic solution

	Instance	Size	Ref.2 HYBA Ref.2 w/initial sol. Total
		10 × 10	0.18	0.33	0.07	0.40
		25 × 25	1.86	0.82	0.39	1.21
	Case study	25 × 75 50 × 100	4.36 53.15	1.54 6.07	0.80 4.24 10.31 2.34
		75 × 125	268.89	20.24	19.73 39.97
		150 × 200 2754.68 196.38	390.15 586.53
		10 × 10	0.14	0.09	0.04	0.13
		25 × 25	1.43	0.12	0.10	0.22
	Random	25 × 75 50 × 100	3.65 35.52	0.19 0.33	0.49 1.15	0.68 1.48
		75 × 125	90.48	0.46	2.22	2.68
		150 × 200 193.48	1.95	7.28	9.23

Acknowledgements

José-Fernando Camacho-Vallejo has been partially supported by the program of professional development of professors with the grant PRODEP/511-6/17/7425 for research stays during his sabbatical year. Justo Puerto has been partially supported by Spanish/FEDER grants number MTM2016-74983-C02-01 (MINECO). Martine Labbé has been partially supported by the Fonds de la Recherche Scientique -FNRS under Grant(s) no PDR T0098.18.