Shohreh Ahvar
email: shohreh.ahvar@telecom-sudparis.eu

Jagruti Sahoo
email: jagrutiss@gmail.com

Ehsan Ahvar
email: ehsan.ahvar@inria.fr

Mouhamad Dieye
email: dieye.mouhamad@gmail.com

Roch Glitho
email: glitho@ece.concordia.ca

Halima Elbiaze
email: elbiaze.halima@uqam.ca

Noel Crespi
email: noel.crespi@mines-telecom.fr

PCPV: Pattern-based Cost-efficient Proactive VNF placement and chaining for value-added services in content delivery networks

Keywords: Content Delivery Networks, Network Function Virtualization, Virtual Network Function, Cost, Placement

come

I. INTRODUCTION

Content Delivery Networks (CDNs) are largely distributed infrastructures of surrogate servers (hereafter reffered to as server) placed in strategic locations [START_REF] Pathan | Content Delivery Networks: State of the Art, Insights, and Imperatives[END_REF] to serve end-users (EUs) with reduced latency. Value-added services (VASs) now play a very important role in CDNs. Now around 47% of Akamais revenues come from its VAS offerings [START_REF]Akamai q2 earnings preview: Value-added services will drive growth, but currency headwinds likely played the minor spoilsport[END_REF]. Some examples are website/application acceleration (e.g., route optimization, TCP optimization, stream splitting) [START_REF] Pathan | Advanced Content Delivery, Streaming, and Cloud Services[END_REF], analytics, content protection, advertisement overlays/tickers and content adaptation (e.g., transcoding, compression).

Network Function Virtualization (NFV) [START_REF] Mijumbi | Network function virtualization: State-of-the-art and research challenges[END_REF] could enable CDNs to provision VAS as chains of Virtual Network Functions (VNFs), which allows VASs to scale in an elastic manner. Besides, the updating of an existing VAS or the introduction of new VASs could be achieved with increased agility. We model the Placement and Chaining of VNFs (PC-VNF) for VAS as an optimization problem where the goal is to find the optimal number, location and efficient chaining of VNFs instances so that the CDN provider cost is minimized and QoS is satisfied. This paper focuses on a proactive placement of VNFs where VNFs are deployed in an optimal manner before any request is received from the EUs to access the service. This type of deployment is triggered when a content provider requests the CDN provider to deploy a set of VASs. The rest of the paper is organized as follows: Section 2 presents PC-VNF problem in CDNs and ILP formulation. Section 3 describes the proposed algorithm. Section 4 portrays the simulation results. Section 5 concludes the paper.

II. VNF PLACEMENT PROBLEM

A. Problem Description

Assume that content providers request the CDN provider to deploy a set of VASs. The request specifies the servicerelated parameters (e.g., the description of functionalities that constitute the services) including the QoS threshold to be satisfied. Given a content X, the servers containing X, a set of EUs requesting the content, their workload and a set of services to be accessed by the EUs, the PC-VNF problem targeted in this paper is to find the number and location of VNF instances along with assigning the Service Function Chains (SCs) to the EUs for minimizing the cost of CDN provider and satisfying the QoS threshold of all services offered to the EUs. The reader should note that for each EU, the delay for viewing the video is the sum of the following delays: 1. Delay from the server (selected to stream the video) to the first VNF of SC and; 2. Delay for transmission and processing the video in SC and; 3. Delay from the last VNF to the EU. Our problem is a slight variation of the well-known Bin Packing problem [START_REF] Monaci | A set-covering-based heuristic approach for binpacking problems[END_REF] and the Hierarchical Facility Location-Allocation problem [START_REF] Farahani | Survey: Hierarchical facility location problem: Models, classifications, techniques, and applications[END_REF] which is an NP-Hard problem, calling for an efficient heuristic

B. ILP Formulation

Let us consider N as a set of servers and U as a set of EUs. The physical topology of the network is represented by a directed graph G = (V, E), where V = N ∪ U is the set of nodes composed of the servers and EUs connected by

h f = G(V f NF , E f NF )
, where V f NF is the set of VNFs that will be installed on nodes in N. E f NF is the set of virtual edges that dictate the head and tail of SC. Table I delineates the inputs and variables used in our ILP formulation. Our objective is to minimize the cost shown in (1) including software license cost per instance, α k and site γ, the sum of operational costs for all deployed VNF instances, the communication cost (i.e., sum of the bandwidth costs amongst each pair of servers u, v ∈ N hosting VNFs of SC in each EU service request

h f ∈ H, f ∈ U, β f (u,v) ).
Furthermore, it includes the bandwidth costs between the servers u ∈ N hosting the tail VNFs of SC in each EU service request h f ∈ H, and EU f ∈ U,

β f (u,v) . Where, β f (u,v) = L f • σ (u,v) • B u is the cost of using bandwidth between u ∈ V and v ∈ V, with hop count σ (u,v) , for load L f for EU f ∈ U. min ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ k ∈K,n∈N, j ∈I k x k,n, j × α k + f ∈U,(p,q)∈E f NF ,(u,v)∈E β f (u,v) • y f , p,q u,v + f ∈U,k ∈V f NF |∀i ∈V f NF (k,i)∈E f NF ,u ∈N, j ∈I k ,(u, f )∈E β f (u, f ) • λ f k,n, j + z n × n∈N (γ + k ∈K, j ∈I k (R k • δ n • x k,n, j )) (1)
VNF Placement:

n∈N j ∈I k λ f k,n, j = 1, ∀ f ∈ U, k ∈ K (2) f ∈U L f • λ f k,n, j P k • x k,n, j ∀k ∈ K, n ∈ N, j ∈ I k (3) k ∈K j ∈I k R k • x k,n, j C n • z n ∀n ∈ N (4)
VNFs chain mapping:

j ∈I p λ f k,u, j • j ∈I q λ f k,v, j = y f , p,q u,v ∀ f ∈ U, (p, q) ∈ E f NF , (u, v) ∈ E (5) (u,v)∈E y f , p,q u,v = 1 ∀ f ∈ U, (p, q) ∈ E f NF (6) f ∈U (p,q)∈E f NF L f • y f , p,q u,v BW (u,v) ∀(u, v) ∈ E (7)
QoS guarantee:

(p,q)∈E f ,(u,v)∈E,∀ f ∈U D (u,v) • L f • y f , p,q u,v + k ∈V f ∀i ∈V f ∃(k,i)∈E f ,u ∈N, j ∈I k ,∀ f ∈U T k,u • L f • λ f k,u, j + k ∈V f ∀i ∈V f (k,i)∈E f , u ∈N, j ∈I k ,(u, f )∈E,∀ f ∈U (T k,u + D (u, f ) ) • L f • λ f k,u, j D h f (8) Variable x k,n, j ∀k ∈ K, n ∈ N, j ∈ I k is used to identify unique instance j ∈ I k of VNF type k ∈ K on server n ∈ N. Variable z n ∀n ∈ N is used to record servers hosting VNFs.
We ensure in (2) that for each EU service request h f ∈ H, f ∈ U, an instance of the requested VNF type must be assigned to a server. However, an instance of a VNF type can cater to multiple users. ( 3), ( 4) and ( 6) ensure that the VNF instances, the host servers and links are not overloaded respectively. The edge (p, q) between two consecutive VNFs in each EU service request must be assigned to a physical edge (u, v) between two servers u and v, in [START_REF] Monaci | A set-covering-based heuristic approach for binpacking problems[END_REF]. As ensured in [START_REF] Farahani | Survey: Hierarchical facility location problem: Models, classifications, techniques, and applications[END_REF], the VNFs and their respective ordered edges are mapped to only one pair of physical servers and their edge. ( 8) is QoS constraint which includes the VNF processing delay on the surrogates, the network communication delay between the VNFs and the delay from the EUs to the last VNFs of SCs.

III. PATTERN-BASED COST-EFFICIENT PROACTIVE VNF PLACEMENT (PCPV)

PCPV is run once for each SC which is accessed by a set of EUs. VNF placements should satisfy QoS and prevent VNF and server from overloading. PCPV has two initial placement and modification phases. Phase 1 simplifies PC-VNF by focusing on QoS satisfaction and relaxing the VNF overloading constraint. To simplify the PC-VNF, two assumptions are considered: (i) users can be located everywhere in the network (ii) VNFs do not have capacity limitations. Considering the simplified problem, phase 1 finds the minimum number of instances of each VNF to cover the entire network with promised service delay, selects appropriate servers to deploy the instances and chains the deployed VNFs. Phase 2, if necessary, modifies the initial PC-VNF (i.e., adding/removing the selected servers and VNFs) based on extra information such as the VNF's capacity, the predicted EU's load and location information to prevent VNF (and server) from overloading.

A. Phase 1: Initial VNF Placement and chaining

To initially place the VNFs, we propose partitioning and patterning techniques. Due to the limited capacity of servers, it may possible a SC cannot be allocated on a server. Therefore, the phase 1 first divides SC to some partitions (or PTs) composed of one or more VNFs and is located on one of the servers. Phase 1 then finds the minimum number of needed instances of each PT to cover the network with promised QoS and selects appropriate servers for placing these PT instances on them. We do both of them at the same time using patterning technique. To this end, for each PT, phase 1 creates a pattern.

By mapping the patterns into the network, the minimum number of needed instances to cover the whole network is found and also location (servers) of these PT instances is estimated at the same time. This pattern-based method provides a scalable method for any size of network where the source of the chain is unknown in advance. The five parts of phase 1 are described as follows.

1) Chain Partitioning: PCPV divides a chain h into several PTs, PartNum. The maximum possible size of PTs equals to the minimum capacity, C, of the servers. To accomplish partitioning, PCPV creates and fills the first PT. If still some VNFs of the chain remain, PCPV fills the second PT. This procedure continues until all VNFs of the chain are selected.

In case of VNF fragmentation during filling a PT, the fragmented VNF is migrated to the next PT. After partitioning, the capacity of the smallest VNF of PT i is defined as PartCap[i] and the total requirement of the VNFs of PT is considered as PT requirement Re [i]. Also, PTs numbered in ascending order from the chain head to the EU. Each PT i receives the content as input and yields the processed output content to the next PT i+1 in SC. The output of the last PT in SC is delivered to the EU. Thus, the EUs are directly connected to the server hosting the last PT of SC. We define the ones (i.e., EUs or another PTs of the chain) who receive the content from an instance of the PT as the customers of that PT instance.

2) Pattern designing: This section first presents the proposed patterns. Then, it describes why and how a pattern size is adapted based on the service characteristics (e.g., the maximum Delay Threshold (DT) and the number of created PTs (chain length)).

Proposed patterns-After finding the number of PTs for SC, PCPV defines a pattern for each created PT. In order to design a pattern, for the sake of QoS fairness, the best idea is to deploy an instance in the center of an area to reduce the number of needed instances for satisfying a requested QoS. However, finding a server exactly in the center of an area is not always possible (due to equipment limitations). To this end, PCPV considers a square shape zone (instead of a point) in the center of the area to deploy an instance (second part of this section explains how we set the zone size to make sure at least one server located inside it). For the sake of unity, we design our patterns based on this square shape zone unit with edge of d. PCPV also takes into account location of PT customers to define pattern for each PT. We present our propose patterns with a simple example. Assume a chain including three PTs. Where the EUs are customers of PT 3 instances, PT 3 instances themselves are customers of PT 2 instances and, finally, PT 2 instances are customers of PT 1 instances. In order to design a pattern for each PT, we first need to know location of its customers. As the locations of PTs 1 and 2 customers are not clear prior to placement of PT 3 instances, we start from the last PT where its customers are EUs.

For the last PT, as EUs (i.e., customers of the last PT) are assumed to be located anywhere in the network in the phase 1, we propose a pattern to cover anywhere in the network with the minimum number of needed last PT instances. To this end, we propose a nine same-size zone pattern for the last PT to find the central zone in the area. As it is a symmetric pattern, deploying the last PT instance in a server in central zone minimizes the variance of the distance (i.e., delay) between the last PT and the EUs located in any place of the nine-zone pattern. Therefore, a fairness in QoS is provided for the PT instance customers. Here, the maximum distance between each instance of PT 3 and its customers,

X 3 , is 2 √ 2.
Likewise, pattern 2 is created based on merging four copies of pattern 3 together. An instance of PT 2 is placed in a zone in the center of the merged squared-shape area (i.e., the sticking point of these four copies of pattern 3) to offer a minimized variance of delay to its customers (i.e., instances of PT 1). The maximum distance between each instance of PT 2 and its customers, X 2 , is 2.5 √ 2. We have the same situation for PT 3 to offer a minimized variance of the delay to its customers (i.e., instances of PT 2). The maximum distance between each instance of PT 1 and its customers, X 1 , is 4 √ 2. This procedure can be continued in situation of exiting more PTs. Finally, the maximum chain length can be calculated based on the three defined patterns (i.e., MaxChainLength=

X 3 + X 2 + X 1 =8.5 √ 2d).
We can generalize the example as follows: for a PT i of a SC with PartNum PTs, its pattern size is 4 times larger than pattern size of PT i+1 (i.e., its edge is two times longer). Therefore, as edge size of the nine-zone pattern of the last PT is 3 * d, PCPV sets a general formula to find edge PW i (in unit of zone with edge d) for PT i pattern using equation [START_REF] Chew | Finding Largest Empty Circles with Location Constraints[END_REF]:

PW i = (3 × d) × 2 PartNum-i (9)
Instance of pattern i is then deployed in the central zone of the pattern. Finally, based on the defined patterns, for a SC with PartNum PTs, its maximum length is calculated in terms of zone size d by 10:

MaxChainLength = Part Num i=1 X i = k Part Num • d (10)
Where X i is the maximum possible distance between PT i instances (anywhere in the central zone) and their customers (anywhere in the pattern). As the chain length is in the unit of d, we can also consider it as coefficient of d. Next part describes how d is calculated.

After defining the patterns for PTs of a SC, PCPV identifies the appropriate pattern size for placement to minimize the number of used servers while satisfying the service DT. To identify optimal patterns size, as patterns are designed based on the same-size zones, we identify the optimal zone size. Zone size has an important effect on both service delay and cost. Increasing size of the zone decreases the operational cost, (due to deploying fewer instances and using fewer servers) and increases both communication cost and service delay (due to increasing the chain length). As the operational cost plays greater role on total cost in comparison to the communication cost, PCPV attempts to increase zone size.

However, PCPV cannot increase size of the zone (which leads to increasing chain length) in a way that violate the requested DT. Hence, PCPV increases the zone size (which leads to increasing the chain length) up to the point that still satisfying requested DT (and communication cost). To find the optimal zone size, PCPV calculates the maximum permitted SC length based on the service DT using (11), inspired by [START_REF] Wehrle | Modeling and Tools for Network Simulation[END_REF].

dist = (φ -δ) • W 1000 (11)
Where dist is the distance (in mi) between two nodes, W is the propagation speed of the network medium (in mi sec ) and φ is the end to end delay (in ms). We reserve ( 1k ) of total end to end delay for delay from content to the first VNF where k is the number of VNFs in the chain and we call it δ. Using (11), the maximum chain length (dist) is calculated by considering φ as the service DT. Finally, putting the maximum permitted SC length calculated from (11) in equation ( 10), PCPV finds the optimal zone size d opt .

Nevertheless, the low values of service delay (φ) may result in a relatively small optimal zone size. Hence, some of the selected zones may have no server to be selected. To overcome this issue, we use the algorithm proposed in [START_REF] Toussaint | Computing largest empty circles with location constraints[END_REF] for solving the Largest Empty Circle (LEC) problem [START_REF] Chew | Finding Largest Empty Circles with Location Constraints[END_REF]. It computes the largest empty circle (d 0 ) such that in any other circle with a greater diameter in the same topology, we will have at least one node. We use this notion to ensure that there will be at least one server inside a zone. The final optimal zone size is therefore determined as the maximum of d opt and d 0 .

3) Patterns mapping: By mapping the defined patterns (as described in section III-A2) to the network, a set of zones are determined,

Z tot Sel = Z i Sel 1 i PartNum . Where, Z i
Sel is a set of selected zones for instances of PT i, Z i Sel = Z i j 1 j |Z i Sel | . Z i j refers to the selected zone for the jth instance of PT i.

Two special cases may arise. The first one consists of having only one copy of a PT in the network. Thus, the succeeding PT is placed in the location closest to the previous one. In the second case, the pattern size is greater than the network size. Accordingly, some patterns cannot be entirely located in the network. Thus, if a server of a pattern is located outside the network, the pattern is shifted into the network until a server is found. 4) Server Selection: For each selected zone Z i j , a server with the least cost (e.g., most energy efficient) is selected to place the instance j of PT i. In this case, S i Sel is a set of selected servers for instances of PT i, S i Sel = S i j 1 j |S i Sel | . We consider S tot Sel = S i Sel 1 i PartNum as a set of all selected servers for all PTs. As a result, for each PT i, a set of instances

I i dep = I i j 1 j |I i dep |
, are deployed in the network (I i j refers to the j th instance of PT i). 5) Chaining: For each instance (e.g., I i j ), its customers in its cover area A i j are assigned to it to receive content. A i j is the indicated area for I i j defined by the pattern of PT i. However, EU assignment to last PT will be done in phase 2 (remind that we do not have EUs location information in phase 1).

B. Phase 2: VNF Placement Modification

Phase 2 includes two following steps.

1) Placement Modification:

As phase 2 has extra information of EUs location, it first assigns the EUs to the instances of the last PT based on the cover area of last PT instances.

Afterwards, phase 2 can calculate load of all instances based on the assigned costumers to them, where ρ(I i j ) is load of instance I i j . Phase 2 then removes all instances without customers and it finds overloaded instances (i.e., instances with overloaded VNFs) starting from the first instance of the last PT to the last instance of the first PT of the chain. When the load of an instance I i j (i.e., ρ(I i j )) exceeds its capacity (i.e., PartCap[i]), a new instance of PT i has to be deployed.

If there is enough capacity on the server S i j where instance I i j is deployed, the new instance is deployed on that server as well. Otherwise, we select the next lowest cost server (called S i new ) in zone Z i j to deploy the new instance. If no server in zone Z i j with enough capacity is available, we expand our search area to the closest zone of the intersection region of A i j and A i-1 j (i.e., CL i j ). j is the index of PT i -1 instance which I i j is already assigned to it as its customer. Deploying the new instance inside CL i j guarantees that the distance between the new instance of PT i and instance I i-1 j does not violate the maximum possible distance (X i-1 ) and service DT. However, the new instance is not located in the center of the A i j . As the new instance is closer to instance I i-1 j than to already deployed I i j , SC length becomes shorter. The reduction of SC gives possibility of supporting customers farther than predefined cover area of the new instance in phase 1 and, therefore, can cover whole A i j . Considering the mentioned example in section III-A2, assume instance I 3 1 is overloaded and there is no server with enough capacity in zone Z 3 1 where instance I 3 1 has been already deployed. As I 3 1 has been assigned to I 2 1 , PCPV places the new instance of PT 3 in the closest zone to instance I 2 1 . While the distance between instance I 3 1 and instance I 2 1 is 2.5 √ 2d, the distance between the new deployed instance of PT 3 and I 2 1 is reduced to 1.5 √ 2d. This reduction allows the new instance of PT 3 to cover its customers located farther than the pre-defined distance of 2 √ 2d (up to 3 √ 2d).

It is noted that, after the placement modification, the initial number of instances for each PT obtained in phase 1 based on the pattern mapping, may be changed.

2) Customer Assignment: After modification, if necessary, phase 2 should chain/re-chain a PT instance to a next-level PT instance and/or assign EUs to the instances of the last PT.

IV. PERFORMANCE EVALUATION

As processing delay and link bandwidths is not considered in the simulations, PCPV is compared with the modified version of ILP implemented in CPLEX. In particular, the constraint ( 7) is relaxed and in constraint [START_REF] Toussaint | Computing largest empty circles with location constraints[END_REF], the processing delay of VNFs is ignored and the links delay is traffic independence. We run simulations for 5 different random mesh topologies and report the average of 5 simulations. Three scenarios are designed for small and large environments.

A. Experimental Set up

Scenario I compares PCPV to ILP. A CDN provider spans on 300*300 mi 2 geographical area including nine states, each state of size 100*100 mi 2 . We consider that the CDN provider owns a server in each state and the servers are placed randomly in the states. For this configuration, we consider 9, 12, 15 and 18 EUs number and 3 VNFs per SC with the Service Delay Threshold (D th h ) of 1.5 ms. Scenario II shows the scalability of PCPV. Twenty various configurations including a different number of states and EUs are defined. We simulate different CDN networks environments with the size of 600*600 mi 2 including 36 states, 700*700 mi 2 including 49 states, 800*800 mi 2 with 64 states, 1200*1200 mi 2 including 144 states finally, an environment with a size of 2500*2500 mi 2 including 625 states. Four various numbers of EUs including 50, 100, 150 and 200 are considered. Scenario III shows the impact of the number of VNFs in a SC. We investigate the performance of the proposed algorithm by considering SCs that consist of 3,6 and 9 VNFs in the environment including 625 states. We also assume that the CDN provider owns at least one server in each state of size 100*100 mi 2 . The number of servers in each state is chosen randomly between 1 and 5. Bandwidth Cost (B), 10 Dollar/Gbps, EU Load (L f ), 1 Gbps, servers capacity (C n ) 32 vCPU, site license cost (γ) 1000 Dollar, servers Operational cost (δ n ) randomly between 5 and 10 Dollar/vCPU and VNF license cost (α k ), 1000 Dollar/vCPU are the other set up setting in this experiment.

B. Results and Discussions

To evaluate the effectiveness of PCPV, execution time and cost components are calculated.

As mentioned in section III-A, PCPV focuses on reducing operational cost to achieve total cost reduction. However, our method is designed for provisioning the services with the QoS constraint which needs to be satisfied. Placing SCs in such a manner that satisfies the service DT for the EUs leads to keep the communication cost lower than a threshold. Therefore, PCPV is not causing a high communication cost while trying to reduce the operational cost sharply. As a result the performance of PCPV is increased where the impact of operational cost is more than communication cost and total cost is determined by almost the operational cost. This tendency is strongly affected by the cost parameters such as VNF license cost, site license cost and communication cost per unit of load. As to execution time, in a small size environment with 9 states and different numbers of Eus, PCPV takes around 0.02 second to run whereas ILP needs almost a day. Also, in the scenario II PCPV has less than a minute execution time even for a highest EUs number and states (29.1522 seconds for 625 states and 200 EUs). In addition, by increasing the number of VNFs per chain from 3, to 6 and 9 in scenario III, the execution time remains low (29.4656 and 28.9628 respectively).

V. CONCLUSION AND FUTURE WORK

This paper proposed an ILP model and a pattern-based proactive VNF placement algorithm to guarantee the service DT in offering VNF based VASs to EUs. The goal was to place VNFs in a way that it leads to the optimum number of VNFs to reduce the cost while still satisfying the service DT.

Fig. 1 . 3

 13 Fig. 1. 3 VNFs per chain -9 states (Scenario I)

Fig. 1

 1 Fig.1shows PCPV performance in the worst case where the communication cost is almost 10 times the operational cost. Fig.1(d)shows the total cost of PCPV and ILP for different numbers of EUs. The results prove that PCPV performance is close to ILP. However, ILP outperforms PCPV in communication cost by improving 30-46% (Fig.1(c)). Whereas PCPV utilizes fewer severs (Fig.1(a)) and improves operational cost by 66-77% in comparison to ILP (Fig.1(b)) to serve different numbers of EU requests. Therefore, PCPV could keep its performance close to ILP by utilizing fewer servers and having less operational cost. Utilizing fewer servers in PCPV is due to its placement policy. Using partitioning and patterning techniques as well as detecting optimal zone enable PCPV to optimize number of VNFs and find appropriate servers.

TABLE I .

 I INPUT PARAMETERS AND VARIABLES

			Network Inputs		Service Inputes
	N	Set of servers in the network, N ⊆ V	H	Set of services hf requested by user f ∈ U
	U	Set of EUs in the network, U ⊆ V	K	Set of VNF types that constitute all services hf ∈
					H
	E	The set of edges (i.e., logical communication	V
			links) in the network
				αk	Software license cost of a VNF instance of type
					k ∈ K
	Bn	Bandwidth cost of unit load (1 Gbps) per hop	Tk,n	Processing delay of VNF instance of type k ∈ K
			server n		on server n ∈ N for unit load (1 Gbps)
	β (u,v) f	f along edge (u, v) ∈ E Bandwidth cost incurred by sending load of EU	Rk	Resource requirement for VNF type k ∈ K
	Cn	Capacity of server n ∈ N in terms of resource	Pk	Processing capacity (Gbps) of VNF type k ∈ K
			units
	γ		Site license cost	Lf	Load of EU f ∈ U
	δn	Operational cost for unit resource (vCPU) for server n	Dh f	QoS (i.e., Service Delay),threshold of service hf ∈ H
				Variables
	xk,n, j	1, if instance j of VNF type k is assigned to server n ∈ N and 0, otherwise	y u,v f , p,q	1, if edge (u, v) hosts VNF edge (p, q) of SC for EU f and 0, otherwise
	λ	f k,n, j	is mapped to its instance j on server n and 0, 1, if VNF type k belonging to SC of EU f	zn	1, if server n is used and 0, otherwise
			otherwise
	directional edges E. Let us also consider K as a set of VNFs
	of different types, such as Video Transcoder, Video Mixer and
	Video Compressor. A service request h f ∈ H, generated by EU
	f ∈ U, is represented by a directed graph

f NF Vf ⊆ K is a set of VNFs that constitute service hf ∈ H, f ∈ U BW(u,v) The bandwidth capacity of edge (u, v) ∈ E E f NF Set of VNF edges of the VNF chain for service hf ∈ H , f ∈ U, D(u,v) Delay of unit load (1 Gbps) for edge (u, v) ∈ E I k Set of VNF instances of type k ∈ K σ(u, v) Hop count of the edge (u, v) ∈ E

IEEE International Conference on Network Softwarization (NetSoft 2018) -Technical Sessions

We showed PCPV is scalable and its performance is close to ILP. As future work, PCPV will be integrated into a scaling algorithm to handle fluctuations in EU workload over time.


Shohreh Ahvar
email: shohreh.ahvar@telecom-sudparis.eu

Mohammad M Mirzaei
email: mehdi.mirzaei@ut.ac.ir

Jeremie Leguay
email: jeremie.leguay@huawei.com

Ehsan Ahvar
email: ehsan.ahvar@inria.fr

Ahmed M Medhat

Noel Crespi
email: noel.crespi@mines-telecom.fr

Roch Glitho
email: glitho@ece.concordia.ca

SET: a Simple and Effective Technique

published or not. The documents may come

I. INTRODUCTION

Nowadays computer network services are considered as a key component in keeping the network running at all times.

To provide various network services, the telecommunication service provider s network includes a number of middleboxes. They can support various types of functions, and for this reason, middleboxes are important for network operators.

By increasing network requirements in both scale and variety, service providers have to add new middleboxes and upgrade already existing middleboxes continuously. One recent research shows that the number of different middleboxes in enterprise network is comparable to the number of physical routers [START_REF] Bari | On orchestrating virtual network functions[END_REF]. However, adding and updating middleboxes comes out with high Capital Expenditures (CAPEX) and Operational Expenditures (OPEX). In addition, the services deployed traditionally are not scalable.

To address these issues, Network Functions Virtualization (NFV) [START_REF] Han | Network function virtualization: Challenges and opportunities for innovations[END_REF] [START_REF] Mijumbi | Network function virtualization: State-of-the-art and research challenges[END_REF] has been proposed to transform middleboxes from specialized hardware appliances to software running on inexpensive, commodity hardware (e.g., x86 servers with 10Gb NICs) [START_REF] Etsi | Etsi gs nfv 001, network function virtualization (nfv) use cases[END_REF].

An NFV-based (i.e., software-based) middlebox is a good solution to reduce CAPEX, as network operators no longer need to buy specialized hardware. Moreover, software-based middleboxes can be deployed and managed dynamically without the necessity of having network administrators which reduces OPEX [START_REF] Mijumbi | Network function virtualization: State-of-the-art and research challenges[END_REF]. In addition to improving OPEX and CAPEX costs, NFV gives an opportunity to network operators to manage their network functions easily. The software-based middleboxes are referred to Virtual Network Functions (VNFs) in NFV terminology.

In order to use VNFs, they can be placed on computational nodes (e.g., servers, switches, data centers) that meet their resource demands. The computational nodes must provide NFV Infrastructure (NFVI) functions to support the execution environment. The placed VNFs can also be chained together to provide a required service.

Placement and chaining of the VNFs can affect on both quality of service (QoS) and cost. For this reason, it has recently attracted attention of both network providers and researchers.

Many methods have been proposed to place chains of VNFs with the goal of minimizing the cost for the service providers [START_REF] Bari | On orchestrating virtual network functions[END_REF]- [START_REF] Mijumbi | Network function virtualization: State-of-the-art and research challenges[END_REF] and [START_REF] Luizelli | Piecing together the nfv provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF]- [START_REF] Benson | Network traffic characteristics of data centers in the wild[END_REF]. However, a few of them considered the effects of ordering of VNFs on minimizing the cost by taking into account the compression/decompression rate, capacity and requirement of VNFs and servers. Compression of traffic happens with firewalls and compression/decompression is a feature of media gateways.

This paper makes the following contributions. First, we question the effects of the VNF ordering on placement algorithms in the achievement of the SET objective (i.e., cost efficiency). To this end, we put ourselves in the shoes of designers and describe a typical VNF placement process through a use case. We then propose a Simple and Effective Technique (SET) which can work in combination with VNF placement and changing algorithms. SET considers different possible ordering for the functions (e.g., a firewall) which may not have ordering constraints.-This flexibility optimizes resource usage as it may reduce the number of activated functions and decrease the network footprint."

Afterwards, by combining SET with CCVP (i.e., a recently proposed VNF placement algorithm), we show that the combination of the proposed method with CCVP can actually outperform the CCVP alone with changing the order of VNFs in the chain.

The rest of the paper is organized as follows. In Section II, we overview related work. A use case is presented in Section III to demonstrate the importance of VNF ordering. Section IV presents the main contribution of this paper which is basically the SET technique. The assessment of SET is covered in Section V. Finally, Section VI concludes the paper.

II. RELATED WORK

Several VNF placement and chaining algorithms which pay special attention to aspect of cost efficiency have been recently proposed. Some of them focused on specific individual costs. Luizelli et al. [START_REF] Luizelli | Piecing together the nfv provisioning puzzle: Efficient placement and chaining of virtual network functions[END_REF] and Fang et al. [START_REF] Fang | Joint spectrum and it resource allocation for efficient vnf service chaining in inter-datacenter elastic optical networks[END_REF] tried to minimize the number of instances, Moens et al. [START_REF] Moens | Vnf-p: A model for efficient placement of virtualized network functions[END_REF] minimized the number of used servers and Qu et al. [START_REF] Qu | Reliability-aware service provisioning in nfv-enabled enterprise datacenter networks[END_REF] targeted minimizing the communication bandwidth usage.

Unlike the above-mentioned works with the simple objective, Ghaznavi et al. [START_REF] Ghaznavi | Elastic virtual network function placement[END_REF] presented a solution where the elasticity overhead and the trade-off between bandwidth and host resource consumption were considered for placement of a same type of VNF. In another study, Ghaznavi et al. [START_REF] Ghaznavi | Service function chaining simplified[END_REF] completed their work for multiple VNF instances placement. Mechtri et al. [START_REF] Mechtri | A scalable algorithm for the placement of service function chains[END_REF] tried to maximize the providers revenue based on the number of accepted CPUs and bandwidth resources.

Riggio et al. [START_REF] Riggio | Virtual network functions orchestration in enterprise wlans[END_REF] minimized the links and nodes utilization to increase the accepted service chain requests in enterprise WLANs. The authors then extended their work in [START_REF] Riggio | Scheduling wireless virtual networks functions[END_REF] where a VNF placement heuristic, called WiNE (Wireless Network Embedding), was proposed. Finally, a few studies have attempted to consider more comprehensive cost models. Lin et al. [START_REF] Lin | Demand-aware network function placement[END_REF] presented an MILP and Game Theory based VNF placement with the goal of minimizing the cost of deploying VNF instances as well as the computing and network cost in optical networks. Zeng et al. [START_REF] Zeng | Orchestrating tree-type vnf forwarding graphs in inter-dc elastic optical networks[END_REF] considered the cost of IT resource and spectrum utilization of fiber links as their objective in addition to the cost of VNF deployment (instantiating) for the VNF placement in optical data centers. Bouet et al. [START_REF] Bouet | Cost-based placement of vdpi functions in nfv infrastructures[END_REF] proposed a solution where the cost includes the network as well as the license cost per site and vCPU in VNF instances. Bari et al. [START_REF] Bari | On orchestrating virtual network functions[END_REF] considered also a penalty cost to be paid to the customer for the Service Level Objective (SLO) violations.

Our previous work [START_REF] Ahvar | Ccvp: Cost-efficient centrality-based vnf placement and chaining algorithm for network service provisioning[END_REF] considered a complete cost model including license, computing and communication cost. They introduced a cost adaptive model for low and high traffics and different units of cost. None of the above mentioned works have considered the effects of VNFs ordering of the chain on cost where modification rate, VNFs capacity and requirement are taken into account in VNF placement and chaining algorithms. By introducing a use case, next section shows this effect.

There are also quite a few works which have considered ordering of VNFs. Allybokus et al. [START_REF] Allybokus | Virtual Function Placement for Service Chaining with Partial Orders and Anti-Affinity Rules[END_REF] considered partial ordering as a constraint and proposed a heuristic solution based on a linear relaxation of the problem. However, their work does not consider the rate modification for the VNFs which makes the effect of ordering more critical. Also, our work is different from their study, because we provide a simple and effective technique to improve the cost efficiency which works with any VNF placement and chaining algorithm. When we have rate modification feature in a VNF, then the next VNF in the chain receives less traffic, so it can be shared between more flows and fewer VNF instances can be used as a result. In this regard, Addis et al. [START_REF] Addis | Virtual network functions placement and routing optimization[END_REF] considered this feature for the VNFs which is called compress/decompress rate. In their work, they introduced tunneling VNFs which require decompression. However, unlike our work, they have not considered the effect of this feature on having different cost for different order of VNFs in the chain.

III. USE CASE

In this section, we present a use case to show the effects of VNF ordering on cost where the VNF placement algorithm takes into account the VNFs requirement and capacity in the placement and chaining of VNFs. To this aim, we define a chain with a partial order constraint for its VNFs. Let us consider a chain including 3 VNFs with the size of 1 vCPU and 2 GRAMs. Assume V1 is a Deep Packet Inspection (DPI), V2 is a traffic shaper (we call it QoS VNF), and V3 is a firewall (i.e., FW). We can also assume that all the flows pass through the same service with the same order. The objective is to find optimal order and placement for VNFs.

After considering the ordering constrains, 3 ordering options are possible : Option 1-DPI, QoS, FW; Option 2-DPI,FW,QoS; and Option 3-QoS.FW,DPI.

We implemented one of the recently proposed cost-efficient VNF placement and chaining algorithms (CCVP [START_REF] Ahvar | Ccvp: Cost-efficient centrality-based vnf placement and chaining algorithm for network service provisioning[END_REF]) which considers VNF capacity and requirement in taking decision for VNF placement and chaining. In CCVP, the same service with the same order is used for all the flows. In order to better evaluate the performance of our SET method, we extend CCVP by considering the rate modification for the VNFs as well.

CCVP first creates a graph G of the network topology. In the second step, based on the graph G characteristics and other received information (e.g., available instances and requested flows), it selects appropriate servers (i.e., server with highest centrality and enough capacity), deploys instances on them and indicates each flow passes through which selected server for all VNF types of the chain. The procedure of placing new instances of VNFs of the chain is continued until the total cost becomes larger than the former ones (and the former one could For the topology, we consider a real platform called Soft-FIRE. Three European testbeds federated in SoftFIRE project are considered in this use case. A computing cost is considered for each testbed based on Azure cost model. Since our VNF size is considered as size of 1 vCPU and 2 GRAMs, we use the computing per instance per month for instance type A1 v2 in Azure. The cost for this instance based on the location of the testbeds in Azure is as follows : Ads (Europe)= 30.50 ,Surrey (UK) = 37.20 , Fokus (Germany): 34.97. For the communication cost (i.e., the sum of the bandwidth used by the chains in the network), we used communication cost for outbound data for each Gb per month in Azure (i.e., 0.696). Based on Azure cost model, the traffic inside the testbeds is free [START_REF] Han | Network function virtualization: Challenges and opportunities for innovations[END_REF]. Also, for VNFs License cost with refer to FortiGate-VM00 (http://www.avfirewalls.com/), we assume 1250 dollar per VNF. Capacity of the VNFs is considered 4 Gb/s, and for the testbeds available capacity we consider 2 vCPU, 4G RAM for each testbed. 2 traffics have been assumed: ADS -→ Surrey (T1), ADS -→ Fokus (T2). Loads of the flows are as follows: T1 = 1000 Gb/month and T2=3000 Gb/month. The CCVP placement decisions, based on this ordering as an input, are shown in figure 1).

With calculating the cost shown in I, we can see the effects of ordering on the cost when VNF placement and chaining algorithm takes into account VNFs capacity and requirement. Option 1 leads to deploy more VNF instances. Therefore it has the highest amount of license and computing cost. Also, the communication cost of the placement by CCVP for different ordering on all 3 options is different. Therefore, by trying differing ordering, we can reduce the cost for the provider. However, it takes time to check all the possible orders especially where the topology is large and the chain is long. Therefore, there is a need of a heuristic to narrow the search space. This section presents our proposed SET technique. As depicted in Section III, ordering of VNFs plays a significant role on the total cost of VNF placement and chaining algorithms. In the proposed method, we applied the Tabu search to find the best order. In fact, the SET technique is a Tabu search which prunes various feasible order of VNFs in the chain by considering the partial order constraints. Then, it selects the best permutation based on its placement costs. SET works with any VNF placement and chaining algorithm which calculates the cost of each permutation.

A. Tabu search

A feasible permutation f p : [f 1 , f 2 , ..., f n ] is an order of VNFs, f i , which satisfies partial order constraints. With n VNFs, there is at most n! feasible orders of VNFs. If we have k out of n! feasible permutations, then, there is a set of all feasible permutation F P : {f p i |i ∈ [1, k]}. The task of the Tabu search is to find the best order of VNFs (i.e., an order with the minimum cost, among all the feasible orders)

We define three actions for our Tabu search which can generate a new permutation np from another ordered VNFs. However, it is possible that the generated permutation is unfeasible np / ∈ F P . In such cases, the permutation is not considered. Our actions can generate all possible orders and are listed below:

• Insertion. This action inserts f i in place of f j and shifts others.

• Swap. This action swap f i and f j .

• Reversion. this action revers a given sub sequence of an order. The Tabu search starts with a random permutation ro ∈ F O as the current order. It then analyzes the neighbors of this order, which can be generated by authorized actions and after that, selects the best neighbor to be considered in the next iteration. When a neighbor is selected, the related action is added to the Tabu list. In each iteration, authorized actions are the actions which are not included in Tabu list and at the end of each step, the most old action is removed from the Tabu list. This steps is repeated in all iterations as long as the cost is under a predefined threshold or a fixed number of iterations is done.

B. The underlying algorithm

Our proposed technique is capable to work with all VNF placement and chaining algorithms. It is worth mentioning that the level of success for SET depends on the items that the VNF placement algorithm takes into account. CCVP already considers the capacity and requirement of VNFs and server properties. However, the VNFs with compression/decompression features were not considered in the algorithm decision. In this paper, we extend CCVP [START_REF] Ahvar | Ccvp: Cost-efficient centrality-based vnf placement and chaining algorithm for network service provisioning[END_REF] by adding the modification rate consideration for VNFs to perfectly show the effect of ordering on minimizing the cost.

The pseudo code of SET integrated with the extended version of CCVP algorithm is illustrated in Algorithm 1.

V. PERFORMANCE EVALUATION

The goal of this section is to show the effect of SET on cost efficiency of VNF placement algorithms.

SET has been implemented in MATLAB (the source code is available here: https://github.com/mehmir/set). The simulation setup and the metrics, which are used for the evaluation, are described in the following.

A. Simulation Setup

Network Topology: We used a real network topology from university data center research network [START_REF] Benson | Network traffic characteristics of data centers in the wild[END_REF] composed of 23 nodes and 42 links.

Traffic data set: the number of flows changes from 4 to 8. The sources and destinations are generated randomly by using MATLAB. It is worth highlighting that, even though the flows configuration is set randomly, once set, it remains fixed across the runs of all tested algorithms ensure comparability of the results. Chains of 5 VNFs have been applied in the simulations. The capacity of the physical servers, physical links, VNF instances and server running costs, which we have used in this section, have been derived from research papers [START_REF] Bari | On orchestrating virtual network functions[END_REF] [START_REF] Bouet | Cost-based placement of vdpi functions in nfv infrastructures[END_REF]. Table II shows the server and VNF data assumption which have been used in the evaluation. The capacity, requirement and compression/decompression rate for the VNFs are selected differently for the each VNF to show the effect of these items in differentiating the cost of different orders. The first, second and third numbers after each VNF in the following list, show the capacity, requirement (i.e., Required CPU cores for the VNF) and compression rates for the VNF respectively. VNF1 5000 1 1; VNF2 4000 2 0.8; VNF3 3000 3 0.9; VNF4 2000 4 1; VNF5 1000 5 1. We also considered a partial order constraint for the VNFs where VNF2 should be before VNF3.

B. Cost Evaluation

Integrating SET with VNF placement algorithms which take into account the server and VNFs properties can help to reduce the cost. We verify this claim by integrating SET with CCVP.

We set the maximum number of iterations to 20 for Tabu search. As we can see in the table III, for 8 flows more orders meet the partial order constraints.

Fig. 1 .

 1 Fig. 1. feasible options: (a) Option 1 placement, (b) Option 2 placement, (c) Option 3 placement

Algorithm 1 : 7 best neighbor cost = ∞ 8 for 15 current order = best neighbor 16 Add

 1781516 Pseudo code of SET integrated with CCVP 1 T abu list = { } 2 choose current order randomly 3 best = current order 4 best cost = ∞ 5 while all iterations done or best cost ≤ threshold do 6 best neighbor = [ ] action i ∈ T abu list do 9 Generate the neighbor of current order called N by action i 10 if N ∈ F P then 11 cost = CCV P (N ) 12 if cost ≤ best neighbor cost then 13 best neighbor cost = cost 14 best neighbor = N actions{best neighbor} to T abu list 17 remove oldest action from T abu list 18 if best neighbor cost ≤ best cost then 19 best cost = best neighbor cost 20 best = best neighbor 21 return (best, best cost)

Fig. 2

 2 depicts the results of the overall cost of the data center topology for 4 and 8 flows. As Fig. 2(a) shows, SET selects ordering ID 10 as input of CCVP. If we consider, without SET, we might select the ordering ID 8 as the input of CCVP, SET can improve cost of CCVP around 50% by selecting the right ordering (i.e., ordering ID 10) as input. Similarly, as Fig. 2(b) shows, SET could improve cost of the CCVP algorithm up to 40% for 8 flows. Both Figs. 2(a) and (b) validate that our proposed SET technique could noticeably improve cost of CCVP and it can

Fig. 2 .

 2 Fig. 2. Simulation results: Overall Cost for different orderings (a) for 4 flows, (b) for 8 flows

TABLE III ORDERING

 III OF VNFS IN CCVP

This paper presented a Simple and Effective Technique (SET), which works based on Tabu search, to consider VNF ordering for improving cost of VNF placement algorithms. We believe SET can be integrated with any VNF placement and chaining algorithm (which takes into account capacity, requirement and compression rate of VNFs) to reduce the cost for the provider. As a proof of concept, we tested SET with integrating it with an extended version of one of the recent cost-efficient VNF placement and chaining algorithms (i.e., CCVP). The results proved our idea. As future work, we plan to integrate SET with more VNF placement and chaining algorithms.


Sanda Lefteriu
email: sanda.lefteriu@imt-lille-douai.fr

Stefano Grivet-Talocia
email: stefano.grivet@polito.it

Broadband modeling of passive components via augmented equivalent circuit models

This paper presents a gray-box modeling technique for obtaining stable and passive circuit models of electrical components. Starting from frequency response measurements and from an initial basic equivalent circuit model, a fully automated procedure performs series and parallel augmentations, the value of the added components being estimated by minimizing the deviation of the model response from measurements. An iterative application of this procedure yields accurate fits, allowing a compact representation of parasitic and second-order effects through suitably and automatically connected circuit elements.

I. INTRODUCTION

The continuous increase in the operating frequencies of integrated circuits requires accurate component models for inclusion in simulators during the design and verification phases. The fundamental (functional) behavior of basic components can be represented by a single element (an inductance L or a capacitance C or a resistance R), which, however, cannot capture parasitics and second-order effects and can only match the response of the hardware component over a limited frequency band. Instead, equivalent circuit models (ECMs) for spiral inductors [START_REF] Cao | Frequency-independent equivalentcircuit model for on-chip spiral inductors[END_REF], [START_REF] Wang | A wide-band 2 π equivalent-circuit model for spiral inductors on silicon[END_REF], ceramic capacitors [START_REF] Sullivan | Improved distributed model for capacitors in high-performance packages[END_REF] or thin-film resistors [START_REF] Li | A high-frequency scalable model for thin-film resistor on gan substrate[END_REF] must take into account parasitics, including skin and proximity effects, that are relevant at high frequencies.

There are two main approaches for ECM extraction. Starting from frequency responses obtained by direct measurement or electromagnetic simulation, the first approach fits a rational model to the measured transfer function [START_REF] Gustavsen | Rational approximation of frequency domain responses by vector fitting[END_REF], [START_REF] Lefteriu | A new approach to modeling multiport systems from frequency domain data[END_REF], including passivity enforcement as a post-processing step [START_REF] Grivet-Talocia | Passive Macromodeling: Theory and Applications[END_REF]. This approach is purely black-box: although a circuit synthesis is always possible, the obtained ECM does not reveal any physical insight. The other approach starts from an initial basic circuit model (which is preserved) and enriches it by adding components that aim at reproducing high-order effects and parasitics. The present paper focuses on this second approach.

Previous approaches for automated broadband model generation did not ensure stability and passivity by construction. In [START_REF] Kolstad | A new circuit augmentation method for modeling of interconnects and passive components[END_REF], the error between the measurements and the response of the nominal circuit is modeled by a small number of pole-zero elements, the final model being found by blackbox fitting. Stability is enforced by pole-flipping and passivity is enforced a-posteriori. The same nonlinear optimization problem that we propose in this paper (see also [START_REF] Lefteriu | Topological fitting: broadband modeling of passive components via augmented equivalent circuit models[END_REF]) is solved in [START_REF] Paul | Broadband modeling of high-frequency microwave devices[END_REF] by expressing the unknown as a rational function, with its numerator and denominator coefficients found via least squares, stability being enforced by constraining positive coefficients. Passivity is then enforced via perturbation of the imaginary eigenvalues of a Hamiltonian matrix pencil [START_REF] Grivet-Talocia | Passive Macromodeling: Theory and Applications[END_REF]. These two approaches still make use of black-box fitting, although applied to incrementally update an initial basic ECM.

Our paper makes a step towards a fully-automated circuit generation for components described by frequency domain measurements, by avoiding black-box fitting. Starting from a nominal circuit (an incomplete ECM of a given component or subsystem), we propose to generate all possible augmentations of this circuit by inserting passive elements in series or in parallel to the existing ones. We treat each added circuit element independently, and we solve a nonlinear optimization problem via a Sanathanan-Koerner iteration combined with an alternating least-squares to find the element values that provide the best match to the given measurements. This process guarantees stability and passivity of the final model by construction, in addition to providing an automated placement of circuit elements required for parasitic representation.

The paper is structured as follows. Section II presents a review of the modified nodal analysis (MNA) used to describe electrical circuits. In Sec. III, we propose to perform all possible series and parallel augmentations of the nominal circuit. Section IV discusses the proposed approach for topological fitting. Lastly, Sec. V presents a numerical example validating our approach and Sec. VI concludes the paper.

II. CIRCUIT MODELING VIA MNA -SHORT REVIEW

MNA [START_REF] Ho | The modified nodal approach to network analysis[END_REF], [START_REF] Vlach | Computer methods for circuit analysis and design[END_REF] is a powerful technique for formulating and solving circuit equations. The reduced incidence matrix A describes the topology of the circuit: its rows correspond to the nodes; each column, corresponding to a single element, is all-zero except for a single 1 and -1 pair at the two rows pointing to the nodes to which this element is connected (the row corresponding to a global reference node is removed). Kirchhoff's current and voltage laws describe the topological equations: Ai b = 0 and v b = A T e, where i b are branch currents, v b , the branch voltages, and e, the nodal voltages.

For the four types of elements: shunt dynamic (capacitors), shunt resistive (conductances), series dynamic (inductances) and series resistive (resistors), the columns of the incidence matrix are partitioned accordingly, as A C , A G , A L , A R , respectively. We also define corresponding diagonal matrices C b , G b , L b , R b storing the values of the elements (including, as a particular case, empty matrices if no elements of a given class are present). The constitutive equations for voltage controlled (shunt) elements are i G = G b v G and i C = sC b v C , while for current controlled (series) elements, they are v R = R b i R and v L = sL b i L , where v ν , i ν collect branch voltages and currents on branches of type ν = {C, G, L, R}.

The unknowns of the MNA system are the node voltages and the branch currents for the series elements: x = e i L i R T , yielding a descriptor-form representation

E dae ẋ(t) = A dae x(t) + Bu(t), y(t) = Cx(t),
with transfer matrix H(s) = C (sE dae -A dae ) -1 B, where B and C are input and output mappings that depend on the adopted input-output representation (Z, Y , or S), and

E dae =   A C C b A T C 0 0 0 L b 0 0 0 0   , A dae =-   A G G b A T G A L A R -A T L 0 0 -A T R 0 R b   .

III. CIRCUIT AUGMENTATION

Given a nominal circuit, augmentations are performed:

1) in parallel with shunt C and G elements between any two nodes, if no such elements exist. New columns are appended to the incidence submatrices A C , A G , and the (unknown) element value is added to matrices C b , G b . 2) in series with L and R elements between any two nodes, if no such elements exist. The branch is split into seriesconnected branches by adding new nodes. The row size of the incidence submatrices A L and A R grows by two, some entries are modified to account for the split, while a new column is appended to account for the new element. Consistently, L b and R b are also augmented.

To minimize user interaction, we generate all possible series and parallel augmentations of a nominal (starting) circuit according to these rules. Initial component values are set to 0.

IV. ELEMENT PERTURBATION

The next step is to identify appropriate values for the newly added elements. In practice, all elements (including the preexisting ones) should be re-optimized due to the modified topology. Our approach is to perturb each element value independently to minimize the error between the response of the model and the given measurements, looping through all elements via an alternating least squares iteration.

A. Single element perturbation

We distinguish between the two different element types

• shunt element: if connected between nodes k and l, we define an all vanishing selector vector ξ with the same size as the MNA matrices E dae , A dae , except entries ξ k = 1 and ξ l = -1.

• series element: we define ξ = 0 n e L 0 n R T for an inductance and ξ = 0 n 0 n L e R T for a resistance, with 0 ν being a null vector of size ν = {n, n R , n L } (number of nodes except reference, of resistors, of inductors, respectively), and where e L is the unit vector with a single 1 at the position of the corresponding inductor, and similarly for e R .

We define G := sE dae -A dae . A change in the element value induces a perturbation G = G 0 +∆αξξ T , where G 0 refers to the nominal circuit, and ∆α is the perturbation in the value of the element. For a dynamic element, ∆α = s∆C i or ∆α = s∆L i . For static elements, ∆α = ∆G i or ∆α = ∆R i . The Sherman-Morrison formula yields

G -1 = G -1 0 - G -1 0 ∆αξξ T G -1 0 1 + ∆αξ T G -1 0 ξ .
Hence, the transfer function is

H(s) = CG -1 B = H 0 (s) -∆α F(s) 1 + ∆αK(s)
.

where

F(s) = CG -1 0 ξξ T G -1 0 B, K(s) = ξ T G -1 0 ξ
, and H 0 (s) = CG -1 0 B. We now define the error between the measurements H i and the response of the perturbed system at frequencies ω i as

E= N i=1 H(jω i ) -H i 2 F = N i=1 ∆H i -F i ∆α 1 + ∆α • K i 2 F , with ∆H i = H 0 (jω i ) -H i , F(jω i ) = F i and K(jω i ) = K i .
Naturally, we wish to minimize this error with respect to ∆α. The optimization parameter ∆α appears in the numerator and denominator, making this problem nonlinear. With the Sanathanan-Koerner (SK) iteration [START_REF] Sanathanan | Transfer function synthesis as a ratio of two complex polynomials[END_REF], we relax the problem by using ∆α k-1 from the previous iteration in the denominator, yielding a linear problem for the current ∆α k min

∆α k N i=1 ∆H i + (∆H i K i -F i ) ∆α k 1 + ∆α k-1 • K i 2 F
which is solved via linear least-squares.

B. Several elements

For perturbing several elements, we employ alternating least squares, by optimizing each element independently using Sec. IV-A. This iterative process is run until convergence or until a maximum number of iterations is reached. The identification is effective if the number of raw samples is much larger than the number of elements that are required for an accurate fit. We illustrate the proposed approach on a 2-port spiral inductor test case, initially characterized by 200 S-parameter samples from 100 MHz up to 10 GHz obtained from a fullwave electromagnetic simulation. The initial nominal circuit is a single inductance L nom connected between the two ports, both referenced to a common ground node, whose value is also to be identified. We assume that this inductance is representative only at low frequency. Therefore, we refine the ECM as follows.

1) An initial ECM with two shunt capacitors (Fig. 1) is identified using data only up to 5 GHz; 2) A second nominal ECM is generated by halving the component values and cascading two identical sections (Fig. 2) to better reproduce proximity effects [START_REF] Cao | Frequency-independent equivalentcircuit model for on-chip spiral inductors[END_REF], [START_REF] Wang | A wide-band 2 π equivalent-circuit model for spiral inductors on silicon[END_REF]. These values are then optimized to match the full frequency behavior.

3) The automated augmentation from Sec III is run on the above ECM without additional user intervention.

The resulting circuit is shown in Fig. 3 (some elements have small values and can be removed without loss of accuracy, however, to illustrate the automated augmentation, we show them nevertheless). Responses for this ECM match the field solver data closely on the whole frequency band (Fig. 4). Similar performance was observed on other examples (not shown due to space limitations).

VI. CONCLUSION

This paper aims to make a step towards automated generation of equivalent circuit models for components. These models are valid in a broad frequency band and are stable and passive by construction. We observed that the scheme works quite well on basic components such as spiral inductors and capacitors, leading to simple topologies with few elements.

The algorithm may run into difficulties if the size of the ECM grows too large, or if the topology of the nominal circuit is not appropriate for the application at hand. In the future, we will consider combining the SK process with a Newton iteration and use regularization for alternate least-squares, to improve and speed up convergence.

Fig. 2 .

 2 Fig. 1. π model in low frequency

Fig. 3 .

 3 Fig.3. Circuit obtained after optimizing the element values for the circuit obtained by all possible series and parallel augmentations of the 2π configuration in Fig.2

Fig. 4 .

 4 Fig.[START_REF] Li | A high-frequency scalable model for thin-film resistor on gan substrate[END_REF]. Fitting the augmented circuit on the entire frequency band; measurements are shown together with the models from Fig.1and Fig.3

  This work was supported by the ACRI foundation through the Young Investigator Training Program and the TD COST Action TD1307 European Model Reduction Network (EU-MORNET) as a Short Term Scientific Mission. Moreover, this work was completed in the framework of CE2I project (Convertisseur d'Energie Intgr Intelligent). CE2I is co-financed by European Union with the financial support of European Regional Development Fund (ERDF), French State and the French Region of Hauts-de-France.


Marc Rebillat
email: marc.rebillat@ensam.eu-nazih.mechbal@ensam.eu

Nazih Mechbal

Damage localization in composite plates using canonical polyadic decomposition of Lamb wave difference signals tensor

Damage localization in composite plates using canonical polyadic decomposition of Lamb wave difference signals tensor

Marc Rébillat, Nazih Mechbal

Introduction

Monitoring in real-time and autonomously the health state of structures is of high interest in the industry, and more specifically in the aeronautic and civil engineering applications fields. Such a process is referred to as Structural Health Monitoring (SHM) [START_REF] Balageas | Structural health monitoring[END_REF][START_REF] Farrar | An introduction to structural health monitoring[END_REF]. To achieve this goal, structures become "smart" in the sense that they are equipped with sensors, actuators, and artificial intelligence that allow them to state autonomously regarding their own health. One can compare smart structures with the human body which, thanks to its various senses and nerves, is able to assess if it has been hurt, where it has been hurt, and to estimate how severe it is. Following this analogy, the SHM process is classically decomposed into four steps: damage detection, localization, classification, and quantification [START_REF] Worden | The fundamental axioms of structural health monitoring[END_REF].

The structure under study is here a composite complex aeronautic nacelle excited by means of Lamb waves. To deploy SHM to such structure, it is equipped with piezoelectric elements that can be used both as sensors and actuators. Each element is actuated one by one using a tone burst at high frequency, produces a Lamb wave that propagates throughout the structure, and the resulting Lamb wave is finally measured by the other piezoelectric elements acting as sensors. If a structure equipped with 𝑁 piezoelectric elements and for which acquisition is performed over K samples is considered, one naturally ends up with a tensor 𝑴 ∈ ℝ 𝑁×𝑁×𝐾 at the end of the SHM process. To monitor the possible apparition of damage, measurements are first performed in a reference state to get a reference tensor 𝑹. Then, during the life cycle of the structure measurements at unknown states are performed and provides the tensor 𝑼. The tensor 𝛅 that corresponds to the difference between 𝑹 and 𝑼 is the basis of the detection, localization, classification, and quantification steps of SHM.

The three-dimensional nature of the difference tensor 𝜹 ∈ ℝ 𝑁×𝑁×𝐾 allows for the use of specific analysis tools [START_REF] Kolda | Tensor decompositions and applications[END_REF][START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF][START_REF] Fanaee-T Et | Tensor-based anomaly detection: An interdisciplinary survey[END_REF]. Even if during the last decade, tensors analysis have been widely applied for signal processing purposes [START_REF] Kolda | Tensor decompositions and applications[END_REF][START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF][START_REF] Fanaee-T Et | Tensor-based anomaly detection: An interdisciplinary survey[END_REF], they have found relatively few applications in SHM and reported applications mainly focused on the detection step. For example: damage detection based on tensors in a civil engineering context [START_REF] Cheema | On Structural Health Monitoring Using Tensor Analysis and Support Vector Machine with Artificial Negative Data[END_REF][START_REF] Prada | Three-way analysis of structural health monitoring data[END_REF], tensor-based damage detection for non-destructive evaluation of composite structures using ultrasounds [START_REF] You | Tensor-based ultrasonic data analysis for defect detection in fiber reinforced polymer (FRP) composites[END_REF], or application of tensors for denoising purposes in composite plates monitored by ultrasonic waves [START_REF] Zheng | PARAFAC Decomposition for Ultrasonic Wave Sensing of Fiber Bragg Grating Sensors: Procedure and Evaluation[END_REF]. Thus, to the knowledge of the authors, the advantages of tensors for damage localization by means of Lamb waves in composite plates has never been investigated.

The focus is thus put in this paper on the localization step of the SHM process. Classical methods for damage localization by means of Lamb waves in composite plates are usually based on a path by path analysis of data [START_REF] Su | Identification of damage using Lamb waves: from fundamentals to applications[END_REF][START_REF] Fendzi | A general Bayesian framework for ellipse-based and hyperbola-based damage localization in anisotropic composite plates[END_REF][START_REF] Fendzi | A data-driven temperature compensation approach for Structural Health Monitoring using Lamb waves[END_REF] (i.e. on only one row of the tensor 𝜹 ∈ ℝ 𝑁×𝑁×𝐾 ). Typically, the time-of-flight for each path "actuator-sensor" is extracted, and then using knowledge of the wave propagation speed and by means of triangulation the damage localization is estimated [START_REF] Su | Identification of damage using Lamb waves: from fundamentals to applications[END_REF][START_REF] Fendzi | A general Bayesian framework for ellipse-based and hyperbola-based damage localization in anisotropic composite plates[END_REF][START_REF] Fendzi | A data-driven temperature compensation approach for Structural Health Monitoring using Lamb waves[END_REF]. These methods thus process each path independently and then integrate all the information together to form a localization map from which damage localization is inferred. As SHM data are three-way, highly redundant, and correlated, a vector-based one-way approach as depicted above cannot capture all these relationships and correlations together. The tensors thus appear as a very promising tool able to improve damage localization by means of Lamb waves in composite plates. The paper is organized as follows: a tensor-based damage localization algorithm is detailed in Section 2 and then validated experimentally in Section 3. A conclusion and a discussion are then drawn in Section 4.

Tensor based damage localization

A simple physical model of Lamb wave propagation

Wave propagation within structures can be as a first approximation thought as a very simple physical phenomenon: waves propagate with a velocity 𝑣 in all directions around their excitation point and are attenuated with an attenuation factor 𝛼(𝑑) ∈ ℝ + that is distance dependent. Then, in the absence of damage, a wave 𝑠(𝑡) is sent by the piezoelectric element 𝑛 and the elements {𝑚 ∈ 1: 𝑁, 𝑚 ≠ 𝑛} receives the signals:

𝑠 𝑛𝑚 𝑅 (𝑡) = 𝛼(𝑑 𝑛𝑚 )𝑠 (𝑡 - 𝑑 𝑛𝑚 𝑓 𝑠 𝑣 ) Eq. 1
where 𝑑 𝑛𝑚 denotes the distance between the elements 𝑛 and 𝑚, 𝑡 the sampled time and 𝑓 𝑠 the sampling frequency. It is important to notice here that isotropy is assumed whereas it may not be perfectly the case in practice. Let's now introduce a damage at position 𝐷 within this structure. When the wave emitted by the element 𝑛 hit the damage, it is reflected, and a new wave is reemitted within the structure. As a first approximation, one can assume that damage acts as a secondary wave source that reemits any incoming wave in all directions with a reflection coefficient 𝛽 ∈ ℝ + . The signal received by the element 𝑚 is then:

𝑠 𝑛𝑚 𝑈 (𝑡) = 𝑠 𝑛𝑚 𝑅 (𝑡) + 𝛽𝛼(𝑑 𝑛𝐷 )𝛼(𝑑 𝐷𝑚 )𝑠 (𝑡 - 𝑓 𝑠 (𝑑 𝑛𝐷 + 𝑑 𝐷𝑚 ) 𝑣 ) Eq. 2
If the focus is now put on the difference signal, one has:

𝛿 𝑛𝑚 (𝑡) = 𝑠 𝑛𝑚 𝑈 (𝑡) -𝑠 𝑛𝑚 𝑅 (𝑡) = 𝛽𝛼(𝑑 𝑛𝐷 )𝛼(𝑑 𝐷𝑚 )𝑠 (𝑡 - 𝑓 𝑠 (𝑑 𝑛𝐷 + 𝑑 𝐷𝑚 ) 𝑣 ) Eq. 3
It is then possible to take the Fourier transform of this signal and one ends up with the following transfer function:

𝐻 𝑛𝑚𝑘 = Δ 𝑛𝑚 (𝑘 -𝑘 0 ) 𝑆(𝑘 -𝑘 0 ) = 𝛽𝛼(𝑑 𝑛𝐷 )𝛼(𝑑 𝐷𝑚 )𝑒𝑥𝑝 (-𝑖2𝜋𝑓 𝑠 (𝑘 -𝑘 0 ) ( 𝑑 𝑛𝐷 + 𝑑 𝐷𝑚 𝐾𝑣 ))
Eq. 4

where 𝑆(𝑘 -𝑘 0 ) denotes the Fourier transform of the input signal 𝑠(𝑡), Δ 𝑛𝑚 (𝑘 -𝑘 0 ) the Fourier transform of the difference signal 𝛿 𝑛𝑚 (𝑡), and 𝐾 the total number of samples. 𝑘 0 and 𝑘 𝑀 stands for the frequency indexes over which the phase analysis starts and stops. 𝐻 𝑛𝑚𝑘 is thus the 𝑘-th coefficient of the Fourier transform for the difference signal on the path "actuator 𝑛sensor 𝑚" It is then possible to compute the phase of each element:

Φ 𝑛𝑚𝑘 = 𝜙[H 𝑛𝑚𝑘 ] = -2𝜋𝑓 𝑠 (𝑘 -𝑘 0 ) ( 𝑑 𝑛𝐷 + 𝑑 𝐷𝑚 𝐾𝑣 ) Eq. 5
The tensor 𝚽 ∈ ℝ 𝑁×𝑁×𝐾 containing the coefficients Φ 𝑛𝑚𝑘 can be interpreted as a threeway tensor [START_REF] Kolda | Tensor decompositions and applications[END_REF][START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF][START_REF] Fanaee-T Et | Tensor-based anomaly detection: An interdisciplinary survey[END_REF].

Canonical polyadic decomposition of the phase tensor

The idea is now to be able to make a structure popping out of the tensor 𝚽 defined in Eq. 5. From tensors literature, it is well known that tensors can be decomposed using the Canonical Polyadic Decomposition (CPD) up to a rank 𝑅 [START_REF] Kolda | Tensor decompositions and applications[END_REF][START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF][START_REF] Fanaee-T Et | Tensor-based anomaly detection: An interdisciplinary survey[END_REF]. Such a decomposition consists in finding a triplet (𝒂 ∈ ℝ 𝑁×𝑅 , 𝒃 ∈ ℝ 𝑁×𝑅 , 𝒄 ∈ ℝ 𝐾×𝑅 ) that allows for a more compact representation of a given tensor (see Fig. 1). According to the notations of Fig. 1, the CPD of the phase tensor 𝚽 defined in Eq. 5 can be expressed as:

Φ 𝑛𝑚𝑘 = -∑ 𝑎 𝑛𝑟 𝑏 𝑚𝑟 𝑅 𝑟=1 𝑐 𝑘𝑟 Φ 𝑛𝑚𝑘 Eq. 6
What is interesting here is that by analyzing Eq. 5 and by exploiting the fact that ∀ 𝑖 ∈ [1, 𝑁] 𝑑 𝑖𝐷 = 𝑑 𝐷𝑖 , the tensor 𝚽 can be exactly expressed as a tensor of rank 𝑅 = 2 by choosing:

𝒂 = [ 𝑑 1𝐷 1 … … 𝑑 𝑁𝐷 1 ] 𝒃 = [ 1 𝑑 1𝐷 … … 1 𝑑 𝑁𝐷 ] 𝒄 = 2𝜋𝑓 𝑠 𝐾𝑣 [ (𝑘 0 -𝑘 0 ) (𝑘 0 -𝑘 0 ) … … (𝑘 𝑀 -𝑘 0 ) (𝑘 𝑀 -𝑘 0 ) ] Eq. 7
By looking in more detail at this tensor decomposition, it is particularly striking to notice that 𝒂 and 𝒃 both theoretically provide direct access to {𝑑 𝑖𝐷 } 𝑖∈ [1,𝑁] that are the distances between each piezoelectric element and the damage position. On the knowledge of these distances, damage localization is thus theoretically possible. Furthermore, 𝒄 is parametrized by 𝑣 the wave velocity within the material and by the signal processing parameters 𝑓 𝑠 and 𝐾 and is thus relatively easy to estimate.

In summary, it is demonstrated that the CPD of rank 𝑅 = 2 of the phase of the difference signals between a healthy and damaged states potentially provides direct access to all the distances between the piezoelectric elements and the damage, which could allow for damage localization.

Managing unicity of CPD

Unfortunately, even if very efficient numerical tools are available to compute a CPD for a given tensor [14], such a decomposition is not unique. The issue is that here not only a decomposition is sought, but a decomposition that can be physically interpreted according to Eq. 7. It is however possible here to exploit the powerful property that the CPD is unique up to a scaling and a permutation of its terms. This guarantees that a decomposition of arbitrary choice can always be found. Therefore, a meaningful way to descale a numerically obtained CPD, or to obtain a unique one that makes sense physically is needed.

Mathematically, for the phase tensor 𝚽, once a first decomposition (𝒂, 𝒃, 𝒄) has been numerically obtained, what is needed is to find {𝜆 𝐴1 , 𝜆 𝐴2 , 𝜆 𝐵1 , 𝜆 𝐵2 , 𝜆 𝐶1 , 𝜆 𝐶2 } such that:

[ 𝑑 1𝐷 1 … … 𝑑 𝑁𝐷 1 ] = 𝒂 [ 𝜆 𝐴1 0 0 𝜆 𝐴2 ] = [ 𝛼 1 𝛽 … … 𝛼 𝑁 𝛽 ] [ 𝜆 𝐴1 0 0 𝜆 𝐴2 ] Eq. 8 [ 1 𝑑 1𝐷 … … 1 𝑑 𝑁𝐷 ] = 𝒃 [ 𝜆 𝐵1 0 0 𝜆 𝐵2 ] = [ 𝛾 𝛿 1 … … 𝛾 𝛿 𝑁 ] [ 𝜆 𝐵1 0 0 𝜆 𝐵2 ]
Eq. 9

2𝜋𝑓 𝑠 𝑣𝐾 [ 𝑘 0 -𝑘 0 𝑘 0 -𝑘 0 … … 𝑘 𝑀 -𝑘 0 𝑘 𝑀 -𝑘 0 ] = 𝒄 [ 𝜆 𝐶1 0 0 𝜆 𝐶2 ] = [ 𝜇 1 𝜖 1 … … 𝜇 𝐾 𝜖 𝐾 ] [ 𝜆 𝐶1 0 0 𝜆 𝐶2 ] Eq. 10
satisfying the followings constraints [START_REF] Sidiropoulos | Tensor decomposition for signal processing and machine learning[END_REF]:

{ 𝜆 𝐴1 𝜆 𝐵1 𝜆 𝐶1 = 1 𝜆 𝐴2 𝜆 𝐵2 𝜆 𝐶2 = 1 Eq. 11
By analyzing the above equations, it can be easily seen that these coefficients can be estimated as:

𝜆 𝐴2 = 1/𝛽 𝜆 𝐶1 = 𝑚𝑒𝑎𝑛 [ 2𝜋𝑓 𝑠 (𝑘 -𝑘 0 ) 𝐾𝑣𝜇 𝑘 ] 𝜆 𝐴1 = 1 𝜆 𝐵1 𝜆 𝐶1 𝜆 𝐵1 = 1/𝛾 𝜆 𝐶2 = 𝑚𝑒𝑎𝑛 [ 2𝜋𝑓 𝑠 (𝑘 -𝑘 0 ) 𝐾𝑣𝜖 𝑘 ] 𝜆 𝐵2 = 1 𝜆 𝐴2 𝜆 𝐶2 Eq. 12
It is important here to notice that to go back from an arbitrary numerical CPD to a CPD that is physically relevant, the knowledge of the velocity 𝑣 is needed. This descaling factor, i.e. the velocity of Lamb waves within the material is here indeed necessary to convert a time-domain information (extracted from phase here) to a distance information as done in any classical localization algorithm [START_REF] Fendzi | A general Bayesian framework for ellipse-based and hyperbola-based damage localization in anisotropic composite plates[END_REF][START_REF] Fendzi | A data-driven temperature compensation approach for Structural Health Monitoring using Lamb waves[END_REF][START_REF] Su | Identification of damage using Lamb waves: from fundamentals to applications[END_REF]. However, the Lamb wave velocity 𝑣 can be very easily estimated from experimental data in the reference state by computing the times of arrivals of the first wave packets and making use of the known distances between piezoelectric elements.

In summary, it is shown here that starting from a numerical CPD of the phase tensor and using knowledge on the velocity 𝑣 of Lamb waves within the material under study derived from input experimental data, it is possible to access to two estimates of the distances between all the piezoelectric elements and the damage {𝑑 𝑖𝐷 𝑎 } 𝑖∈ [1,𝑁] and {𝑑 𝑖𝐷 𝑏 } 𝑖∈ [1,𝑁] .

Even if theoretically ∀𝑖 ∈ [1, 𝑁] 𝑑 𝑖𝐷 𝑎 = 𝑑 𝑖𝐷 𝑏 , this may not be the case in practice due to several factors (experimental noise, numerical issues, …) and it has thus been chosen to introduce and to keep the two notations.

Damage localization imaging

The last step of the damage localization algorithm consists now in drawing a map able to highlight the most probable damage localization from the two sets of distances {𝑑 𝑖𝐷 𝑎 } 𝑖∈ [1,𝑁] and {𝑑 𝑖𝐷 𝑏 } 𝑖∈[1,𝑁] estimated previously. Let's consider a structure under study over which coordinates of a current point 𝑃 can be defined. It is then possible to compute for any current point on the structure the distances {𝑑 𝑖𝑃 } 𝑖∈ [1,𝑁] between this current point and the 𝑁 piezoelectric elements. The point of the structure that will most probably be the damage location should thus in theory satisfy:

∀𝑖 ∈ [1, 𝑁] 𝑑 𝑖𝑃 = 𝑑 𝑖𝐷 𝑎 = 𝑑 𝑖𝐷 𝑏 Eq. 13
Consequently, a very intuitive damage localization index (DLI) can be defined as:

DLI(𝑃) = 1 ∑ (2𝑑 𝑖𝑃 -𝑑 𝑖𝐷 𝑎 -𝑑 𝑖𝐷 𝑏 ) 2 𝑁 𝑖=1
Eq. 14 Finally, the damage imaging algorithm simply consists in plotting DLI(𝑃) over the structure and in searching for its maximum value. One can notice that the choice done here is arbitrary and that other potentially more optimal choices could have been made (minimizing the sum of squares, separate optimization procedures using 𝑑 𝑖𝐷 𝑎 and 𝑑 𝑖𝐷 𝑏 , …). The damage localization algorithm proposed here can thus be summarized as follows:

Algorithm overview

-Step 1: Compute the difference tensor 𝜹 between a reference and an unknown state. -Step 2: Compute the phase of the Fourier transform for the difference signal on the path "actuator 𝑛sensor 𝑚" and build the tensor 𝚽 on this basis (see Sec.

Experimental results

Experimental setup

The geometrically complex aeronautics structure under study consists here in the fan cowl part of a nacelle of an Airbus A380. This structure is 1.5 m in height for a semi circumference of 4 m and is made of composite monolithic carbon epoxy material. It has been equipped with 30 piezoelectric elements manufactured by NOLIAC (diameter of 25 mm) and possesses many stiffeners delimiting various areas; as shown in Figure 3. This part is particularly challenging because it contains several stiffeners. The damage is simulated using two 35 mm Neodymium magnets placed on both faces of the structure at the position indicated the yellow "X". The excitation signal sent to the PZT element is a "5 cycles burst" with a central frequency of 𝑓 0 = 200 kHz and with an amplitude of 10 V. This signal has been chosen here to maximize the propagation of the 𝑆 0 mode [START_REF] Su | Identification of damage using Lamb waves: from fundamentals to applications[END_REF][START_REF] Fendzi | A general Bayesian framework for ellipse-based and hyperbola-based damage localization in anisotropic composite plates[END_REF]. In each phase of the experimental procedure, one PZT is selected as the actuator and the other act as sensors. All the PZTs act sequentially as actuators. Resulting signals are then simultaneously recorded by the others piezoelectric element and consist of 1500 data points sampled at 1 MHz. The Lamb wave propagation speed within the material is estimated around 5200 m/s for the 𝑆 0 mode. Signals were acquired 10 times in both the healthy (reference) and damaged (unknown) states.

As pre-processing steps, the measured signals are first denoised by means of a discrete wavelet transform up to the order 4 using the "db40" wavelet. Those signals are then filtered around their center frequency 𝑓 0 using a continuous wavelet transformation based on "morlet" wavelets and with a scale resolution equals to 20. The diaphonic part present in the measured signals (i.e. the copy of the input signals that appears on the measured signal due to electromagnetic coupling in wires) has been previously eliminated based on the knowledge of the geometrical positions of the PZT and of the estimated waves propagation speed 𝑣 in the material.

CPD of the phase tensor in practice

To build the tensors 𝚽 and 𝑨 (see Sec. 2.1), the matrix 𝜹 containing the differences of signals between the healthy and damaged states have been built. To remove undesirable reflections from the differences signals, only the first wave packets have been retained in the difference signals. In Figure 4 [Left], the resulting normalized pre-processed difference signals are plotted in blue as well as the corresponding input signals (in red). From this figure it can be seen that the simplified underlying hypothesis leading to Eq. 3 is well satisfied after the pre-processing steps (denoising and first wave packet isolation): indeed, the difference signals contain a single reflection that arrive to sensor with variable delays. It should also be noted that in practice due to the data acquisition system being used, the piezoelectric elements can be considered either as actuators or as sensors, but not as both. Thus, nothing is measured on the "diagonal" part of the tensor (i.e. when 𝑛 = 𝑚, see the diagonal of Figure 4). In practice the matrices 𝜹 and 𝚽 are thus only partially known. The next step consists in computing first numerical CPDs of the phase and amplitude tensors. These CPDs are computed using the TensorLab toolbox running in a Matlab environment [14]. Obtaining a numerical CPD is nothing else than solving an optimization problem satisfying some constraints associated with the particular form of the decomposition being sought. Thus, initial values of the triplets (𝒂 ∈ ℝ 𝑁×𝑅 , 𝒃 ∈ ℝ 𝑁×𝑅 , 𝒄 ∈ ℝ 𝐾×𝑅 ) must be provided to the optimization algorithm before running it. Here 𝒄 is initialized according to Eq. 7 as information regarding Lamb waves velocity 𝑣 has been previously estimated. The matrices 𝒂 and 𝒃 are initialized by considering that an initial guess localization for the damage is the barycenter of all the positions of the piezoelectric elements. The CPDs are then obtained through a nonlinear least squares algorithm. Once the two numerical CPDs are obtained, they are descaled to allow for their physical interpretation as explained in Sec. 2.3. At that moment it is then possible to access to two estimates of the distances between all the piezoelectric elements and the damage {𝑑 𝑖𝐷 𝑎 } 𝑖∈ [1,𝑁] , and {𝑑 𝑖𝐷 𝑏 } 𝑖∈[1,𝑁] (see Sec. 2.3) and to compute the damage localization index defined by Eq. 14 for all points 𝑀 in the area of interest where the damage could be located. As for each case under study 10 repetitions in the healthy and damaged states are available, the damage localization index is computed 100 times. Localization results are provided in Figure 5. On this figure, the damage localization index is plotted for the case under study as well as the positions of the different piezoelectric elements. In addition to the resulting damage localization maps, the position of the actual damage as well as the position of the best estimation and of the median estimation provided by the algorithm are provided. Furthermore, the associated "Best errors" and "Median errors' are reported in the title and are both around 27 cm. The median estimated localization for each case under study lies relatively close to the actual damage position when considering the overall structure dimension. It can also be seen that the damage localization maps provided by the damage localization index graphically indicate the correct area and that the best localization results are really close to the actual damage position. Thus, the relative localization error (i.e. the localization error divided by the larger portion of the structure expressed in %) is ≃ 6 -7% for the case being considered. Another noticeable point with respect to the proposed algorithm is that it takes as inputs directly the signals measured by all the piezoelectric elements without any additional preprocessing or selection steps. This makes the algorithm relatively robust in the present case.

Localization results

CONCLUSION

Monitoring in real-time and autonomously the health state of aeronautic structures is referred to as Structural Health Monitoring (SHM) and is a process decomposed in four steps: damage detection, localization, classification, and quantification. Structure under study was here a complex aeronautic nacelle and the focus was put on the localization step of the SHM process. The fact that SHM data are naturally three-way tensors has been investigated for this purpose. It is demonstrated in this paper that under classical assumptions regarding wave propagation, the canonical polyadic decomposition of rank

Fig. 1 :

 1 Fig. 1: Schematic representation of the CPD of the tensor 𝚽

Fig. 2 :

 2 Fig. 2: Overview of the damage localization algorithm

  2.1). -Step 3: Compute the CPD of rank 𝑅 = 2 of the tensor 𝚽 (see Sec.

Figure 3 :

 3 Figure 3: Fan cowl of an A380 Nacelle made of composite materials and equipped with 30 piezoelectric elements and possessing many stiffeners. Damage position is indicated with the yellow "X".

Figure 4 :

 4 Figure 4: [Left] Example of input signals (red) and of pre-processed difference signals (blue) for several "actuator-sensor" paths. Amplitudes have been normalized. [Right] Part of the phase tensor 𝚽 for several "actuator-sensor" paths. Phase is plotted in the range [0.9𝑓 0 , 1.1𝑓 0 ] with 𝑓 0 = 200 kHz.The phase and amplitude are then computed from the discrete Fourier transform of these signals. As input signals are band limited around their central frequency 𝑓 0 , only the phase and amplitude in the range [0.9𝑓 0 , 1.1𝑓 0 ] is considered here. The phases, which are parts of the phase tensor 𝚽 are plotted in Figure4[Right]. From this figure as expected from Eq. 5, the phases decrease linearly with the frequency and the slopes of these phases are different.The next step consists in computing first numerical CPDs of the phase and amplitude tensors. These CPDs are computed using the TensorLab toolbox running in a Matlab environment[14]. Obtaining a numerical CPD is nothing else than solving an optimization problem satisfying some constraints associated with the particular form of the decomposition being sought. Thus, initial values of the triplets (𝒂 ∈ ℝ 𝑁×𝑅 , 𝒃 ∈ ℝ 𝑁×𝑅 , 𝒄 ∈ ℝ 𝐾×𝑅 ) must be provided to the optimization algorithm before running it. Here 𝒄 is initialized according to Eq. 7 as information regarding Lamb waves velocity 𝑣 has

Figure 5 :

 5 Figure 5: Obtained localization results for the different damage localization indexes and for the full structure under study.

  2.2). -Step 4: Estimate Lamb wave velocity 𝑣 to descale the previous numerically obtained CPD and extract damage to piezoelectric elements distances estimates {𝑑 𝑖𝐷 𝑎 } 𝑖∈[1,𝑁] and {𝑑 𝑖𝐷 𝑏 } 𝑖∈[1,𝑁] (see Sec. 2.3). Compute the damage localization index and draw a damage localization map in order to estimate the most probable damage localization (see Sec. 2.4).

-Step 5:

damaged states provides direct access to the distances between the piezoelectric elements and the damage. This property is used here to propose an original and robust tensor-based damage localization algorithm. This algorithm is successfully validated on experimental data coming from the structure under study. Thanks to the tensor formalism this algorithm appears to be very robust when integrating all the information together in the localization process through the canonical polyadic decomposition. Future work will now focus on a methodological comparison between this tensor-based algorithm and classical algorithms thanks to numerical simulations and experimental data.