B Recur

H Balacey

J Bou Sleiman

J B Perraud

J.-P Guillet

A Kingston

P Mounaix

Ordered subsets convex algorithm for 3D terahertz transmission tomography

Keywords:

We investigate in this paper a new reconstruction method in order to perform 3D Terahertz (THz) tomography using a continuous wave acquisition setup in transmission mode. This method is based on the Maximum Likelihood for TRansmission tomography (ML-TR) first developed for X-ray imaging. We optimize the Ordered Subsets Convex (OSC) implementation of the ML-TR by including the Gaussian propagation model of THz waves and take into account the intensity distributions of both blank calibration scan and dark-field measured on THz detectors. THz ML-TR reconstruction quality and accuracy are discussed and compared to other tomographic reconstructions.

Introduction

Terahertz technology (between 0.1 and 4 THz) is now a well-established tool to achieve contactfree and non-destructive testing (NDT) [START_REF] Kawase | Non-destructive terahertz imaging of illicit drugs using spectral fingerprints[END_REF][START_REF] Shen | Detection and identification of explosives using terahertz pulsed spectroscopic imaging[END_REF][START_REF] Chan | Imaging with terahertz radiation[END_REF][START_REF] Huang | Tissue characterization using terahertz pulsed imaging in reflection geometry[END_REF][START_REF] Fukunaga | Terahertz spectroscopy applied to the analysis of artists materials[END_REF][START_REF] Caumes | Terahertz tomographic imaging of xviiith dynasty egyptian sealed pottery[END_REF][START_REF] Ospald | Aeronautics composite material inspection with a terahertz time-domain spectroscopy system[END_REF]. In the field of 3D imaging, X-Ray computed tomography (CT) is an omnipresent technique which provides 3D visualization of dense materials such as human bodies, biological tissues or industrial samples. Nevertheless, this powerful technique cannot be easily applied to soft materials owing to the low absorption of the X-Ray radiation. Alternatively, THz radiation offers a good penetration depth through various soft materials such as plastics, papers or paintings, allowing direct applications in non-destructive inspection. Most demonstrations to date have been performed in 2D, however 3D THz CT is an emerging technique which has been investigated during the last decade [START_REF] Ferguson | T-ray computed tomography[END_REF][START_REF] Wang | T-ray imaging and tomography[END_REF][START_REF] Wang | Pulsed terahertz tomography[END_REF][START_REF] Awad | Transmission terahertz waveguide-based imaging below the diffraction limit[END_REF][START_REF] Yin | Wavelet based local tomographic image using terahertz techniques[END_REF][START_REF] Brahm | Volumetric spectral analysis of materials using terahertz-tomography techniques[END_REF][START_REF] Recur | Investigation on reconstruction methods applied to 3D terahertz computed tomography[END_REF][START_REF] Recur | Propagation beam consideration for 3d thz computed tomography[END_REF][START_REF] Recur | THz radiations for tomographic inspection[END_REF][START_REF] Jp | Review of Terahertz tomography techniques[END_REF].

Low NEP (Noise Equivalent Power) and fast detectors have been developed and demonstrated in laboratory but transfer to industry is quite limited. Thus detectors are usually monoto few-pixels large, resulting in a long raster scanning process to acquire 2D images. Additionally, in order to perform 3D THz CT, several radiographs have to be acquired around the sample at different viewing angles.

Tomographic algorithms allowing accurate reconstruction from a limited number of projections are essential in order to overcome these current technical limitations. An important consideration in 3D THz CT concerns the choice of the reconstruction method to be able to visualize the different cross-sectional images and the final 3D volume of the sample. Iterative reconstructions developed for X-ray CT, such as the Simultaneous Algebraic Reconstruction Technique (SART) [START_REF] Andersen | Simultaneous algebraic reconstruction technique (SART) : A superior implementation of the ART algorithm[END_REF] and the Ordered Subsets Expectation Maximization (OSEM) [START_REF] Shepp | Maximum likelihood reconstruction for emission tomography[END_REF][START_REF] Hudson | Accelerated image reconstruction using ordered subsets of projection data[END_REF], have been investigated for 3D THz CT in [START_REF] Recur | Investigation on reconstruction methods applied to 3D terahertz computed tomography[END_REF][START_REF] Recur | Propagation beam consideration for 3d thz computed tomography[END_REF] . These methods are especially known for their ability to provide a high quality reconstruction even if the acquired data consists of a few number of projections. However they do not take into account a realistic transmission model.

Thus in this paper we propose to investigate the statistical reconstruction method denoted Maximum Likelihood for TRansmission (ML-TR) tomography [START_REF] Fessler | Hybrid poisson/polynomial objective functions for tomographic image reconstruction from transmission scans[END_REF][START_REF] De Man | A study of four minimization approaches for iterative reconstruction in X-ray CT[END_REF]. This method is based on Poisson distribution model of transmitted radiation and allows the introduction of some a priori knowledge about the imaged object [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF][START_REF] Lange | A theoretical study of some maximum likelihood algorithms for emission and transmission tomography[END_REF][START_REF] Lange | Globally convergent algorithms for maximum a posteriori transmission tomography[END_REF]. Specifically, we focus on the implementation denoted the Ordered Subsets Convex (OSC) algorithm since it has an efficient convergence rate despite the noise level and sparsity of acquired data [START_REF] Kamphuis | Accelerated iterative transmission ct reconstruction using an ordered subsets convex algorithm[END_REF][START_REF] Erdogan | Ordered subsets algorithms for transmission tomography[END_REF].

The paper is organized as follows: in the next section, we describe the OSC algorithm and specific optimizations to take into account THz radiation properties. In particular we include the Gaussian beam propagation model proposed in [START_REF] Recur | Propagation beam consideration for 3d thz computed tomography[END_REF][START_REF] Recur | THz radiations for tomographic inspection[END_REF]. Section 3 describes the experimental setup used to perform 3D THz acquisitions and measurements needed by the OSC algorithm. Before concluding, we discuss the reconstruction quality in section 4, by comparing the new algorithm with other standard reconstructions.

Ordered subsets convex algorithm for 3D THz tomography

Expectation maximization (EM) refers to a wide class of iterative reconstructions now well established for X-ray transmission tomography [START_REF] Lange | EM reconstruction algorithms for emission and transmission tomography[END_REF][START_REF] Lange | A theoretical study of some maximum likelihood algorithms for emission and transmission tomography[END_REF][START_REF] Lange | Globally convergent algorithms for maximum a posteriori transmission tomography[END_REF]. These methods reconstruct the 3D structural volume of the sample from a set of radiographs acquired around the object. Assuming Poisson distribution of the acquired data, the probability to observe the measurements R according to the sample μ, denoted p(R|μ), is defined by:

p(R|μ) = ∏ i e -R(i) ( R(i)) R(i) R(i)! ( 1
)
where i is a projection line index and R is the expected intensity measurement, proportional to the number photon counts, at detector position i. R is modeled by the Beer-Lambert law as follows:

R(i) = I 0 (i)e -∑ j w i j μ( j) + bg(i) [START_REF] Shen | Detection and identification of explosives using terahertz pulsed spectroscopic imaging[END_REF] where I 0 corresponds to the maximal photon counts (blank calibration scan), bg is the background noise (dark-field) and w i j is a weight coefficient defining the voxel j contribution on the detector measurement i. w i j is usually proportional to the distance traveled by the line i in the voxel j. I 0 and bg are determined experimentally (cf. section 3.3).

Based on the likelihood p(R|μ) and the transmission model R, the method consists of maximazing the log likelihood log p(R|μ), ie. finding the solution μ that minimizes the partial derivative ∂ log p(R|μ) ∂ μ . The solution can not be determined directly since the partial derivative leads to a transcendantal equation owing to the exponential part in the transmission model. Thus several approximations have been proposed based on series expansion or Newton numerical analysis [START_REF] De Man | A study of four minimization approaches for iterative reconstruction in X-ray CT[END_REF][START_REF] Lange | Globally convergent algorithms for maximum a posteriori transmission tomography[END_REF]. The latter in particular leads to the ordered subsets convex (OSC) algorithm [START_REF] Kamphuis | Accelerated iterative transmission ct reconstruction using an ordered subsets convex algorithm[END_REF]. We define and adapt the OSC algorithm for THz tomography in the following.

Ordered subsets convex (OSC) algorithm

The ordered subsets convex (OSC) algorithm [START_REF] Kamphuis | Accelerated iterative transmission ct reconstruction using an ordered subsets convex algorithm[END_REF] consists of iterating in t and subsets s + 1 in order to update each voxel j of the volume μ until convergence of the solution. The volume obtained by update using the subset s is used as starting volume of the next subset s + 1. A main iteration t is completed when all subsets have been processed. The OSC algorithm updates each voxel as follows:

μ t s+1 ( j) = μ t s ( j) + α t μ t s ( j) ∑ i∈S(s) w i j Rt s (i) -R(i) ∑ i∈S(s) w i j Rt s (i) ∑ k w ik μ t s (k) (3)
where Rt s (i) is the expected photon counts computed from μ t s using Eq. ( 2), α t is a relaxation parameter and S(s) are the radiographs in the subset s. α t = 1 for all t in the following.

THz transmission including Gaussian beam model

In X-ray CT transmission tomography, the beam shape can be considered constant since the wavelength is small with respect to the object size so that there is no diffraction effect. For the same reason, the intensity distribution is uniform over the beam cross-section. This property allows transmission model Eq. ( 2) to be used directly in all reconstruction methods.

In THz CT imaging, the propagation beam is close to a Gaussian distribution (cf. Fig. 1) which is both determined by the THz wave properties and the lens used to focus the beam. The radius of the beam (from the beam axis) has its minimum value w 0 at beam waist. According to the wavelength λ , the radius at the position z from the beam waist is:

w(z) = w 0 1 + z z R 2 ( 4
)
where z R = πw 2 0 λ is the Rayleigh range. Moreover, the intensity distribution over the crosssection is given by:

G( j) = I max w 0 w(z) 2 exp -2r 2 w 2 (z) (5)
where j is the 3D position vector related to the distance r from the beam axis and z (cf.

μ( j) = f ( j) G( j) (6)
where denotes the convolution operator. It results that the volume reconstructed by the OSC is no more μ but f , leading to the following update step:

f t s+1 ( j) = f t s ( j) + f t s ( j) ∑ i∈S(s) w i j Rt s (i) -R(i) ∑ i∈S(s) w i j Rt s (i) ∑ k w ik μ t s (k) (7)
where μ is now defined by Eq. ( 6) in the denominator as well as in the computation of Rt s . The overall algorithm, denoted THz ML-TR in the following, processes as follows: 1) compute the projection values Rt s (i) from current estimate of f for all i in the subset s using Eq. ( 2), where μ is given by Eq. ( 6); 2) update volume f using Eq. ( 7); 3) repeat steps 1-2 for all subsets; 4) iterate steps 1 to 3 until convergence of the solution. In our reconstructions, we have stopped the iterations if, i) the projection quadratic error

∑ i Rt s (i) -R(i) 2 is less than 0.5% of the initial backprojected error ∑ i [R(i)] 2 , or if, ii
) a maximum of 10 iterations have been performed (to avoid long reconstruction time if i) is hard to achieve). We have usually observed that the algorithm stopped between 7 and 9 iterations, with a residual quadratic error less than 0.5 %.

Experimental acquisition setup

We first describe in this section the experimental setup used in order to perform THz tomography. Since the beam propagation follows a Gaussian distribution, we detail how to model this distribution from measurements of the beam waist. Then the tomographic acquisition of the sample is described. Finally, we measure and define the blank scan, I 0 , and dark-field, bg, models used in the THz ML-TR reconstruction.

Description of the experimental setup

The experimental setup of the 3D millimeter wave tomographic scanner is a Gunn diode and a tripler coupled with a horn antenna (cf. Fig. 2). The output power is 12 mW at 287 GHz. The THz beam is then collimated and focused with a pair of PTFE (Polytetrafluoroethylene) lenses ( f = 100 mm focal length and D = 50.8 mm). The sample is positioned on a three-axes motorized stage comprising the X,Y and θ movements, respectively. The detection is performed with a Schottky diode and the beam is modulated at 350 Hz by a mechanical chopper. The amplitude of the transmitted THz signal is acquired with a lock-in amplifier (time-constant: 30ms).

Then we checked the 2D transversal profile of the THz beam at the beam waist and outside the Rayleigh zone. At the sample position, the beam profile is homogeneous with a Gaussian circular shape (2.3 mm beam diameter, measured at full width half maximum (FWHM)) in agreement with the theoretical values obtained from the propagation of Gaussian beam models. Using Eq. ( 5), a simulation of the 3D ray profile can be performed according to the wavelength and the measured w 0 ≈ FW HM at beam waist. Fig. 1(a) shows the 3D profile of simulated 287 GHz source focused with w 0 = 2.3 mm and Fig. 1(b) illustrates the corresponding Gaussian intensity distribution at beam waist. This model is used in the THz ML-TR algorithm to take the Gaussian propagation Eq. ( 5) into account during reconstruction. According to the wavelength and beam waist, the Rayleigh range of this setup, positionated between z 1 and z 2 along the propagation axis such that w(z 1 ) = w(z 2 ) = 2w 0 is about 55.13mm.

Sample acquisition

The sample studied in this paper is a head spray shown on Fig. 3(a). Several regions of this object are out of the Rayleigh range since its largest dimension is greated than 90mm. A 2D transmission radiograph of the sample is obtained by moving the object in the X and Y directions, with a scan step of 1mm in both directions. For 3D tomography, the sample is then rotated in order to provide different radiographs of the object. The overall acquisition of the head spray is processed for 36 projections uniformly distributed around 180 • with an angle step of 5 • and each projection size is 156 × 100 pixels (acquisition time is about 8min per projection). The measured projections are shown on Fig. 3(b) and Media 1.

Measurements of blank scan and dark-field

In order to process the line acquisition Eq. ( 2), the blank scan, I 0 , and dark-field, bg, have to be determined. The blank scan corresponds to the maximal measured intensity when the source is lighting the detector. The background corresponds to the residual ambient noise received by the detector when the source is inactive.

First, we have measured the blank scan for 5 projections (the same size as for sample acquisition) by illuminating the detector with the THz source (without a sample). Thus according to the projection size, we have measured 5 × 156 × 100 = 78000 values, in order to estimate the source intensity distribution. The histogram provided in Fig. 4(a) depicts the blank-scan distribution (Y-Axis) according to the measured intensity (in mA). This distribution can be modeled by a normal distribution N ( Ī0 , σ I 0 ), where Ī0 = 7.0979 and σ I 0 = 0.02056 (cf. green curve on Fig 4(a)). From this first approximation, a fitting algorithm is applied to refine the blank scan distribution parameters. We obtain finally Ī0 = 7.086 and σ I 0 = 0.0165 (cf. blue curve on Fig. 4(a)). Second, we have processed in a similar fashion the dark-field distribution. The measurements (cf. red histogram on Fig. 4(b)) follow a normal distribution approximated by bg = -0.00803 and σ bg = 0.00044. We obtain bg = -0.00780 and σ bg = 0.00035 after

refinement.

Thus efficient models of blank scan and dark-field in Eq. 2 can be the distributions N ( Ī0 , σ I 0 ) and N ( bg, σ bg ) determined experimentally. In that case, each value I 0 (i) and bg(i) is chosen such that a value sets (for all i in one projection) follow these distributions, respectively. However, due to the very low observed standard deviations, we simplify these models by defining I 0 (i) = Ī0 and bg(i) = bg for all i.

Discussion

In this section, we discuss the quality of THz ML-TR reconstruction compared to backprojection of filtered projections (BFP) and the simultaneous algebraic reconstruction technique developed for THz tomography (THz-SART) in [START_REF] Recur | Propagation beam consideration for 3d thz computed tomography[END_REF][START_REF] Recur | THz radiations for tomographic inspection[END_REF]. BFP reconstructs each voxel by averaging projection values crossing it. To assume that a voxel is a linear combination of projection values, these values have to represent a linear combination of volume data. Similarly, SART is based on the Karczmarz algorithm to solve iteratively the linear system of equations between voxels and projections. Then, contrary to the THz ML-TR, these methods are based on the absorbance

A(i) = -log R(i)-bg(i) I 0 (i)
to recover the solution f from all A(i) = ∑ j w i j [ f ( j) G( j)]. Projection preprocessing before BFP or THz-SART is done using the same I 0 (i) = Ī0 and bg(i) = bg values than those used in the THz ML-TR method.

Quality and accuracy observed on reconstructed 2D slices

Fig. 5 shows a set of 6 cross-sectional images along the Y-Axis obtained by the standard BFP, the THz-SART reconstruction and the new THz ML-TR method, respectively. The SART result is obtained after 2 iterations. THz ML-TR convergence is obtained after 10 iterations using 2 subsets. Even if reconstruction time is quite larger, it is negligible compared to the acquisition time.

First we can notice that the reconstructed object is surrounded by many artifacts on BFP results. This method is particularly sensitive to the noise and it is unable to reconstruct accurately from a few number of projections. Indeed, it is common to acquire at least N projections sized N × M to reconstruct accurately a volume of M slices sized N 2 . Thus in our study case, a correct BFP reconstruction of the head spray could be achieved from 156 acquired radiographs. Even if SART reconstruction quality seems better, one can remark the presence of few streak artifacts (high intensity lines). Moreover, several parts such as the base (left image) and outer surface of the acquired object are not reconstructed whereas the intensity of the other ones is very high.

The artifacts owing to both the BFP and SART methods are clearly identified compared to the THz ML-TR reconstruction in the difference images given in Fig. 6(a-b). The transmission model used in the THz ML-TR method leads to a better uniformity of reconstructed intensities: both the inner and outer parts of the object are recovered with high quality. Since transmis- sion value rapidly approaches zero when the beam is going through a high density material, it is close to a nil measurement (nothing detected on the sensor) experienced when the beam is completely absorbed by the sample. Thus THz ML-TR is less sensitive to full absorption effect. Conversely, full absorption provides very high and locally non-proportional values on the absorbance projections. Thus it leads to streak artifacts in the reconstructions, as visible on the BFP, and, to a lesser extent, on the SART results. Moreover, a better accuracy is obtained using the Gaussian beam model in the reconstruction. As an illustration of the improvements, Fig. 6(c) shows ML-TR without Gaussian model, and Fig. 6(d) the difference image between ML-TR and THz ML-TR. We can notice a blur in ML-TR result, even at the centre of the object, owing to the beam propagation considered as linear in the Rayleigh length. Thus optimized reconstruction corrects for blur effect whatever the object position with respect to the Rayleigh range so that Rayleigh length is no longer a limiting factor for 3D THz tomography. As an illustration we can notice an equivalent accuracy of the object shapes which are inside and outside the Rayleigh range in images Fig. 5(c).

Quality and accuracy discussed from 3D rendering images

In order to confirm the new method efficiency compared to the others, we discuss now the feasability of a 3D segmentation in order to provide a 3D rendering of the sample. This study is realized using a new software developed in our group which performs processing sequences to analyze 3D THz images. We will describe the software capabilities in a future paper.

Usually, one applies prefiltering to reduce noise and preprocessing to enhance contrast before 3D rendering. Even if these processings/filterings lead to a better 3D visualisation, they can reduce or erase relevant data about the sample, leading to a biaised analysis of the object. In other words, the inspection performance is directly linked to the tomographic reconstruction quality. Thus in our study, segmentation and rendering are performed directly from reconstructed images (without any preprocessing and prefiltering), allowing us to discuss tomographic reconstruction quality objectively.

The 3D rendering corresponding to the outer sample surface is shown on Fig. 7. Noise of BFP result makes 3D visualisation absolutely wrong and impractical for any outer shape analysis or measures (cf. Fig. 7(a)). 3D rendering of THz-SART volume (cf. Fig. 7(b)) is slightly better since we can distinguish partially the object structure. However, as we have already noticed in Fig. 5, outer surface disappears on the bottom of the object and few streak artifacts penalize 3D visualisation. Conversely, 3D rendering of the THz ML-TR volume ((cf. Fig. 7(c)) allows a 3D visualisation without any noise or streak artifact. From such a 3D rendered volume, ob-

Conclusion

In this paper, the OSC algorithm developed to perform ML-TR reconstruction has been optimized for 3D THz tomography. Optimisations consist of i) including the Gaussian beam model in the computation of expected photon counts, ii) modeling the blank scan and dark-field from measurements and, iii) taking into account these models in the algorithm. The new method, called THz ML-TR, has been compared to BFP and THz-SART reconstructions.

Comparisons have demonstrated that the THz ML-TR method outperforms other reconstructions. Including a more realistic transmission model, it provides a better accuracy of outer and inner structures of the object. 3D rendering, usually processed to perform internal structure inspection, has also demonstrated the quality of THz ML-TR reconstruction. Conversely, other methods are hampered by noise or streak artifacts, making structure inspection hard to achieve.

Future works will focus on the development of THz ML-TR for reconstructing from multienergy and/or multi-intensity acquisitions. We expect that such an approach could solve the hole artifacts observed on the THz ML-TR results. Thus the simplification made here to include dark-field and blank scan models into the reconstruction process will have to be reconsidered and a dual reconstruction update step will be investigated.

Fig. 1 .

 1 Fig. 1. (a) Simulated propagation of the beam for a 287 GHz source (position in meters), w 0 = 2.3 mm at beam waist according to previously measured beam. (b) Gaussian intensity distribution at beam waist (Intensity (a.u)).

  Fig 1(a)), and I max is the intensity at the center of the beam waist. Such a propagation model is included in Eq. (2) by defining:

Fig. 2 .

 2 Fig. 2. Experimental setup: THz beam, delivered by a Gunn diode and a tripler (power 12 mW at 287 GHz) is collimated and focussed with a pair of PTFE lenses. Acquired sample is positionated on a three-axes motorized stage comprising the X,Y and θ movements.

Fig. 3 .Fig. 4 .

 34 Fig. 3. (a) Picture of the acquired head spray. (b) 36 projections of the sample acquired uniformly in the angle range [0 • , 180 • ] with an angle step of 5 • between two consecutive projections (Media 1).

Fig. 5 .

 5 Fig. 5. Few cross-sections along Y-Axis of the imaged object obtained respectively by BFP, THz-SART and THz ML-TR methods (from the first to the last line, respectively). From left to right: y = 15, 30, 45, 60, 75, and 90.

Fig. 6 .

 6 Fig. 6. (a) Difference image BFP at THz ML-TR, and (b) SART and THz ML-TR at position y = 80. (c) Slice at position y = 15 obtained by ML-TR without Gaussian beam model. (d) Difference image between ML-TR and THz ML-TR of cross-sectional images at y = 15.

Acknowledgments

The authors acknowledge the Agence Nationale de la Recherche (ANR) for their financial support in the InPoSec project (www.inposec.org). Thanks are also due to the referees for their insightful comments.


Ahmed Frikha
email: ahmed.frikha@irisa.fr

Samer Lahoud
email: samer.lahoud@irisa.fr

Bernard Cousin
email: bernard.cousin@irisa.fr

Besides

A Hybrid End-to-End QoS Path Computation Algorithm for PCE-based Multi-Domain Networks

Keywords: QoS routing, inter-domain routing, path computation element (PCE), crankback mechanisms, pre-computation, on-demand computation

The Inter-domain Quality of Service (QoS) routing is a challenging problem for today's Internet. This problem requires the computation of paths that cross multiple domains and meet the different QoS constraints. In addition, the methods of computation that are used must meet the constraints of confidentiality and autonomy imposed by the domains of different operators. The path computation element (PCE)-based architecture offers a promising solution for the inter-domain QoS routing. It ensures the computation of endto-end QoS paths while preserving the confidentiality and the autonomy of the domains.

In this paper, we propose a novel hybrid end-to-end QoS path computation algorithm, named HID-MCP, for PCE-based networks. HID-MCP is a hybrid algorithm that combines the advantages of pre-computation and ondemand computation to obtain end-to-end QoS paths. Moreover, it integrates crankback mechanisms for improving the path computation results in a single domain or in multiple domains based on the PCE architecture.

The simulation results show that our algorithm has an acceptance rate of the requests very close to the optimal solution. Moreover HID-MCP outperforms BGP in terms of success rate and the difference is up to 30% in a realistic network. Detailed analysis are provided to assess the performance of our algorithm in terms of success rate and computational complexity.

Introduction

Nowadays, diverse advanced applications are provided over IP-based networks (e.g. IPTV, video-on-demand, and VoIP). Guaranteeing the Quality of Service (QoS) to such applications is a difficult problem, especially when service delivery requires crossing heterogeneous domains under the responsibility of different operators. In such a case, the problem becomes more complex. Moreover, operators adopt different policies which consolidate their economic interests and confidentiality clauses. Thus, cooperation between operators for providing QoS is possible only if it satisfies their policies. Routing is one of the primary mechanisms for providing QoS. It consists in the computation of an end-to-end path which ensures the delivery of the service while meeting the QoS constraints. Several recent works studied the multi-domain QoS routing. Yannuzzi et al [START_REF] Yannuzzi | On the challenges of establishing disjoint QoS IP/MPLS paths across multiple domains[END_REF] presented the challenges of finding disjoint QoS paths in multi-domain networks. One of the most challenging problems is the domain visibility. In fact, information about domain topology and QoS metrics are confidential and cannot be exchanged between domains. This makes finding an inter-domain QoS path a hard task. Some other studies proposed to aggregate QoS information and exchange them between domains. Uludag et al [START_REF] Uludag | Analysis of topology aggregation techniques for QoS routing[END_REF] provided an analysis of the different techniques of topology aggregation for QoS routing.

The path computation element (PCE)-based architecture [START_REF] Farrel | Path Computation Element (PCE)-Based Architecture[END_REF] offers a promising solution for the inter-domain QoS routing. It enables computing end-to-end QoS paths in multi-domain networks while preserving confidentiality clauses and autonomy of the domains by distributing the computations over the domains. This architecture supposes that each domain has one or more PCE responsible of the intra-domain computation and can cooperate with other PCE using the PCE-to-PCE protocol in order to compute an end-to-end path [START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF]- [START_REF] Torab | On cooperative inter-domain path computation[END_REF]. The backward-recursive procedure (BRPC) is introduced by Vasseur et al [START_REF] Vasseur | A Backward-Recursive PCE-Based Computation (BRPC) Procedure to Compute Shortest Constrained Inter-Domain Traffic Engineering Label Switched Paths[END_REF] for the multi-domain path computation. This procedure relies on the PCE architecture to compute an end-to-end QoS path that meets the QoS requirement. This procedure allows PCE to exchange aggregated paths with other PCE using a compact structure, called virtual shortest path tree (VSPT). The drawback of this procedure is that it does not allow domains to exchange all aggregated paths but only one single path per entry border node. Geleji et al [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF] provided a comparison of different schemes for the end-to-end path computation based on the PCE architecture.

In this paper, we propose a novel inter-domain QoS routing algorithm relying on the PCE-based architecture, named HID-MCP. HID-MCP is based on a hybrid computation scheme. The hybrid computation scheme combines the advantages of the pre-computation scheme [START_REF] Orda | QoS routing: the precomputation perspective[END_REF]- [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF] such as the low computational time and the advantages of the on-demand computation such as the high acceptance rate of the requests [START_REF] Ma | Influence of Link State Updating on the Performance and Cost of QoS Routing in an Intranet[END_REF]. Our algorithm consists of two phases: An offline phase and an online phase. In the offline phase, HID-MCP pre-computes a set of QoS paths, as well as look-ahead information for each domain. Look-ahead information gives a measure of the best QoS performance that can be provided by the domain. In the online phase, HID-MCP combines the pre-computed paths to obtain an end-to-end path that fulfills the QoS constraints. Combining the pre-computed paths does not lead always to an end-to-end path. In such a case, a crankback mechanism is executed to perform on demand computations. Combining pre-computation and on-demand computation using a crankback mechanism improves the computation results and allows computational complexity to be reduced.

We distinguish two different crankback mechanisms: intra-domain crankback and inter-domain crankback. The intra-domain crankback is solicited to locally improve the computation results by executing an on-demand computation in the failing domain, while the inter-domain crankback performs a global improvement by executing an on-demand computation starting from one of the downstream domains in the crossed domain set. Note that this solution extends our recently published work in [START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF] by enhancing the algorithm side and simulation side. We enhance the algorithm by adding a new parameter, named hops. This parameter represents the number of backward hops in term of domain to be done before executing the on-demand computation. Precisely, in [START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF], the inter-domain crankback signaling is sent to the destination domain. While in this paper, the inter-domain crankback signaling can be sent to one of the downstream domains in order to execute the on-demand computation starting from this domain, and not only starting from the destination domain. This parameter allows us to tune the performance of the algorithm, and to ensure a trade-off between the rapidity and the success rate of the algorithm. We detail the operation performed by HID-MCP using some examples. We improve the simulation side by considering larger and realistic networks. We also compare the average computational time of the path combination and the on-demand computation in order to show the advantage of using pre-computation. Finally, we study the performance of our algorithm by comparing it with two algorithms. The first algorithm is an exact one, and it represents the optimal solution. The second algorithm is based on the Border Gateway Protocol (BGP), and it represents the inter-domain routing protocol employed in today's Internet. The simulation results show that our algorithm has an acceptance rate of the requests very close to the optimal solution. Moreover HID-MCP outperforms BGP in terms of success rate and the difference is up to 30% in a realistic network. Detailed analysis are provided to assess to the performance of our algorithm in terms of success rate and computational complexity. Besides, our solution relies on the PCE architecture to overcome the limitations related to inter-domain routing such as domain autonomy and confidentiality.

The Inter-Domain QoS Routing Problem

Inter-domain QoS routing, also known as Inter-Domain Multi-Constraint Path (ID-MCP) computation problem, consists in computing a path subject to multiple QoS or cost constraints between a source and a destination node of a multi-domain network. Computing this path requires the knowledge of the QoS metrics on the network links. The QoS metrics can be classified into three types: bottleneck metrics such as bandwidth, additive metrics such as delay, and multiplicative metrics such as loss rate. The multiplicative metrics can be translated into additive metrics using the logarithm function. The bottleneck metrics can be resolved by omitting all links which violate the constraints and then computing the path on the residual graph. Therefore, we consider in the following, only additive metrics.

Let us introduce some notations to formally define the ID-MCP problem. Let G(N, E, D) denote a graph of a multi-domain network, N is the set of nodes, E is the set of links and D is the set of domains. The graph G is composed of D domains. Let m be the number of QoS constraints. An mdimensional weight vector is associated with each link e ∈ E. This vector consists of m non-negative QoS weights w i (e), i = 1..m. Let p be a path in the graph G(N, E, D) and w i (p) be the weight of p corresponding to the metric i. As metrics are additive, w i (p) is given by the sum of the weights of the i th metric of the links of the path p: w i (p) = ej ∈p (w i (e j )). Let The ID-MCP problem belongs to N P-Complete problem [START_REF] Wang | Quality-of-Service Routing for Supporting Multimedia Applications[END_REF] and may have zero, one or multiple solutions (feasible paths). Beside, information about the internal topology or the QoS metrics on the links is confidential since the operators can be in competition. Therefore, computing such a path while meeting all of these constraints is a challenging task.

Currently, the universal inter-domain routing protocol is BGP. Although BGP has proved its efficiency to enable inter-domain routing with confidentiality and autonomy concern, this protocol cannot solve the ID-MCP problem since it does not take into account QoS constraints. One major limitation of BGP is the unique path propagation. BGP cannot propagate several paths for the same destination which provides from QoS-aware inter-domain routing. Another problem is that BGP performs the inter-domain route selection without any coordination between the domains. To overcome these limitations many extensions for BGP are proposed to support QoS routing [START_REF] Knoll | BGP Extended Community Attribute for QoS Marking[END_REF]- [START_REF] Griffin | Interdomain routing through QoS-class planes[END_REF]. However, the QoS capabilities of these propositions remain limited, due to the fact that only one path per destination can be propagated from a domain to another neighbor domain. Therefore, it is not possible to find an end-to-end path that meets the QoS constraints using these propositions.

The research community has recently been exploring the use of distributed architecture to solve the ID-MCP problem. The IETF defines in [START_REF] Farrel | Path Computation Element (PCE)-Based Architecture[END_REF] an alternative architecture based on Path Computation Elements (PCE) that enable the cooperation between domains to find a path subject to multiple constraints. In this architecture, the domains designate one or more PCE, able to perform advanced routing computations within the domain, to cooperate with external PCE and then offer the inter-domain cooperation needed for QoS end-to-end routing. Moreover, the PCE-based architecture preserves confidentiality and the autonomy of each domain by distributing the computation over the domains. Each domain computes locally the segment of the end-to-end path that traverses it.

To our knowledge, few works have been proposed to solve the ID-MCP problem using the PCE architecture. The algorithm proposed in [START_REF] Bertrand | A Distributed Exact Solution to Compute Inter-domain Multi-constrained Paths, The Internet of the Future[END_REF] extends the exact algorithm SAMCRA [START_REF] Van Mieghem | Concepts of exact QoS routing algorithms[END_REF] to an inter-domain level to solve the ID-MCP problem. The drawback of this algorithm is its high complexity. The ID-PPPA algorithm [START_REF] Frikha | Pre-computation Based Heuristic for Inter-Domain QoS Routing[END_REF] and the ID-MEFPA algorithm [START_REF] Frikha | Performance Evaluation of Pre-computation Algorithms for Inter-Domain QoS Routing[END_REF] rely on the PCE architecture and attempt to solve the ID-MCP algorithm using a pre-computation scheme. These solutions are efficient in terms of rapidity and have a good success rate. The pre-computation scheme ensures a very low computational time, but requires a periodical update of the network link state information to maintain a high success rate. Work in [START_REF] Esmaeili | Enhanced Crankback Signaling for Multi-Domain IP/MPLS Networks[END_REF] proposes also a promising distributed solution with crankback mechanisms for inter-domain routing. However this solution cannot take into account several QoS metrics.

3 Computation Schemes for Inter-domain QoS Routing Different computation schemes have been proposed in the literature to solve the QoS routing problem [START_REF] Ma | Influence of Link State Updating on the Performance and Cost of QoS Routing in an Intranet[END_REF]. The on-demand computation scheme attempts to find a feasible path for each request using network state information. The computation is triggered upon the reception of a QoS request. This computation scheme provides a high probability of finding a feasible path since the network state information is up-to-date [START_REF] Ma | Influence of Link State Updating on the Performance and Cost of QoS Routing in an Intranet[END_REF]. This scheme is widely used in existing networks but presents some serious limitations with the emerging applications in the Internet. In fact, the ID-MCP problem belongs to N P-Complete problem, consequently, the performance of on-demand routing in terms of response time is severely affected. In contrast with the on-demand computation scheme, the pre-computation scheme allows the QoS routing problem to be solved while speeding up the response time [START_REF] Orda | QoS routing: the precomputation perspective[END_REF]- [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF]. The pre-computation proceeds in two phases: It prepares in advance a set of paths satisfying predetermined QoS requests. Then, at the reception of a QoS request, it attempts to rapidly provide a feasible path among the pre-computed paths. However, the deployment of a pre-computation scheme in networks can run into problems.

One problem is that the stored paths are computed based on a snapshot of the network state taken before the reception of the QoS request [START_REF] Ma | Influence of Link State Updating on the Performance and Cost of QoS Routing in an Intranet[END_REF]. Consequently, the pre-computed paths do not necessarily satisfy the QoS requests after a change in the network state [START_REF] Ma | Influence of Link State Updating on the Performance and Cost of QoS Routing in an Intranet[END_REF]. In addition, a pre-computation scheme cannot predict all possible QoS requests even if the snapshot remains valid. This decreases the success rate of the pre-computation scheme. For these reasons, a hybrid computation scheme is proposed. This scheme combines the two aforementioned computation schemes and is performed in two phases. As with the pre-computation scheme, the first phase consists in preparing in advance a set of paths or segment of paths. The second phase is more sophisticated. When the pre-computed paths do not lead to a path which satisfies the request, on-demand computation is triggered. Therefore, the hybrid computation scheme ensures a high acceptance rate for requests while reducing the computation time. In this paper, we rely on a hybrid computation scheme to solve the ID-MCP problem.

The HID-MCP Algorithm

In this paper, we propose a novel inter-domain QoS routing algorithm based on a hybrid computation scheme and named HID-MCP (Hybrid ID-MCP). The present section details the operations performed by this algorithm. Fig. 1 illustrates the building block architecture of the algorithm. This architecture is implemented in each PCE of each domain. Based on this architecture, each PCE performs autonomous computations and cooperates with other PCE to compute the end-to-end path. This reinforces domain autonomy and confidentiality. As shown in Fig. 1, the HID-MCP algorithm consists of two phases. In the first phase, named the offline path computation phase, the algorithm executes a path segment computation procedure and computes look-ahead information in each domain. The pre-computed path segments and look-ahead information are stored in a database for later use. The second phase, named online path computation phase, is triggered upon the reception of a QoS request. In this phase, HID-MCP computes an end-to-end path that spans multiple domains and fulfills the QoS constraints. The algorithm attempts to find such a path, first by combining the stored pre-computed intra-domain path segments, second by executing an on-demand path computation procedure that takes benefits from the stored look-ahead information to speed up the computational time of the algorithm.

The Offline Path Computation Phase

The offline computation phase consists of computing in advance a set of intradomain paths subject to multiple predetermined QoS constraints. It also computes look-ahead information at the level of each entry border node of the corresponding domain. In the following, we detail the operations involved in these two computations.

Offline path computation phase

Path segment computation procedure Look-ahead information computation procedure Path combination procedure On-demand path computation procedure

Online path computation phase

Path segments stored Look-ahead information stored

Fig. 1 The building block architecture of the HID-MCP algorithm.

The Path Segment Computation Procedure

This procedure pre-computes a set of paths from each entry border node of the domain toward the other nodes of this domain as well as the entry border nodes of the neighbor domains. These paths satisfy a set of predetermined additive QoS constraints. In practice, some QoS metrics are more critical for certain applications, such as the delay for the VoIP-based applications. Therefore, our procedure pre-computes for each single QoS metric the path which minimizes the weight corresponding to this metric. For example, it pre-computes the path which minimizes the delay; this path can be useful for the VoIP-based applications.

Let D q be the considered domain, n 1 be a border node of D q , n 2 be a node of D q or an entry border node of a neighbor domain, and m be the number of the QoS metrics, our procedure computes m shortest paths from n 1 to n 2 . Each shortest path minimizes a single QoS metric. Hence, from each entry border node n 1 of D q , this procedure computes m shortest path trees. Each shortest path tree is computed using the Dijkstra algorithm and considering a single metric. Therefore, our procedure executes Dijkstra m times per border node. The complexity of the path segment computation procedure is given by equation [START_REF] Yannuzzi | On the challenges of establishing disjoint QoS IP/MPLS paths across multiple domains[END_REF], where B is the number of the entry border nodes of the domain.

O(B * m(N log(N ) + E)) (1)
The complexity of this procedure depends on the number of constraints m. For one border node, this procedure is in O (m(N log(N ) + E)) corresponding to m times the complexity of Dijkstra, which is O((N log(N ) + E)). Considering the B entry border nodes of the domain, the global complexity is then in:

O(B * m(N log(N ) + E)).

Look-Ahead Information Computation Procedure

During the offline phase of HID-MCP, we propose the computation of lookahead information in each domain. This information gives a measure of the best QoS performance that can be provided by the domain. Particularly, it allows the computation search space of a potential on-demand path computation procedure to be reduced. For instance, this information allows infeasible paths to be discarded from the search space of the procedure before exploring these paths. Therefore, look-ahead information reduces the computational complexity of the online phase and contributes to maintain a reasonable response time.

Look-ahead information is inferred from the result of the pre-computation algorithm. Let n 1 be a border node of the domain, and n 2 be a node of the domain or an entry border node of a neighbor domain, and p * n1 →n2;i denotes the pre-computed shortest path between node n 1 and node n 2 considering the metric i. The weight w i (p * n1 →n2;i ) is the lowest possible path weight between n 1 and n 2 . Similarly, let us denote by

→ W * n1 →n2 = (w * 1 , .
., w * m ) the vector where

w * i = w i (p * n1 →n2;i ). Then, → W * n1 →n2
represents the lowest weights to reach n 2 from n 1 for each single metric. We note that a path does not necessarily exist with these lowest weights for all the metrics simultaneously. However, this vector can be used in the online path computation phase to discard infeasible paths from the search space.

The complexity of the look-ahead information computation procedure is given by equation [START_REF] Uludag | Analysis of topology aggregation techniques for QoS routing[END_REF].

O(m * N * B) (2)
Look-ahead information is inferred from the result of the path segment computation. At each entry border node of the domain, there are at most m * N stored pre-computed paths. Hence, at the level of an entry border node n the complexity of computing the N vectors

→ W * n →nj , where n j ∈ N , is in O(m * N ).
Therefore, the complexity of computing the look-ahead information for all the entry border nodes of the domain is in O(m * N * B).

The Online Path Computation Phase

The online path computation consists in finding a feasible end-to-end path using the pre-computed paths and taking advantage of the look-ahead information. Upon the reception of a QoS request, the source and the destination domains are determined. According to the cooperation policy, the service provider computes the best domain sequence that links the source and the destination domain using the PCE-to-PCE protocol [START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF]. Operators can also base their selection for the best domain sequence according to the number of traversed domains, i.e. by selecting the shortest domain sequence. Furthermore, selecting the best domain sequence will reduce the complexity of the problem by limiting the path computation in the domain sequence and not in all domains. The path computation is triggered in the destination domain toward the source domain following the selected domain sequence. Note that, without loss of generality, we rely on the backward-recursive procedure (BRPC), introduced by Vasseur et al [START_REF] Vasseur | A Backward-Recursive PCE-Based Computation (BRPC) Procedure to Compute Shortest Constrained Inter-Domain Traffic Engineering Label Switched Paths[END_REF], for the end-to-end path computation.

Algorithm 1 s

The source node d

The destination node Seq

The selected domain sequence Seq = {D 1 , D 2 , .., Dr} r

The number of domains in Seq D 1

The destination domain Dr

The source domain Hq

The VSPH received in the domain Dq hops

The number of backward hops in terms of domains for the crankback Algorithm 2

Pq

The set of pre-computed paths in domain Dq Eq

The set of egress nodes in domain Dq Iq

The set of ingress nodes in domain Dq m

The number of constraints Algorithm 3 l

The maximum number of shortest paths selected from a VSPH k

The parameter of TAMCRA Ψq

Set of look-ahead information in domain Dq Lq

Set of look-ahead information in domain Dq from Iq to Eq Table 1 The used notations in all the algorithms of the online computation Algorithm 1 Online Phase of HID-MCP (Seq,s,d,hops) if H q+1 = φ then 5:

1: q ← 1; H 1 ← φ;
q ← q + 1; 6:

else if hops == 0 then 7:

H ← On demand computation(Dq, Hq, s, d); 8:

if H = φ then 9:

q ← q + 1; 10: else 11:

reject request ← true; 12:

end if 13: else 14:

H ← φ; q ← max{1, q -hops}; 15:

while (q ≤ r) and not (reject request) do 16:

H q+1 ← On demand computation(Dq, Hq, s, d); 17:

if H q+1 = φ then 18:

q ← q + 1; 1 illustrates the used notations in all the algorithms of the online computation. As shown in this table Seq = {D 1 , D 2 , .., D r } denote the selected domain sequence, where D 1 is the destination domain and D r the source domain. d is the destination node and s is the source node. Algorithm 1 illustrates the operations performed in the online phase of HID-MCP. First, our algorithm attempts to compute an inter-domain path by combining the pre-computed paths in each domain D q in Seq starting from the destination domain D 1 : the path combination procedure is called (line 3). Operations performed by this procedure are detailed in section 4.2.1. The result of the combination procedure in each domain D q is a set of sub-paths linking the destination node to the entry border nodes of the up-stream domain D q+1 . These sub-paths are sent to domain D q+1 to combine them with the pre-computed segments in domain D q+1 . To preserve domain confidentiality, sub-paths are communicated between domain under a novel compact structure named VSPH (Virtual Shortest Path Hierarchy1 ). This structure contains only the end nodes of the paths (the destination node and the entry border nodes of the up-stream domain) as well as the weight vector of each path. The VSPH received in domain D q is denoted by H q in algorithm 1. A virtual path p d→n is represented in Combining the pre-computed paths in each domain can lead to an endto-end path, as detailed in section 4.2.1. However in some cases, no feasible path is found, i.e. the returned VSPH is empty. We introduce in the following two novel approaches using crankback mechanisms in order to overcome this limitation. The first approach executes an intra-domain crankback while the second approach executes an inter-domain crankback. Both of these approaches perform an on-demand path computation. The aim of the on-demand path computation procedure is to provide better results than the pre-computed ones. Operations performed by this procedure are detailed in section 4.2.2.

The parameter hops in the algorithm represents the number of backward hops in terms of domains to be done before starting the on-demand computation. When hops = 1 the on-demand computation should start from the downstream domain D q-1 , and when hops = q -1 the on-demand compuation starts from the destination domain D 1 , etc. For the intra-domain crankback approach hops = 0. The intra-domain crankback approach (lines 6-12) executes the on-demand path computation procedure in the current domain D q , i.e. where the combination has failed. Then, if a feasible path is found in the current domain, this path is sent to the up-stream domain which will resume the path combination procedure. Otherwise, if the algorithm does not find a solution in the current domain, i.e. H q+1 = φ, the request is rejected.

The inter-domain crankback approach (lines 13-23) executes the on-demand path computation procedure starting from the domain D max{1,q-hops} . When qhops < 1, the computation starts from the destination domain D 1 . Each domain executes the on-demand path computation procedure and sends the computed VSPH to the up-stream domain. The computation stops when an end-to-end path is found or when the on-demand path computation procedure does not find a solution, i.e. the returned VSPH is empty. In the latter case, the request is rejected. VSPH: Hq Pre-computed paths Fig. 2 Crankback signaling when the combination fails in domain Dq Fig. 2 illustrates the crankback signaling when the combination fails in the domain D q . In this example, the number of constraints m = 2, and the constraint vector equals [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF][START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF]. Domain D q receives the VSPH H q from domain D q-1 . It combines the pre-computed paths in D q , with the aggregated paths in H q . The combination procedure is detailled in the next section with an example. The result of the combination procedure is the VSPH H q+1 . However, H q+1 does not contain any feasible paths. In this case, the crankback mechanism must be executed. According to the value of hops, the on-demand computation is executed in either the current domain D q (hops = 0) or starting from a downstream domain in {D q-1 , D q-2 , .., D 1 } (hops > 0).

Path Combination Procedure

The aim of this procedure is to combine the paths in the received VSPH with the internally pre-computed one. Algorithm 2 illustrates the operations performed by the path combination procedure. First, the combination procedure selects the pre-computed paths linking nodes in the set I q to nodes in the set E q , where I q is the ingress node set and E q the egress node set (lines 1-13). Then, these paths are combined with the aggregated paths received in the VSPH (line 17). Finally, feasible paths are aggregated and added to the new VSPH which will be sent to the upstream domain. Note that, at the level of the destination domain D 1 there is no received VSPH (H 1 = φ), the procedure selects the feasible pre-computed paths linking the destination d to the entry border nodes of domain D 2 , and aggregates them in a VSPH to be sent to domain The new VSPH: Hq+1 Fig. 3 illustrates an example of path combination with two constraints (m = 2) in an intermediate domain D q . In this example the constraint vector Algorithm 2 Path combination procedure (D q ,H q ,s,d) [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF][START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF]. Domain D q pre-computes two paths corresponding to the weight vectors (5, 2) and (2, 3). The VSPH H q contains two aggregated paths from the destination d to the domain D q . The weight vectors of these aggregated paths are respectively (3, 7) and [START_REF] Torab | On cooperative inter-domain path computation[END_REF][START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF]. We combine the aggregated paths with the pre-computed ones and we obtain two paths corresponding to the weight vectors (8, 9) and [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF][START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF]. These two paths are aggregated in a new VSPH H q+1 and then sent to domain D q+1 .

The complexity of the pre-computed path combination procedure at the level of an intermediate domain D q ∈ {D 2 , .., D r-1 } is given by equation [START_REF] Farrel | Path Computation Element (PCE)-Based Architecture[END_REF], where B max denotes the maximum number of border nodes between two domains.

O(m q * B 2 max )

There are at most m q-1 paths from the destination to each entry border node of the domain D q . In addition, at each entry border node, there are at most m * B max stored pre-computed paths to reach the upstream domain D q+1 . Hence, the complexity of combining the pre-computed paths and the received paths at the level of an entry border node is in O(m q * B max ). This operation is performed at each entry border node between the domain D q and the downstream domain D q-1 . Therefore, the global complexity of this procedure at each domain is in O(m q * B 2 max ).

The On-demand Path Computation Procedure

When the pre-computed path combination procedure does not lead to a feasible path, the on-demand path computation procedure is called in the current domain or starting from the domain D max{1,q-hops} according to the two aforementioned approaches. We propose a modified version of the TAMCRA algorithm to perform the on-demand computation. Work in [START_REF] Van Mieghem | Hop-by-hop quality of service routing[END_REF] shows that TAMCRA is an efficient tunable heuristic for the MCP problem. TAMCRA introduces a new parameter k that limits the maximum number of stored paths at each intermediate node when searching for a feasible path. This parameter allows TAMCRA's performance to be tuned: the success rate can be improved by increasing k at the expense of increased computational complexity. Algorithm 3 illustrates the operations performed by the on-demand path computation procedure. First of all, our proposed procedure computes a prediction for the lowest weight vector to reach domain D q+1 through each path in the received VSPH (lines 11-19). We define for each aggregated path d, n, → W (p d→n ) in the VSPH and for each node n k in E q , a weight vector → W * (d → n j → n k ) that represents the sum of the weight vector of the computed segment p d →nj and the lowest weight vector to reach n k from n j (line 13). Therefore, . This score represents the lowest score to reach d through p d →nj . A path p, that has a score S(p) > 1, is infeasible since it cannot lead to any node in E q while meeting the QoS constraints. Then, we classify the remaining paths in VSPH according to the score S. We select the l shortest paths having the l lowest scores, where l is a parameter of HID-MCP. The parameter l of HID-MCP is very important to reduce the computational complexity of the on-demand computation procedure and to decrease the number of paths exchanged between domains. After that, for each selected shortest path p d →nj , we initialize the node n j by the corresponding weight vectors → W (p d →nj ) and we execute the TAMCRA algorithm starting from node n j to reach the nodes in E q . The parameter l represents the maximum number of executing the TAMCRA algorithm. We note that at the destination domain, i.e. where the computations start, there is no received VSPH. Hence, the on-demand procedure executes Algorithm 3 On demand computation procedure (D q ,H q ,s,d) TAMCRA starting from the destination node. Finally, we aggregate the feasible paths computed by TAMCRA in a new VSPH.

→ W * (d → n j → n k ) = → W (p d →nj )+ → W * nj →n k ,
Fig. 4 illustrates an example of on-demand computation in the intermediate domain D q with m = 2, k = 2 and l = 2. In this example the constraint vector equals → C= [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF][START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF]. The VSPH H q contains four aggregated paths from the destination d to D q . Two aggregated paths from d to n 1 corresponding to the weight vectors [START_REF] Torab | On cooperative inter-domain path computation[END_REF][START_REF] Farrel | Path Computation Element (PCE)-Based Architecture[END_REF] and [START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF][START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF]. And two aggregated paths from d to n 2 with the weight vectors (5, 3) and [START_REF] Farrel | Path Computation Element (PCE)-Based Architecture[END_REF][START_REF] Torab | On cooperative inter-domain path computation[END_REF]. We combine the aggregated paths with The new VSPH: Hq+1 The received VSPH: Hq Look-Ahead information The l shortest paths in Hq (l=2) the look-ahead information in order to discard infeasible paths from H q and select the l shortest paths from H q . The two shortest aggregated paths in H q are p d →n1 and p d →n2 corresponding to the weight vectors [START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF][START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF] and [START_REF] Torab | On cooperative inter-domain path computation[END_REF][START_REF] Farrel | Path Computation Element (PCE)-Based Architecture[END_REF], respectively. Finally, we initialize node n 1 with (4, 4) and node n 2 with (5, 3), and we execute the TAMCRA algorithm twice. One time starting form n 1 and the another time starting from n 2 . We obtain three feasible paths from d to n 5 with the weight vectors: (8, 8), [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF] and [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF]. [START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF] and [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF] are dominated by [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF], so we discard them from the VSPH H q+1 . We also obtain three feasible paths from d to n 6 with the weight vectors: (7, 10), [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF] and [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF][START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF]. [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF][START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF] is dominated by [START_REF] Geleji | A Performance Analysis of InterDomain QoS Routing Schemes Based on Path Computation Elements[END_REF][START_REF] Frikha | Hybrid Inter-Domain QoS Routing With Crankback mechanisms[END_REF], so we discard it from H q+1 .

The complexity of the on-demand path computation procedure at the level of an intermediate domain is given by equation [START_REF] Ash | A path computation element (PCE) communication protocol generic requirements[END_REF].

O(l(k * N log(k * N ) + k 2 * m * E)) (4)
The most significant point that determines the complexity of the on-demand path computation procedure is the number of executions of the TAMCRA algorithm. Knowing that the number of initialized node is less or equal to l, the complexity of this operation is in

O(l(k * N log(k * N ) + k 2 m * E)),
corresponding to l times the complexity of TAMCRA.

Simulation and analysis

In this section, we evaluate the performance of our algorithm HID-MCP in two steps. First, we perform an experimental study of the complexity of the two proposed approaches of HID-MCP: HID-MCP with intra-domain crankback mechanism, and HID-MCP with inter-domain crankback mechanism. We assess the performance of HID-MCP in terms of computational time, and percentage of requests requiring either an intra-domain or an inter-domain crankback mechanism. Second, we compare HID-MCP with some existing well-known inter-domain QoS routing algorithms. We run simulations using two different topologies:

-The first topology is a network of five domains. Each domain is built based on Waxman's model with 50 nodes in each one. The probability that two nodes of the network are connected by an edge is expressed in [START_REF] Calvert | Modelling Internet Topology[END_REF]. -The second topology has a symmetrical backbone and is called SYM-CORE. SYM-CORE contains five interconnected domains and is taken from the work in [START_REF] Guichard | Definitive MPLS Network Designs[END_REF]. This topology is used to assess the algorithms in a realistic case.

For these two topologies, we associate with each link two additive weights generated independently following a uniform distribution [START_REF] Orda | Precomputation Schemes for QoS Routing[END_REF]1023]. The QoS constraints are randomly generated according to the following: Let p 1 and p 2 denote the two shortest paths which minimize the first and the second metrics, respectively. Let Z = [w 1 (p 1 ), w 1 (p 2 )] × [w 2 (p 2 ), w 2 (p 1 )] be the constraint generation space. The problem does not belong to N P-Complete problem outside Z, i.e. either infeasible or trivial. As shown in figure 5, we divide this space into ten zones Z i , i = 1..10 and we browse the space from the strictest constraint zone Z 1 to the loosest constraint zone Z 10 . Then, we assess the performance of the algorithms according to these zones. In the following, each figure measures the variation of a performance metric according to the constraint generation zones Z i , i ∈ 1..10, with a 95% confidence interval. 5.1 Performance evaluation of HID-MCP Firstly, we compare the average computational time of the on-demand computation procedure and the combination procedure. This comparison will show us the importance of the pre-computation phase in reducing the computational time. In Table 2, we see the average computational time of the combination procedure and the on-demand computation procedure per domain according to the constraint generation zones Z 3 , Z 5 , Z 8 and Z 10 , in the SYM-CORE topology. The combination procedure has a compuational time lower than 2 × 10 -3 seconds, while the on-demand compuation has a computational time higher than 1.8 × 10 -1 seconds. For the Waxman's model-based topology, Table 3 illustrates the average computational time of the combination procedure and the on-demand computation procedure in the constraint generation zones: Z 3 , Z 5 , Z 8 and Z 10 . We see that the average computation time of the combination procedure is very low (almost lower than 2.6 × 10 -3 seconds) compared with that of the on-demand computation procedure (almost around 1.3 seconds). It is clear that the combination procedure is very fast compared with the on-demand computation. This is thanks to the pre-computation phase which computes in adavance multiple QoS paths and reduces considerably the compuational time of algorithm.

Secondly, we evaluate the success rate of the combination procedure and the need to execute either an intra-domain crankback or an inter-domain crankback in each domain. This study can help operators to take the best decision to execute either an intra-domain crankback or an inter-domain crankback, whenever the combination procedure fails to find an end-to-end feasible path. It also allows us to deduce the percentage of executions of the on-demand computation procedure in each domain for intra-domain crankback based HID-MCP. To determine if the intra-domain crankback may lead to a feasible path or not, we use the look-ahead information stored in the domain where the combination procedure has failed. Precisely, we combine the received VSPH with the look-ahead information, if all the obtained paths are infeasible, the intra-domain crankback cannot lead to a feasible path, and in such a case the inter-domain crankback is required. Otherwise, an intra-domain crankback may lead to a feasible path in the failed domain.

We first consider the Waxman's model-based topology. In this topology, the domain sequence is composed of five domains. The combination procedure is always successful in the three first domains of the domain sequence: the destination domain (i.e. D 1 ), D 2 and D 3 . The success rate of the combination procedure equals 100%. The percentage of requests requiring either inter-domain crankback or intra-domain crankback equals zero. Therefore, the crankback mechanisms are never called from these domains. Fig. 6, and Fig. 7 illustrate these percentages respectively in the intermediate domain D 4 and in the source domain D 5 . The success rate of the combination procedure in D 4 is still high especially when the QoS constraints are not very strict (from Z 3 to Z 10 ). However, when constraints are strict, the combination can fail. For instance, in D 4 the success rate of the combination equals 38 % in Z 1 . The intra-domain crankback mechanism can lead to a fea- sible path in this domain and in Z 1 with a percentage that equals 40%. For the remaining 22% of the cases, the inter-domain crankback is required.

In the source domain, we note that the success rate of the combination is low when the constraints are strict. Nonetheless, this procedure performs well when the constraints are less strict. In this domain, the percentage of requests requiring the inter-domain crankback equals 100% in Z 1 and Z 2 . While, from Z 4 to Z 10 more than 77% of the requests do not require crankback mechanisms. Now, we consider the SYM-CORE topology. We note that this topology is composed of five domains. The length of the domain sequence is three domains as SYM-CORE has a symmetrical backbone (i.e. all the domains are interconnected through one central domain). In the destination and the intermediate domains of the domain sequence, the percentage of success of the combination procedure equals 100% in all constraint generation zones. For the source domain, Fig. 8 shows the percentage of success of this procedure as well as the percentage of requests requiring the inter-domain crankback, and percentage of requests for which the intra-domain crankback may lead to a feasible path. The success rate of combination is good and increases when constaints are less strict. The inter-domain crankbak is required with a percentage lower than 43% from zone Z 5 to Z 10 . All of these presented results prove that the global empirical complexity of HID-MCP remains reasonable especially when the constraints are not very strict, as the inter-domain crankback is required rarely in these zones.

Finally, we focus on the importance of the look-ahead information to discard infeasible paths before executing the on-demand computation procedure. We define the percentage of infeasible paths (P IP ) as the number of infeasible paths in the VSPH received by a specific domain over the total number of paths in this VSPH. Fig. 9 illustrates the percentage of infeasible paths (P IP ) in the VSPH received respectively, in the source domain (i.e. D 5 ), and in the intermediate domains (i.e. D 2 , D 3 , D 4 ). For the destination domain (i.e. D 1 ), there is no received VSPH as the computation starts from this domain. In Fig. 9, we see that the (P IP ) equals zero in D 2 and is lower than five percent in D 3 . In fact, the weights of the paths in D 2 and D 3 are still low compared to the constraint vector. Therefore, detecting infeasible paths using the look-ahead information is unlikely to happen. However, in the source domain (i.e. D 5 , the last domain in the computation process) the percentage of the discarded infeasible paths is very high, and equals 100% in the constraint generation zones Z 1 and Z 2 . This means that no path in the received VSPH can lead to a feasible path. Therefore, in D 5 when the constraints are generated in Z 1 or Z 2 , there is no need to execute the on-demand computation locally because it cannot find a feasible path. In this case operators have the choice between executing an inter-domain crankback or rejecting the request, depending on their policies and on the importance of request (e.g. priority of the request). In D 4 , the P IP is high and equals 60% in the constraint generation zone Z 1 . Thus, 60% of the paths will be discarded from the VSPH, and this will reduce the computational time of the on-demand computation procedure. These results prove the importance of look-ahead information in discarding infeasible paths, especially in D 4 and D 5 . Nonetheless, the look-ahead information remains important in all crossed domains in order to select the l shortest paths. As explained previously, the parameter l is required to limit the number of executions of the TAMCRA algorithm when calling the on-demand computation procedure. Therefore, the lower is l the faster is the on-demand computation.

Comparaison of HID-MCP with existing approaches

In this section, we evaluate the performance of the HID-MCP algorithm by comparing it with two different kinds of inter-domain routing algorithms: practical algorithms employed in today's Internet and exact QoS algorithms based on the PCE architecture. For the first type of algorithms, we select the Border Gateway Protocol (BGP), as it is the current inter-domain routing protocol employed in the Internet. This protocol has a low computational complexity. However, it does not allow domains to exchange information about the QoS metrics. Thus, computing an end-to-end path that meets the QoS constraints using BGP is not always possible, even if such a path exists. Furthermore, with BGP only one single path per destination can be propagated between two neighbor domains. These two limitations prevent BGP to deal with QoS routing.

The second type of algorithms are the exact ones. An algorithm is called exact if it can find a feasible end-to-end path, when such a path exists. As an exact algorithm, we choose the ID-MCP algorithm introduced by Bertrand et al [START_REF] Bertrand | A Distributed Exact Solution to Compute Inter-domain Multi-constrained Paths, The Internet of the Future[END_REF]. This algorithm has the best success rate as it is exact. However, the complexity of executing ID-MCP in each domain corresponds to the complexity of the SAMCRA algorithm given by: O

((K max * N log(K max * N ) + K 2 max m * E)), where K max = min(exp(N -2)!, m i=1 ci maxj cj ) [21]
. This complexity is very high compared with that of BGP and HID-MCP algorithms.

We define the global success rate (GSR) of the algorithms by the ratio of the number of requests for which a feasible path is found to the total number of requests. For the inter-domain crankback based HID-MCP, we choose the parameter hops = q -1. Thus, the crankback signaling will be sent to the destination domain D 1 . We select this value for hops in order to study the performance of HID-MCP in the two extremity cases: intra-domain crankback and inter-domain crankback starting from the destination domain. For the other values of hops, the GSR of HID-MCP will be lower than that of HID-MCP with hops = q -1 and higher than that of HID-MCP with hops = 0. Fig. 10 illustrates the variation of the global success rate (GSR) of HID-MCP with intra-domain crankback mechanism, HID-MCP with inter-domain crankback mechanism, the exact algorithm ID-MCP, and the BGP algorithm, in the Waxman's model-based topology. The BGP algorithm has the lowest success rate. This is due to the fact that only one path can be propagated from a domain to another neighbor domain using BGP, while with the other algorithms multiple paths can be propagated using the VSPH structure. The ID-MCP algorithm is an exact algorithm, while the others are heuristic ones. This is why the success rate of ID-MCP is the highest. We remark that the success rate of HID-MCP with inter-domain crankback when l = 2 and k = 2 is very close to the success rate of the ID-MCP. As explained in section 4, k is a parameter of TAMCRA and l is the maximum number of paths selected from the VSPH for executing the TAMCRA computation. Of course, the higher l and k, the closer the success rate to the optimal solution given by ID-MCP. Even so, with low values of l and k (l = 2 and k = 2), the success rate of the HID-MCP with inter-domain crankback is close to that of ID-MCP. As expected, the success rate of HID-MCP with intra-domain crankback when l = 2 and k = 2 is lower than that of HID-MCP with inter-domain crankback with the same parameters, especially in the middle constraint generation space. In fact, when the combination procedure fails, HID-MCP with inter-domain crankback executes the on-demand computation procedure starting from the destination domain, while HID-MCP with intra-domain crankback executes it only in the current domain. Consequently, the quality of the paths computed by the inter-domain crankback approach in each domain is better than the ones computed by the intra-domain crankback approach. Thus, the probability that a feasible path is found using HID-MCP with inter-domain crankback is higher. However, its computational complexity is high compared to HID-MCP with intra-domain crankback, but remains acceptable compared to ID-MCP. Finally, we evaluate the success rate of algorithms in the SYM-CORE topology. In Fig. 11, we see that the inter-domain crankback based HID-MCP, the intra-domain crankback based HID-MCP, and the ID-MCP algorithm have almost the same success rate, while BGP has the lowest success rate. The differ-ence of success rate between HID-MCP and BGP equals 30% in the constraint generation zone Z 10 . The fundamental result deduced from these figures is that HID-MCP has a global success rate very high compared to that of BGP and very close to that of the exact algorithm.

Conclusion

In this paper, we studied the inter-domain QoS routing problem. We proposed a novel inter-domain QoS routing algorithm based on a hybrid computation scheme, named HID-MCP. Our algorithm relies on the PCE architecture. This architecture allows end-to-end path computation subject to multiple QoS constraints, while preserving the confidentiality clauses and the autonomy of the domains. We introduced two different mechanisms for improving the success rate of the HID-MCP algorithm. The first mechanism performs local improvement using an intra-domain crankback, while the second one executes a global improvement using an inter-domain crankback.

The simulation results showed that HID-MCP, with both of the aforementioned approaches, has a success rate very close to the optimal solution. Moreover, HID-MCP has a success rate very high compared with that of BGP and the gap is up to 30% in the realistic network SYMCORE. We also compared the average computational time of the path combination and the on-demand computation. This comparison study showed us the advantage of using precomputation in reducing the computational time of our proposed algorithm.

As future prospects of this work, the stability of the pre-computed paths should be analyzed when they deal with dynamic changes of link state. This allows inferring the validity of the pre-computed paths after an eventual change of the network conditions.

→W 1

 1 (p) = (w 1 (p), w 2 (p), ..., w m (p)) denote the weight vector of the path p. Definition Given a source node s, a destination node d and a set of constraints given by the constraint vector → C= (c 1 , c 2 , .., c m ), the Inter-Domain Multi-Constraint Path (ID-MCP) computation problem consists in finding a path p which satisfies w i (p) ≤ c i , ∀i ∈ 1..m. Such a path p is called a feasible path.

  the VSPH by d, n, → W (p d→n ) , where d is the destination node, n is an entry border node of the upstream domain D q+1 , and → W (p d→n ) is the weight vector of p d→n .

Fig. 3

 3 Fig. 3 Combining the pre-computed paths in domain Dq with the VSPH.

Fig. 4

 4 Fig. 4 The on-demand computation procedure in domain Dq with m = 2, k = 2 and l = 2

Fig. 5

 5 Fig. 5 Constraint generation zones for m = 2

Fig. 6

 6 Fig. 6 Domain D 4 : Success rate of the combination procedure, percentage of requests requiring intra-domain crankback, and percentage of requests requiring the inter-domain crankback.

Fig. 7

 7 Fig. 7 Domain D 5 : Success rate of the combination procedure, percentage of requests requiring intra-domain crankback, and percentage of requests requiring the inter-domain crankback.

Fig. 8

 8 Fig. 8 Domain D 3 : Success rate of the combination procedure, and percentage of requests requiring intra-domain, and percentage of requests requiring the inter-domain crankback.

Fig. 9

 9 Fig. 9 Percentage of infeasible paths in the received VSPH in the Waxman's model-based topology.

Fig. 10

 10 Fig. 10 Comparison of the global success rate of the algorithms in the Waxman's modelbased topology.

Fig. 11

 11 Fig. 11 Comparison of the global success rate of the algorithms in the SYM-CORE topology.

Table

	19:	else
	20:	reject request ← true;
	21:	end if
	22:	end while
	23:	end if
	24: end while
	25: Return reject request == f alse

1 :

 1 Pq ← {p/p pre-computed path in domain Dq}; 2: H q+1 ← φ; 3: if Dq == D 1 then Selected paths ← pn j →n k /pn j →n k ∈ Pq, n j ∈ Iq, n k ∈ Eq ; 14: if Dq = D 1 then 15:for p d→n j ∈ Hq do 16:for pn j →n k ∈ Selected paths do

	4:	Iq ← {d};
	5: else
	6:	Iq ← {n j /n j leaf node in Hq};
	7: end if
	8: if Dq == Dr then
	9:	Eq ← {s};
	10: else
	11:	Eq ← {n k /n k entry border node of domain D q+1 };
	12: end if
	13: 17:		→ W (pn j →n k );
	18:	if p d→n k is feasible then
			→
	19:	Add d, n k ,	W (p d→n k ) to H q+1 ;
	20:	end if
	21:	end for
	22:	end for
	23: else
	24: 25:	for p d→n k ∈ Selected paths do if p d→n k is feasible then
	26:	Add d, n k ,
	27:	end if
	28:	end for
	29: end if
	30: Return H q+1 ;
	equals

→ W (p d→n k ) ← → W (p d→n j )+ → W (p d→n k ) to H q+1 ; → C=

  look-ahead information in domain Dq ; 3: if Dq = D 1 then 4:Iq ← {n j /n j leaf node in Hq}; n j →n k ∈ Ψ, n j ∈ Iq, n k ∈ Eq ;

	1: temp paths ← φ; f easible paths ← φ;
	2: Ψq ←
	5:	if Dq == Dr then
	6:	Eq ← {s};
	7:	else
	8:	Eq ← {n k /n k entry border node of domain D q+1 };
	9:	end if
	10:	Lq ←
	11:	for d, n j ,
				→
				W * n j →n k ;
	14:	end for
	15:	S(p d→n j ) ← min n k ∈E max i∈1..m	w * i (d →n j →n k ) c i	;
	16:	if S(p d→n j ) ≤ 1 then
			→
	17: 18:	Add [n j , end if	W (p d→n j ), S(p d→n j )] to temp paths;
	19:	end for
	20:	if temp paths = φ then
	21:	Selected paths ← l shortest paths having the lowest S in temp paths
	22:	else
	23:	H q+1 ← φ;
	24:	Return H q+1
	25:	end if
	26:	for
		D q+1
	29:	Add the obtained feasible paths to f easible paths
	30:	end for
	31: else
	32:	Execute TAMCRA in D 1 starting from d
	33:	Add the obtained feasible paths to f easible paths
	34: end if
	35: Extract H q+1 from f easible paths
	36: Return H q+1

→ W * / → W * → W * → W (p d→n j ) ∈ Hq do 12: for → W * n j →n k ∈ Lq do 13: → W * (d → n j → n k ) ← → W (p d→n j )+ → W (p d→n j ) ∈ Selected paths do

27:

Initialize n j with the weight vector → W (p d→n j ) 28:

Execute TAMCRA in Dq starting from n j toward every border node of domain

Table 2

 2 Comparison of the average computational time of combination procedure and the on-demand computation procedure in the SYM-CORE topology.

	Constraint zones	Z 3	Z 5	Z 8	Z 10
	Combination Procedure	4.26 × 10 -3 s	4.37 × 10 -3 s	4.7 × 10 -3 s	5.81 × 10 -3 s
	On-demand Procedure	1.73 × 10 -1 s	2.09 × 10 -1 s	2.34 × 10 -1 s	2.97 × 10 -1 s

Table 3

 3 Comparison of the average computational time of combination procedure and the on-demand computation procedure in the Waxman's model-based topology.

The hierarchy is a structure which enables the storage of multiple paths between any two nodes[START_REF] Molnar | Hierarchies for Constrained Partial Spanning Problems in Graphs[END_REF].