Marco António

Do Rosário

Helder Santos Cruz

Roberto De

Oliveira Rocha

Marcia Helena

Moreira Paiva

Marcelo Eduardo Vieira Segatto

Eglantine Camby
email: ecamby@ulb.ac.be

Gilles Caporossi
email: gilles.caporossi@hec.ca

Roberto De Oliveira Rocha

Marco A R S Cruz

Helder R O Rocha
email: helder.rocha@ufes.br

Marcia H M Paiva
email: marcia.paiva@ufes.br

Marcelo E V Segatto
email: segatto@ele.ufes.br

An algorithm for cost optimization of PMU and communication infrastructure in WAMS

Keywords: Phasor Measurement Unit (PMU), optimal PMU placement, Variable Neighbourhood Search (VNS), communication infrastructure, power system observability

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The impacts of increasing renewable power, free access to transmission networks and deregulation of the electricity sector are changing the nature of power systems. This fact brings more stress to electrical power system operation, leading to an increasing need to improve the control and protection systems offered by Wide Area Measurement Systems (WAMS) based on Phasor Measurement Units (PMU) [START_REF] Bobba | Enhancing grid measurements: Wide area measurement systems, NASPInet, and security[END_REF].

State Estimation (SE) processes a redundant set of measurements in order to obtain in real time a reliable estimate of the operating state of the power system. Conventionally, network state monitoring is done with conventional meters (magnitude of bus, active/reactive power flows and injections along with some voltage magnitude) from which measurements are transmitted to the Supervisory Control and Data Acquisition Systems (SCADA). With the introduction of PMUs in the 1990's, the conventional meters with Remote Terminal Units (RTUs) are now being complemented by PMUs which transmit magnitude and phase measurements to the SCADA system. The PMUs are synchronized through a Global Positioning System (GPS), which allows WAMS to be in control and operation of power systems. Another advantage is the PMU high data sampling combined with SE process that ensures quick voltage control compared to conventional measurements. However, PMUs and all related synchrophasor infrastructure deployment can be very costly [START_REF]Factors affecting PMU installation costs[END_REF] in large power systems, thus provoking the Optimal PMU Placement (OPP) problem on transmission networks.

Many algorithms have been proposed to solve the OPP problem by taking into account criteria such as observability, contingencies (single PMU outage or single branch failure events), critical measurements and critical assemblies.

Most of these algorithms focus on minimizing PMU costs leading to objective functions which simply optimize the number of PMUs and their placements. A few recent studies have called attention to other aspects such as Communication Infrastructure (CI), which makes a greater contribution to monitoring system costs than PMUs. In such studies, the problem is formulated in a more comprehensive way, seeking to optimize the allocation of PMUs and communication infrastructure for a minimum total cost. Centralized WAMS architecture (Figure 1) is very suitable due to an efficient use of control elements, lower cost, good coordination of alarms and event management and lower latency. This performance is a consequence of a smaller number of Phasor Data Concentrators (PDCs) although it is exposed to a greater probability of failure than a distributed architecture [START_REF] Gore | Analysis of Wide Area Monitoring System architectures[END_REF].

Visualization

Central PDC Analytics Among the transmission media for communication in smart grids, dependent media, i.e., those which are part of a power system and are owned by independent system operators (ISO) such as Power Line Communication (PLC), All-Dielectric Self Supporting (ADSS) and Optical Power Ground Wire (OPGW), are preferred due to lower latency compared with other medias [START_REF] Shahraeini | Co-optimal placement of measurement devices and their related communication infrastructure in wide area measurement systems[END_REF]. However, OPGW stands out due to the advantages of high channel capacity and transfer rate, low transmission losses and immunity to electromagnetic interference [START_REF] Mohammadi | A new approach for optimal placement of PMUs and their required communication infrastructure in order to minimize the cost of the WAMS[END_REF][START_REF] Kansal | Bandwidth and latency requirements for smart transmission grid applications[END_REF]. Furthermore, the CI can be co-optimally designed in conjunction with power system planning problems.

CONTROL CENTER

PMU

Contribution

This paper proposes new methods and algorithms for the OPP problem based on a topological observability analysis and for simultaneously optimizing PMU allocation and communication network leading to a minimum total cost in WAMS. Other costs besides PMU and CI are outside the scope of this work. The algorithm for the latter problem is derived from the first and include Dijkstra's single-source shortest-path and Kruskal minimum spanning tree algorithms to achieve a minimum total cost. To construct both algorithms, a

Variable Neighbourhood Search (VNS) heuristic [START_REF] Hansen | Variable neighborhood search[END_REF] is implemented. To the best of our knowledge, it is the first time that a VNS heuristic is applied to the OPP problem.

The problems above are related to graph theory concepts such as Dominating Set and Connected Dominating Set [START_REF] Liao | Hybrid search for the optimal PMU placement problem on a power grid[END_REF][START_REF] Meghanathan | Use of Centrality Metrics to Determine Connected Dominating Sets for Real-World Network Graphs[END_REF]. These concepts have motivated the construction of a new metric, named here Dominance, which plays a fundamental role in improving the efficiency of the proposed algorithms in combination with other known graph theory metrics. Our algorithms have proved to be flexible as they incorporate various cost parameters and contingencies. Additionally, they are scalable when generating results from small networks to large and complex electrical transmission networks and show good efficiency when the results are compared with other works.

This paper proposes an optimization of costs on a centralized WAMS supported by OPGW transmission media, in considering the costs of the switches.

These are understood as electro-optical and/or optic-optical converters integrated into an Add Drop Multiplexer (ADM), which enable PMU connection to the OPGW communication network or allow the various OPGW branches arriving on a bus to form a communication network node.

Related work

There are several works on Optimal PMU Placement (OPP) which have been developed over approximately 25 years [START_REF] Noureen | Phasor measurement unit integration: A review on optimal PMU placement methods in power system[END_REF]. These algorithms are classified according to the type of observability analysis, numerical or topological [START_REF] Sodhi | Optimal PMU placement method for complete topological and numerical observability of power system[END_REF].

Despite the fact that topological methods are faster and more suitable for large systems, some numerical algorithms have been shown to be more efficient in allocating PMUs [START_REF] Korres | Optimal phasor measurement unit placement for numerical observability in the presence of conventional measurements using semi-definite programming[END_REF].

All these algorithms can still choose to use deterministic [START_REF] Huang | Optimal PMU placement considering controlled islanding of power system[END_REF] or metaheuristic methods [START_REF] Nazari-Heris | Application of heuristic algorithms to optimal PMU placement in electric power systems: An updated review[END_REF]. Although metaheuristic algorithms have advantages in computational time (especially in large systems) and in the ability to deal with conditions such as single PMU outage, topological variation, critical measurements and areas of interest, trapping in local minima is still a possible problem [START_REF] Monti | Phasor measurement units and wide area monitoring systems[END_REF].

The consideration of realistic costs and practical constraints related to substations have either led to different approaches beyond the OPP problem -such as the Optimal Substation Coverage Algorithm [START_REF] Pal | General optimal substation coverage algorithm for phasor measurement unit placement in practical systems[END_REF] -or have forced the OPP problem to consider the cost of upgrading substations for PMU placement [START_REF] Pal | A PMU placement scheme considering realistic costs and modern trends in relaying[END_REF].

Communication Infrastructure (CI) is also another factor which is influencing the OPP problem beyond conventional PMU cost-minimization paradigm. In this sense, the problem is modeled as a total cost minimization problem [START_REF] Shahraeini | Co-optimal placement of measurement devices and their related communication infrastructure in wide area measurement systems[END_REF],

taking into account both PMU and CI costs and using a genetic algorithm to confirm that the lowest total cost does not always correspond to the smallest number of PMUs. The paper [START_REF] Dubey | Cost Effective Wide Area Measurement Systems for Smart Power Network[END_REF] has a similar approach and proposes a multiobjective algorithm to address PMU, PDC and CI placement, simultaneously.

In [START_REF] Mohammadi | A new approach for optimal placement of PMUs and their required communication infrastructure in order to minimize the cost of the WAMS[END_REF], this previous approach is explored in greater depth and other aspects and devices involved in communication infrastructure are detailed. Additionally, that work considers contingency conditions and pre-existence of some PMUs and communication cables in certain parts of the network. In the same way, [START_REF] Singh | Optimal cost wide area measurement system incorporating communication infrastructure[END_REF] considered CI cost and several contingency conditions in a binary gravitational search algorithm, which showed better performance than in the previous work.

Still, the same approach to total cost optimization is highlighted in [START_REF] Rather | Realistic approach for phasor measurement unit placement: Consideration of practical hidden costs[END_REF] and a realistic cost-effective model is presented for optimal PMU placement which considers practical and unaccounted cost implications based on a reallife project. The results confirmed that a minimal number of PMUs does not necessarily indicate minimum financial implications of the OPP project.

As an alternative to the existing fiber-based communication infrastructure, microwave communication technology was considered in an OPP problem [START_REF] Appasani | Co-Optimal Placement of PMUs and Their Communication Infrastructure for Minimization of Propagation Delay in the WAMS[END_REF] with focus on minimizing the propagation delay.

The remaining sections of this paper are organized as follows. Section 2

gives an introduction to the fundamental concepts used throughout the paper.

Section 3 presents the OPP problem and a VNS-based heuristic to solve it.

Section 4 extends the problem to a total cost objective-function and details the algorithm implementation. Section 5 provides simulation results and finally, the conclusions are drawn in Section 6.

Technical Background

This section covers the theoretical fundamentals used throughout this paper: the basis of Metaheuristic Variable Neighbourhood Search, graph theory concepts and algorithms, and topological observability analysis.

Metaheuristic Variable Neighbourhood Search

The Variable Neighbourhood Search (VNS) is a metaheuristic introduced by [START_REF] Hansen | Variable neighborhood search[END_REF] that counts on several successful applications in solving combinatorial and global optimization problems [START_REF] Caporossi | Variable neighborhood search for extremal 31 graphs. Three ways to automate finding conjectures[END_REF]. VNS relies on systematic changes in the vicinity for local search. This metaheuristic is characterized by simplicity and efficiency.

VNS is based on a variable neighbourhood method which exploits neighbourhoods progressively further away from the current solution and moves to another solution if it is better than the previous one. This metaheuristic has three stages to search for new solutions: improve the solution prior to a perturbation, then the new solution is adopted and the perturbation is resumed with value k = 1. This process is iterative and continues until a stopping criterion is reached.

The basic structure of the VNS is presented below although there are some variations of the algorithm [START_REF] Hansen | Variable neighborhood search[END_REF]:

Basic VNS Algorithm Function VNS (x, k max , t max ) 1: Repeat 2: k ← 1 3: Repeat 4:
x ← Shake(x, k)

5:

x ← FirstImprovement(x ) 6:

x ← NeighbourhoodChange(x, x , k) To find shortest paths between vertices in a graph G = (V, E) with nonnegative edge weights, the single-source shortest path Dijkstra algorithm is used.

7: Until k = k max ; 8: t ← CpuTime() 9: Until t > t max ;

Graph Theory Concepts and Algorithms

A graph G = (V, E) is a set of |V |
The running time of this algorithm is O(|V |log|V |) [START_REF] Cormen | Introduction to algorithms[END_REF].

The spanning tree of a graph G = (V, E) consists of a connected acyclic subgraph T = (V, E ) which contains all the vertices of G and a subset E ⊂ E of its edges. To find a minimum-spanning-tree, that is, a spanning tree minimizing the sum of edge weights, the Kruskal Algorithm is used. It can easily be made to run in time O(|E|log|V |) [START_REF] Cormen | Introduction to algorithms[END_REF].

The problem proposed in this paper is closely related to the minimum connected dominating set problem [START_REF] Brueni | The PMU placement problem[END_REF]. A dominant set (DS) of a graph

G = (V, E) is a subset D ⊂ V such that each element v of V is either in DS or is connected
to a vertex in DS. In the latter case, v is said to be covered by a vertex in DS.

A connected dominant set (CDS) is a connected subgraph induced by vertices of D. The problem related to these concepts is the cardinality minimization of DS or CDS, with some applications in wireless network problems [START_REF] Ding | Distributed Construction of Connected Dominating Sets with Minimum Routing Cost in Wireless Networks[END_REF].

Topological observability analysis

In a topological observability analysis, a power system is represented by a topological graph, where the vertices represent the network buses and the edges represent the network branches connecting the buses. Topological observability analysis is defined as the existence of at least one spanning measurement tree of full rank in a network [START_REF] Abiri | Optimal PMU placement method for complete topological observability of power system under various contingencies[END_REF]. This tree connects all observable buses and branches which can be observed by direct measurements or calculations. The existence of this tree is commonly assessed by the following rules [START_REF] Hajian | Optimal placement of PMUs to maintain network observability using a modified BPSO algorithm[END_REF]:

1. Assign with direct voltage phasor measurement and direct current phasor measurement of incident lines the buses with PMUs;

2. If voltage and current phasors at one end of a branch are known, then the unknown voltage phasor at the other end of the branch can be calculated (called a pseudo-measure) through Ohm's Law;

3. A current phasor of a branch can be calculated (pseudo-measure) if voltage phasors of both ends of this branch are known;

4. If all line current phasors incident to a zero-injection bus (ZIB) are known except one, the current phasor of the unknown one can be calculated through Kirchhoff's Circuit Law equations. ZIB is a bus in which the net power injection is zero. This is used as pseudo information to make a system observable with a lesser number of PMUs compared to the case when information of ZIBs is not considered;

5. If there is an unknown ZIB bus and voltage phasors of its adjacent buses are all known, then it is possible to calculate the voltage phasor of the ZIB bus.

Proposed Method for Optimal PMU Placement Problem

This section presents the OPP problem formulation considering only the PMU cost, and the proposed method to solve it. The formulation takes into account the following constraints: observability in normal condition and in N-1 contingency conditions, such as single PMU or branch outage.

A new metric designed specifically for the OPP problem, called Dominance Centrality, is introduced here, and a VNS-based algorithm is proposed to solve this problem.

PMU Placement Problem Formulation

For an N bus system, the optimal PMU placement problem is formulated as follows:

Min N i=1 cp i • x i (1)
s.t. A • X ≥ b, (2)
where:

cp i is the cost of placing a PMU at bus i;

A is the adjacency matrix;

b = [1 1 1 . . . 1] T is a unit vector of length N ;
X is the binary decision variable vector for PMU placement, which entries

x i =1 if a PMU is placed at bus i, and zero if otherwise.

Observability requirement is attained by Constraint (2) under normal conditions. However, this constraint has to be replaced by Equations ( 3) or (4) if observability is imposed under single PMU or single branch outage conditions, respectively.

A • X ≥ [2 2 2 . . . 2] T (3)
A j • X ≥ [1 1 1 . . . 1] T , j = 1, 2, . . . , N top (4)
In the latter equation, A j is a new adjacency matrix, slightly different from the original matrix A reflecting topology change inflicted by the branch outage.

N top is the number of possible scenarios of single branch outage.

The optimal PMU placement problem can also be formulated in terms of minimizing the dominating set (DS) of the graph representing the power network [START_REF] Brueni | The PMU placement problem[END_REF]. The graph theory terms relate to system observability in the following way:

• The vertices of a DS define the buses where PMUs are allocated;

• A bus that is covered (observable) is a bus with an allocated PMU or adjacent to a bus with a PMU; and

• The vertices of a DS compose a feasible solution for PMU placement thereby guaranteeing system observability.

Therefore, the problem of finding the buses to allocate PMUs on a power network ensuring system observability corresponds to the problem of finding a minimum DS for the graph that represents the power network.

Dominance Centrality, a new metric for DS

Degree Centrality is largely used in topological algorithms as the main searching criterion for OPP problems [START_REF] Meghanathan | Use of Centrality Metrics to Determine Connected Dominating Sets for Real-World Network Graphs[END_REF][START_REF] Roy | An optimal PMU placement technique for power system observability[END_REF] . However, this strategy is not very efficient in achieving a minimal number of PMUs due to the fact that some vertices with high degrees could cover (or share) the same set of adjacent vertices. To avoid or diminish multiple covering probability, the searching process should consider vertices with high degree but connected to neighbourhood vertices with low degrees. These informations are embedded together in a new metric called here Dominance Centrality (dom). The Dominance Centrality of a vertex v i , dom(v i ), is defined as follows:

dom(v i ) = deg(v i ) -adn(v i ), (5)
where deg(v i ) is the degree of vertex v i , and adn(v i ) is the average degree of all neighbourhood vertices of v i , which is formally defined in Eq. ( 6) as follows:

adn(v i ) = deg(vi) j=1 deg(v j ) • a ij deg(v i ) , (6)
where each v j is a vertex adjacent to v i .

Proposed Optimal PMU Placement Algorithm

This section describes the VNS-based heuristic proposed to solve the OPP problem defined in the previous section. To allow comparisons with other works, PMU unit cost cp i is taken as a constant value in the objective-function (1) stated previously, with the same constrained conditions as stated before.

The VNS is divided into three stages: (1) Initial solution, (2) Shaking and

(3) Local search (descent). These stages can be adapted to each formulation of the proposed problems.

Normal operating condition

Initial Solution: The adjacency matrix of the network, and the degree and Dominance of all its buses (vertices) are calculated and 33% of the buses are selected with the highest Dominance value for PMU allocation (dominant vertex). This percentage value is based on the work [START_REF] Baldwin | Power system observability with minimal phasor measurement placement[END_REF], in which the authors state that PMUs should be allocated to about one-quarter to one-third of the system buses to meet observability criteria. Although this is not true for many graphs, this range is used here as a reference to build an initial solution. Additionally, For branch failure conditions, the implementation looks very similar but the condition to remove PMUs in the local search stage is more relaxed since it depends on the branch failure scenario, i.e., a bus should be observed by two PMUs if a branch failure could affect connectivity to a bus with an allocated PMU, otherwise one PMU is sufficient (4).

PMUs

Proposed Method for PMU and Communication Infrastructure Cost Optimization Problem

This section presents the problem formulation and the method proposed to optimize both PMU placement and related communication infrastructure (CI) for a minimum total cost associated. The proposed method is derived from the previous OPP algorithm discussed in Section 3.

Problem Formulation

The PMU and Communication Infrastructure Cost Optimization Problem is stated in Equations ( 7)-( 10):

Min   N i=1 cp i • x i + N i=1 N j=1 (cf • d ij + cr • x i ) • y ij   (7)
s.t. A • X ≥ b (8) ij∈E y ij = N -1 (9)
ij∈E: i,j∈S

y ij ≤ |S| -1, ∀ S ⊆ V (10)
The objective function [START_REF] Hansen | Variable neighborhood search[END_REF] aims to minimize the total cost related to PMUs (first term) and communication infrastructure (second term). The CI cost has two components. One is related to active elements, such as switches and routers, and the other one is related to OPGW deployment, which is proportional to the length of cables deployed.

Constraint (8) ensures power system observability under normal conditions if b is a unit vector of length N . For N-1 contingency conditions, i.e., single PMU or branch outage, this constraint is replaced by respective Equations ( 3) or (4).

Constraints ( 9) and ( 10) ensure a tree topology for PMUs communication network. In this sense, Equation ( 9) imposes a number of links or edges (y ij ) that is one unit less than the number of vertices N , and Equation ( 10) ensures there is no cycle by imposing that any subset of S vertices must have at most (|S| -1) edges connecting them.

In this optimization, a centralized WAMS architecture is considered and a communication network with OPGW is optimized over the electric network.

In addition, the variables involved are defined as follows:

N is the number of buses in the network; cp i is the cost of placing a PMU at bus i; cr is the cost of active devices, such as routers and switches, allocated in all buses with PMUs;

cf is the OPGW cost per km;

d ij represents the distance between buses i and j; y ij ∈ {0, 1} assumes value 1 if an existing line (i, j) connecting buses i and j is selected to integrate the OPGW communication network;

X is the binary decision variable vector for PMU placement whose elements

x i are defined as one if a PMU is placed at bus i, and zero if otherwise;

A is the adjacency matrix of the network.

Proposed Method for Total Minimum Cost

The proposed problem involves allocating a minimum set of PMUs and also minimizing the distances involved in the fiber cable communication network.

This two-component objective is related to the Minimum Connected Dominating Set (MCDS) problem addressed in Graph Theory, which is in general NP-hard.

Some methods solve the MCDS problem directly whereas other methods propose a two-phase strategy by first building the DS and then connecting it to get a CDS [START_REF] Liu | The first constant factor approximation for minimum partial connected dominating set problem in growth-bounded graphs[END_REF]. This two-phase strategy is not appropriated to the optimization problem proposed here because PMU placement and communication aspect of the synchrophasor infrastructure are intimately connected. In this way, a dual optimization solution is implemented to address two simultaneous optimizations, PMU placement and CI design, in order to minimize the overall cost.

In the first stage of the algorithm, the VNS solves these simultaneous optimizations, allocating PMUs for a global cost optimization. The optimal CI design is guaranteed by Dijkstra Shortest path computations in the local search process. The second stage of the algorithm finds the shortest communication network that links all the PMUs previously found, which means finding a MCDS.

This MCDS is determined in a weighted graph, which is somewhat different from MCDS's problems in the literature that aim to minimize the CDS cardinality.

Therefore the proposed algorithm constructs a MCDS respecting this specificity.

Although there are two stages, the MCDS is predefined in the first stage and the second stage is only necessary for further processing Dijkstra shortest path calculations already done in the first stage. The two stages of the proposed algorithm are described as follows.

Stage I -Optimal PMU Placement

The objective of this stage is not only to minimize costs with PMUs but also to choose the best PMU configuration that can impact more effectively the CI optimization and globally achieve the lowest total cost. Following reference [START_REF] Meghanathan | Use of Centrality Metrics to Determine Connected Dominating Sets for Real-World Network Graphs[END_REF], The flow chart of this algorithm is represented in Figure 2.

which

Stage II -Communication Network Computation

After optimal PMU placement, a communication network interconnecting the buses with PMUs is computed. This is related to constructing a CDS from a pre-determined DS. 1. With the PMUs (DS vertices) and the single-shortest-paths previously calculated, the vertices which do not require PMUs are identified. A subgraph of the original network is obtained by removing these vertices.

2. Kruskal algorithm is applied on this subgraph to set a minimum spanning tree structure to the PMU communication network.

The flow chart of these steps is presented in Figure 3.

Use Kruskal Algorithm to set a Minimum Spanning Tree for G2.

CI End

From Dijkstra shortest path tree, find all vertices not necessary to connect buses with PMUs.

Derivate a new subgraph G2=(V' , E') by suppressing those vertices from the original graph.

Results and Discussion

This section shows results of the proposed algorithms for some of the most common test transmission networks, such as IEEE 14-bus, 24-bus, 30-bus, and 118-bus. Additionally, it shows results performed for IEEE 300-bus and in 5804-bus Brazilian transmission network in order to verify the scalability of the algorithms.

OPP Optimization Results

Results of simulations for an optimal PMU placement problem under normal conditions (without considering ZIB) are presented in Table 1. The optimal number of PMUs required for different systems and related bus locations are shown. The minimum number of PMUs allocated reached the best results found in the literature [START_REF] Roy | An optimal PMU placement technique for power system observability[END_REF], which is the optimal solution in this condition.

Simulations are also performed on larger networks such as IEEE 300-bus and Brazilian 5804-bus. For IEEE 300-bus, the method allocated 87 PMUs which is compared to the best result found in the literature [START_REF] Wen | Optimal multistage PMU placement for wide-area monitoring[END_REF]. Related to the Brazilian network, 34% of buses were equipped with PMUs which is almost at the same rate for IEEE 30-bus (33%). Simulations were also extended to PMU loss and optimal PMU placements considering ZIB (Table 4). The number of PMUs to assure a more robust monitoring system demands much more buses equipped with PMUs compared to normal conditions. Table 5 compares these results to other methods, showing good results in IEEE 14-bus and 24-bus networks but some less efficient results in the others (IEEE 30-bus, 57-bus and 118-bus).

Optimization of Total Cost

After validations of the first algorithm, simulations are also extended to total cost optimization problems in WAMS, and results are presented in this section. To be able to compare our results with related works, two cases (A and B) used in previous works are presented in Table 6, and used as input data for simulations. The same distance matrix is also assumed between buses from previous works, i.e., that all transmission lines have the same conductors with the same configurations and relative distances between system buses were extracted from the system admittance matrix [START_REF] Washington | Power system test case archive[END_REF]. For each case, the total length of the transmission lines is shown in the same table.

Case A

Simulations for this case are summarized in Table 7. It can be noticed from these results that the optimal PMU configurations which result in lower costs and assure system observability require in general more PMUs when compared to simple OPP algorithms (Table 1). Additionally, having more PMUs is better for redundancy measurement and state estimation. This conclusion is also stated in previous works. However, when compared to related works in Table 8,

Optimization of Total Cost in Large Networks

The algorithm was also tested on large networks, such as IEEE 300-bus and BR 5804-bus. Distances between buses in IEEE 300-bus were set following the same procedure stated for smaller IEEE networks in the previous section.

For the larger Brazilian Network BR5804-bus the distance matrix was built in MATLAB R using a random number generator (rgn) seeded to a unit value and a randi function with mean equal to fourteen in order to get a total network length of 100.000 km, which is near the real value [START_REF] Ons | Power and Transmission system of Brazil[END_REF]. Results obtained in Table 12 have confirmed the influence of communication cost on resulting number of PMUs, which is in general higher for cost optimization problems than for simple minimization OPP approaches (Table 2).

It must be stated that different PMU and communication equipment costs could lead to different results. This first algorithm proved to be flexible, efficient and scalable and in which a new graph theory metric proposed has played a fundamental role.

Subsequently, a total cost optimization problem was introduced and the previous algorithm was adapted and extended to solve this problem.

The algorithm proved its simplicity and flexibility in considering several operating conditions and real-life cost parameters for PMU and communication infrastructure. It was able to manage simulations in large and complex networks and resulted in good solutions compared to previous works.

Figure 1 :

 1 Figure 1: Centralized WAMS Architecture.

( 1 )

 1 Initial solution, (2) Shaking (or perturbation), which allows a progressive varying of the neighbourhood of the current solution; and (3) Local search (or descent), which is a local search procedure. The shaking stage is necessary to move the current solution away from a local optimum, and it is parameterized by the shaking amplitude k, which changes the current solution more and more as it increases. Each value of k corresponds to a different neighbourhood N k . The search can become more global and diverse when there is no progress in the solution. If the local search can

  vertices along with a set of |E| edges connecting those vertices. A graph with N = |V | vertices can be represented by an adjacency matrix A = [a ij ] N ×N , where a ij = 1 if there is an edge connecting vertex i to vertex j, and 0 if otherwise. A subgraph is a subset obtained by removing any set of vertices or edges from the original graph. The degree of a vertex i is the number of edges incident to i. Betweenness Centrality (BWC) or simply Betweenness is a measurement of how significant a vertex is in facilitating communication between any two vertices in the network by means of the shortest paths.

  showed the advantage of Betweenness Centrality in the determination of a MCDS, the present methodology proposes to combine both metrics, Dominance and Betweenness for PMU placement. Dominance is used to minimize the number of PMUs, while Betweenness Centrality seeks, at the same time, to influence the best PMU configuration that provides the lowest communication network costs. A VNS-based heuristic is also used, which is slightly different from the one proposed for minimizing the number of PMUs in Section 3.3. The operation of the VNS is detailed as follows: Previous Calculations: From matrix A, the degree and the Dominance values are computed for each vertex; Betweenness of each vertex is computed from network distance and cable unit cost. Dijkstra algorithm is used to compute the single-shortest-paths from the bus with higher Dominance value to all remaining buses. Initial Solution: 33% of the buses with the highest Dominance value for PMU allocation are selected. Additionally, PMUs are assigned to all buses that are adjacent to terminal buses. Finally, a check is performed to identify which buses are not yet observed by PMUs and in such cases PMUs are allocated to those buses. For single PMU outage and branch failure requirements, the initial solution is formed by placing PMUs in all buses in both cases. Shaking: The process to disturb the current solution is the same as described in Section 3.3.1. Local Search: The local search is performed after each shaking and seeks to reduce the number of PMUs. PMUs are sequentially eliminated from buses with lowest to highest Dominance or Betweenness. These two metrics are used alternately from one iteration to another. A PMU is only removed from a bus if it does not violate the observability constraint A • X ≥ b. The new solution is assumed whenever its related total cost is lower or at least the same from the previous solution. The total cost comprises costs with PMUs, switches and OPGW. This OPGW cost is determined using the shortest paths previously calculated. With a new solution a new perturbation restarts with k = 1.

Figure 2 :

 2 Figure 2: Flow chart of Stage I.

Figure 3 :

 3 Figure 3: Flow chart of Stage II.

  are assigned to all buses that are adjacent to terminal buses. Finally, a check is performed to identify which buses are not yet observed by PMUs and in such cases, PMUs are allocated to these buses.

	A • X ≥ b. The new solution is assumed whenever it has at least the same
	number of PMUs and, in this case, perturbation restarts with k = 1.
	3.3.2. Single PMU loss and branch failure conditions
	Supporting single PMU loss is one of the robustness requirements for a mon-
	itoring system. In this condition it is ensured that in case of any single PMU
	failure the system remains observable. To meet this requirement each bus must
	be observed (or covered) by two PMUs, i.e., be adjacent to two PMU buses, or by
	Local Search: The local search is performed after each shaking and seeks
	to reduce the number of PMUs. PMUs are sequentially eliminated from buses
	with lowest to highest Dominance or inversely in alternating iteration. A PMU
	is only removed from a bus if it does not violate the observability constraint

Shaking: This stage consists of disturbing the current solution and moving it away from local minima by allocating more PMUs to other buses. In this way a new solution is obtained which can be improved by a subsequent local search. The number of PMUs to be introduced is determined by a factor that is proportional to the size of the network and the selection of the buses is done in parameterized (k) and probabilistic ways. If the value of k is low, there is a greater probability that buses with higher Dominance values are selected to allocate new PMUs. If k increases, buses with lower Dominance value are selected, i.e., the solution (the neighbourhood) is increasingly varied. The increase of k always happens when the solution can not be improved in the previous local search iteration. When the local search can improve the solution, the disturbance resumes with k = 1. one when the bus itself has another PMU allocated to it. This requirement does not introduce substantial changes in the base algorithm (observability under normal conditions) previously seen. For the initial solution, PMUs are allocated on all buses to ensure observability for single PMU loss. The shaking stage remains the same from the previous algorithm but it is not performed in the first iteration (the initial solution does not allow any perturbation). The local search is performed in the same way but it only accepts PMUs removal if all buses remain covered (or observable) by at least two PMUs (A • X ≥ [2 2 2 . . . 2] T ).

Table 1 :

 1 Optimal PMU Placement Results for normal operating conditions without ZIB.PMUs to meet observability requirements. Table2summarizes the percentage of PMUs allocated in the networks versus the average degree of the graph, which is one of the indicators of graph density. It can be noticed that the Brazilian Network has the highest percentage of buses with PMUs compared to other instances but it is also the least dense network (average degree 2.44).

	System	Optimal PMU locations	No. PMUs
	IEEE-14	2,6,7,9	4
	IEEE-24	2,8,10,16,21,23,24	7
	IEEE-30	2,4,6,9,10,12,18,23,25,29	10
	IEEE-57	1,5,9,12,15,17,21,23,28,30, 36,40,44,48,49,52,56	17
	IEEE-118	1,5,9,12,15,17,20,23,28,30, 36,40,44,46,50,52,56,62, 63,68,71,75,77,80,85,86, 90,94,102,105,110,115	32
	The minimum percentage of buses equipped with PMUs needed for observ-
	ability depends on network topology. Generally, a dense graph (i.e., a graph
	whose number of edges is high compared to the number of vertices) needs less
	For single branch outage operating conditions, simulations on IEEE 30-bus
	and 118-bus were considered with several branch failures scenarios. Results are
	presented in

Table 3 .

 3

Table 2 :

 2 Minimum number of PMUs versus graph average degree.

	System	IEEE-30 IEEE-57 IEEE-118 IEEE-300 BR-5804
	No. PMUs	10	17	32	87	1987
	% bus with PMUs	33%	30%	27%	29%	34%
	Average degree	2.73	2.74	3.03	2.73	2.44

Table 3 :

 3 Optimal PMU Placement Results for branch outage.

	System	IEEE-30	IEEE-118
	Branch	(2,4) (2,6) (12,15)	(4,5) (11,12) (17,30) (48,49)
	outage	(21,22) (23,24)	(54,59) (56,59) (61,62) (70,74)
			1, 5, 9, 11, 12, 15, 17, 20, 22, 25,
	Optimal PMU locations	1, 5, 6, 9, 10, 12, 15, 19, 24, 25, 27	28, 30, 32, 36, 37, 41, 43, 46, 49, 52, 54, 56, 59, 62, 64, 68, 70, 71, 75, 77, 80, 85, 86, 90, 92, 96,
			100, 105, 110, 115
	No. PMUs	11	40

Table 4 :

 4 Optimal PMU Placement results for single PMU loss condition considering ZIB.

	System	Optimal PMU locations	No. PMUs
	IEEE-14	1, 2, 4, 6, 9, 11, 13	7
	IEEE-24	1, 2, 7, 8, 9, 10, 14, 16, 18, 19, 20, 21	12
	IEEE-30	2, 3, 4, 6, 7, 10, 12, 13, 15, 17, 18, 20, 24, 25, 27, 30	16
	IEEE-57	2, 3, 5, 7, 9, 10, 12, 15, 17, 19, 21, 23, 27, 28, 31, 35, 37, 38, 42, 46, 49, 50, 52, 53, 55, 56, 57	27
		1, 2, 6, 7, 8, 9, 12, 13, 15, 17, 18, 20, 21, 23, 24, 27,
		28, 29, 32, 34, 36, 37, 40, 41, 44, 45, 46, 49, 51, 52,
	IEEE-118	54, 56, 57, 60, 62, 66, 68, 70, 72, 75, 77, 79, 80, 83,	63
		85, 86, 87, 89, 90, 92, 94, 96, 100, 102, 105, 107,
		109, 110, 111, 112, 115, 117, 118

Table 5 :

 5 Comparison of optimal PMU placement results for single PMU loss condition considering ZIBs.

	Methods IEEE-14 IEEE-24 IEEE-30 IEEE-57 IEEE-118
	Proposed	7	12	16	27	63
	[28]	7	13	15	26	64
	[32]	7	N/A	14	22	61

Table 6 :

 6 Parameters considered in simulations.

Table 7 :

 7 Total Cost (PMU and CI) Optimization for Normal Operating Conditions -Case A.

	System	No. PMUs	PMUs	No. Switches	km (OPGW)	US$ (x10 6 )
	IEEE-14	5	4, 5, 6, 7, 9	6	199	2.2
	IEEE-30	10	2, 4, 6, 9, 10, 12, 15, 20, 25, 27,	11	625	6.7
			1, 4, 9, 15, 18, 21, 24,
	IEEE-57	18	28, 29, 31, 32, 36, 38,	35	1481.2	15.5
			39, 41, 46, 51, 54
			2, 5, 9, 11, 12, 17, 21,
			27, 30, 31, 32, 34, 37,
	IEEE-118	36	40, 45, 49, 50, 51, 54, 59, 65, 66, 68, 70, 71,	51	1991.5	21.4
			75, 77, 80, 83, 86, 89,
			92, 96, 100, 105, 110

more PMUs for normal conditions, which is better for measurements redundancy and state estimation process. The communication network and the number of nodes involved are smaller in our method.

Table 8 :

 8 Comparison of results under normal conditions -Case A.

	System	Methods	No. PMUs	No. Switches	km (OPGW)	US$ x10 6
	IEEE-30	Proposed	10	11	625	6.7
		[4]	10	15	804.6	8.4
	IEEE-118	Proposed	36	51	1991.5	21.4
		[4]	39	76	3012.6	31.7

Table 9 :

 9 Total Cost (PMU and CI) Optimization for Normal Operating Conditions, considering ZIB -Case B.

	System	No. PMUs	PMUs	No. Switches	km (OPGW)	US$ (x10 6 )
	IEEE-14	4	4, 5, 6, 9	5	199	1.0
	IEEE-30	7	2, 4, 10, 12, 15, 20, 27	10	546.3	2.5
			1, 4, 10, 15, 20, 23,
	IEEE-57	16	28, 29, 31, 32, 36,	33	1476.7	6.5
			39, 41, 47, 49, 54
			3, 8, 11, 12, 17, 21,
			27, 31, 32, 34, 37,
			40, 45, 49, 51, 54,
	IEEE-118	32	56, 61, 66, 70, 71,	49	1983.3	9.7
			75, 77, 80, 83, 86,
			89, 92, 96, 100,
			105, 110

Table 12 :

 12 Total cost optimization in large networks.A new methodology for total cost optimization in WAMS considering PMUs allocation and communication infrastructure costs was presented in this paper.Firstly, a simple optimal PMU placement algorithm based on Variable Neighbourhood Search metaheuristic was which has never been used before in related PMU placement problems.

	System	No. PMUs	No. Switches	km (OPGW)	US$ (x10) 6
	IEEE-300	91	152	7862.4	35.9
	BR-5804	2004	2887	18220	156.7
	6. Conclusion

Acknowledgment This work was partially supported by CAPES-Project 88881.062156/2014-01, CNPq-Projects 462477/2014-2 and 304564/2016-8, FAPES-80607594/17NIDA Project (Brazil) and by the Foundation HEC Montréal (Canada).

$ This work was partially supported by CAPES-Project 88881.062156/2014-01, CNPq-Projects 462477/2014-2 and 304564/2016-8 (Brazil) and by the Fondation HEC Montréal (Canada).

Parameters

Case ). The resulting communication network is also smaller.

Case B

This case considers different cost structure for PMU and communication media. Simulations were done for normal operating conditions and PMU loss requirement and are presented in Tables 9 and10.

The results confirm previous arguments that the number of PMUs for a lower overall cost does not necessarily correspond to a minimal number of PMUs.

Comparing these results to two related works with same simulation conditions (case B) in Table 11, it is possible to see that the proposed method is more effective to get a lower total cost for a meter planning system in IEEE 118-bus.

The cost obtained by the algorithm managed to lower costs by 14.2% under constraint of normal operating condition and 20.5% for PMU loss requirement when compared to previous related works. Yet, the proposed algorithm obtained


Romain Absil

Eglantine Camby

Alain Hertz
email: alain.hertz@gerad.ca

Hadrien Melot

Hadrien Mélot

A sharp lower bound on the number of non-equivalent colorings of graphs of order n and maximum degree n -3

Keywords: non-equivalent colorings, graphical Bell number, extremal theory, fixed maximum degree

teaching and research institutions in France or abroad, or from public or private research centers.

A sharp lower bound on the number of non-equivalent colorings of graphs of order n and maximum degree n -3

Introduction

A coloring of a graph G is the assignment of a color to each vertex of G such that no two adjacent vertices share the same color. A subset of vertices assigned to the same color is called a color class. Two colorings of G are said equivalent if they induce the same partition of the vertex set into color classes. The graphical Bell number B(G) is the number of non-equivalent colorings of G. This invariant has been studied by 1 several authors in the last few years [START_REF] Duncan | Bell and Stirling numbers for disjoint unions of graphs[END_REF][START_REF] Duncan | Bell and stirling numbers for graphs[END_REF][START_REF] Galvin | Stirling numbers of forests and cycles[END_REF][START_REF] Kereskényi-Balogh | Stirling numbers of the second kind and Bell numbers for graphs[END_REF]. It is related to the σ-polynomial of a graph G [START_REF] Brenti | Expansions of chromatic polynomials and log-concavity[END_REF][START_REF] Brenti | Location of zeros of chromatic and related polynomials of graphs[END_REF][START_REF] Korfhage | Sigma-polynomials and graph coloring[END_REF] defined as the polynomial in x such that the coefficient of x k is the number of non-equivalent colorings of G with exactly k non-empty color classes. Hence, B(G) is the value of the σ-polynomial at x = 1. B(G) is also related to the standard Bell number B n (sequence A000110 in OEIS [START_REF] Sloane | The on-line encyclopedia of integer sequences[END_REF]) that corresponds to the number of partitions of a set of n elements into non-empty subsets, and is thus obviously the same as the number of non-equivalent colorings of a graph of order n and without any edge.

Extremal properties of B(G) were recently studied in [START_REF] Hertz | Counting the number of non-equivalent vertex colorings of a graph[END_REF]. In particular, the authors give a sharp upper bound on B(G) for graphs of order n and any fixed maximum degree. Also, given two positive integers n and r with 0 ≤ r ≤ n -1, let G n r denote the class of graphs of order n and maximum degree r. A graph G ∈ G n r is said extremal if B(G) = min H∈G n r B(H). Extremal graphs in G n r are characterized in [START_REF] Hertz | Counting the number of non-equivalent vertex colorings of a graph[END_REF] for r = 1, 2, n-2 and n -1. We continue along the same path by considering the case r = n -3. In other words, we characterize the graphs with maximum degree n-3 that minimize B(G), and thus determine a sharp lower bound on the number of non-equivalent colorings B(G) for graphs G in G n n-3 . The results presented here were first conjectured with the help of the system GraPHedron [START_REF] Mélot | Facet defining inequalities among graph invariants: The system GraPHedron[END_REF] and later outputted again by Digenes [START_REF] Absil | Digenes: genetic algorithms to discover conjectures about directed and undirected graphs[END_REF]. Both of these tools are conjecture-making systems helping discovery in graph theory.

The paper is organized as follows. In the next section, we fix some notations, give additional definitions, and recall some basic properties of B(G). In Section 3, we give a new proof of the sharp lower bound on B(G) for graphs G in G n n-2 . It is based on similar arguments as those that we will use for proving our main theorem. The sharp lower bound on B(G) for graphs G of order n and maximum degree n -3 is established in Section 4, together with a characterization of the extremal graphs in G n n-3 .

Notations and basic properties

We refer to the book of Diestel [START_REF] Diestel | Graph Theory[END_REF] for basic notations and definitions in graph theory. Let G = (V, E) be a simple undirected graph, we denote by n = |V | and m = |E| the order and size of G, respectively. We write G H if two graphs G and H are isomorphic. For a subset W ⊆ V of vertices, we denote by

G[W ] the subgraph of G induced by W , while G \ W denotes the subgraph of G induced by V \ W (i.e., G \ W = G[V \ W ]). The disjoint union of two graphs G 1 = (V 1 , E 1 ) and G 2 = (V 2 , E 2 ) with V 1 ∩ V 2 = ∅ is the graph G 1 ∪ G 2 = (V 1 ∪ V 2 , E 1 ∪ E 2 ), while the join of G 1 and G 2 is the graph G 1 + G 2 = (V 1 ∪ V 2 , E) such that E contains all edges in E 1 ∪ E 2 as
well as all edges linking a vertex in V 1 with a vertex in V 2 . We denote by pG the disjoint union of p copies of G.

Let u and v be two vertices of a graph G = (V, E), we denote by G |uv the graph obtained by identifying (merging) vertices u and v, and removing the edge linking them, if any. If u and v are two adjacent vertices in G, we denote by G -uv the graph obtained by removing the edge that links u with v in G. Similarly, if u and v are two non-adjacent vertices in G, we denote by G + uv the graph obtained by adding an edge between u and v.

We write G for the graph complement of G, while K n denotes the complete graph of order n, P n the path of order n, and C n the cycle of order n. Also, we write K a,b for the complete bipartite graph where a and b are the cardinalities of the two sets of vertices of the bipartition, and S n denotes the star on n vertices, that is K 1,n-1 .

The

degree d(v) of a vertex v in G = (V, E) is the number of vertices adjacent to v. Vertex v is isolated if d(v) = 0, it is dominating (or universal ) if d(v) = |V | -1,
and ∆(G) denotes the maximum degree of G.

A clique in a graph G is a set of pairwise adjacent vertices, while a stable set is a set of pairwise non-adjacent vertices. A coloring of G is the assignment of a color to each vertex of G such that no two adjacent vertices share the same color. A color class is a set of vertices with the same color. Each color class is therefore a stable set, and every coloring of G induces a partition of its vertex set into color classes. Two colorings of G are equivalent if they induce the same partition of the vertex set into color classes. The graphical Bell number B(G) of G is the total number of non-equivalent colorings of G. As mentioned in Section 1, given two positive integers n and r with 0 ≤ r ≤ n -1, G n r denotes the class of graphs of order n with ∆(G) = r, and a graph G in

G n r is extremal if B(G) ≤ B(H) for all H ∈ G n r .
The following deletion-contraction rule (also called the Fundamental Reduction Theorem [START_REF] Dong | Chromatic polynomials and chromaticity of graphs[END_REF]) is a well known method to compute B(G) [START_REF] Duncan | Bell and stirling numbers for graphs[END_REF][START_REF] Kereskényi-Balogh | Stirling numbers of the second kind and Bell numbers for graphs[END_REF].

Property 1. [START_REF] Duncan | Bell and stirling numbers for graphs[END_REF][START_REF] Kereskényi-Balogh | Stirling numbers of the second kind and Bell numbers for graphs[END_REF] Let u and v be two vertices of a graph G: (

B(G) = B(G -uv) -B(G |uv ) if u and v are adjacent vertices (1)
) 2
The following property is a direct consequence of the fact that a dominating vertex in a graph G is alone in its color class for all colorings of G.

Property 2. [7, 9] If v is a dominating vertex of G, then B(G) = B(G \ {v}).
The next property is very useful for recognizing extremal graphs. The following simple property shows how to compute B(G) for a join

G = G 1 + G 2 .
It follows from the fact that none of the vertices of G 1 can share a color with a vertex of G 2 .

Property 4. If G = G 1 + G 2 then B(G) = B(G 1 )B(G 2 ).
Since the complement of a join is the disjoint union of two graphs, we have the following corollary.

Corollary 5. Let G be a graph, and let X 1 , X 2 , . . . , X k denote the vertex sets of the connected components of G. Then

B(G) = k i=1 B(G[X i ]).
The next two properties link the graphical Bell number with Fibonacci and Lucas numbers F n and L n (sequences A000045 and A000204 in OEIS [START_REF] Sloane | The on-line encyclopedia of integer sequences[END_REF]). They also give the values of the graphical Bell number of complements of paths and cycles, which will be very useful for the proofs contained in the next sections. Property 6. [START_REF] Duncan | Bell and stirling numbers for graphs[END_REF][START_REF] Kereskényi-Balogh | Stirling numbers of the second kind and Bell numbers for graphs[END_REF]

B(P n ) = F n+1 for all n ≥ 1. Property 7. [7, 10] B(C n ) = L n for all n ≥ 4. 3 Extremal graphs in G n n-2
The following theorem characterizes the extremal graphs in G n n-2 , and thus gives a sharp lower bound on B(G) for graphs of order n and maximum degree n -2. It was proved in [START_REF] Hertz | Counting the number of non-equivalent vertex colorings of a graph[END_REF], but with a different scheme than that proposed here below. The new proof is based on arguments similar to those we will use in the next section for proving Theorem 14 that characterizes extremal graphs in G n n-3 . Theorem 8. [START_REF] Hertz | Counting the number of non-equivalent vertex colorings of a graph[END_REF] Let G = (V, E) be a graph of order n ≥ 2 and maximum degree

∆(G) ≤ n -2. Then B(G) ≥ n with equality if and only if G K 1 ∪ K n-1 when n = 4, and G K 1 ∪ K 3 or G C 4 otherwise.
Proof. We can assume G ∈ G n n-2 since B(G) > B(G + uv) for every pair of non-adjacent vertices u, v in G. We proceed by induction on n. The result is clearly valid for n = 2 since G then contains two isolated vertices, which means that G K 1 ∪ K 1 , and B(G) = 2. So assume n > 2 and suppose G is extremal in G n n-2 . Notice first that B(K 1 ∪ K n-1 ) = n because either the isolated vertex of K 1 has its own color, or it uses one of the n -1 colors in K n-1 . Hence, B(G) ≤ n. Since B(C 4 ) = 4, it remains to prove that G K 1 ∪ K n-1 when n = 4, and G K 3 ∪ K 1 or G C 4 otherwise. Equivalently, we have to prove that G S n when n = 4, and G S 4 or 2K 2 otherwise.

Let A be the set of vertices of degree 1 in G and let

B = V \ A. Then G[A] is the disjoint union of edges and isolated vertices, that is G[A] = pK 2 ∪ qK 1 with 2p + q = |A|. Moreover, every edge in G[A] is a connected component of G. Also, since G is extremal, we know from Property 3 that G[B]
n ≥ B(G) = B(G[{u, v}])B(G \ {u, v}) ≥ 2(n -2).
We therefore have n = 4, which means that G 2K 2 .

• If p = 0 then G[A] is a stable set, which means that G is the disjoint union of stars S n 1 , . . . , S nr with r i=1 n i = n and n i ≥ 3 (i = 1, . . . , r). If r = 1 then G S n . So assume r > 1 and let W be the vertex set of S n 1 . It then follows from Corollary

5 that B(G) = B(G[W ])B(G \ W ). Since G[W ]
has n 1 vertices and maximum degree n 1 -2 while G \ W has n -n 1 vertices and maximum degree (n -n 1 ) -2, we get a contradiction, using the induction hypothesis:

n ≥ B(G) ≥ n 1 (n -n 1 ) = n + (n 1 -1)n -n 2 1 ≥ n + (n 1 -1)(n 1 + 3) -n 2 1 = n + 2n 1 -3 > n. 4 Extremal graphs in G n n-3
We begin this section with some simple observations about graphs in G n n-3 :

• if n = 3, then C 3 is the only graph in G n n-3 ; • if n ≥ 4, then both C n and K 2 ∪ K n-2 belong to G n n-3 .
Hence, for an extremal graph

G ∈ G n n-3 , we have B(G) = B(C 3 ) = 5 if n = 3, and B(G) ≤ min{B(C n ), B(K 2 ∪ K n-2 )} if n ≥ 4. We know from Property 7 that B(C n ) = L n for n ≥ 4. The next lemma gives the value of B(K 2 ∪ K n-2 ) for n ≥ 4. Lemma 9. If n ≥ 4, then B(K 2 ∪ K n-2 ) = n 2 -3n + 3.
Proof. Let u and v be the two vertices of the K 2 . If vertex u gets the same color as one vertex in K n-2 , then v can get any of the n -3 other colors used in K n-2 , or a new color, which gives a total of (n -2) 2 non-equivalent such colorings. If vertex u does not share a color with a vertex in K n-2 , then v can have any of the colors used in K n-2 , or a new color different from that used by u, which gives a total of n -1 non-equivalent colorings. Hence, B(K

2 ∪ K n-2 ) = (n -2) 2 + n -1 = n 2 -3n + 3. Since L n > n 2 -3n + 3 for n ≥ 8 while L n ≤ n 2 -3n + 3 for 4 ≤ n ≤ 7, we have shown that if G is an extremal graph in G n n-3 , then B(G) ≤ h(n) for all n ≥ 3
, where h(n) is defined as follows:

h(n) =    5 if n = 3 L n if 4 ≤ n ≤ 7 n 2 -3n + 3 if n ≥ 8.
We will show in Theorem 14 that B(G) = h(n) for all extremal graphs G ∈ G n n-3 . We first prove some properties of function h(n).

Lemma 10. Let n and r be two integers such that n ≥ r + 3 ≥ 5. Then

h(n) < F r+1 h(n -r) + 2F r (n -r -1).
Proof. Consider any fixed r ≥ 2, and let f (n) = F r+1 h(n -r) + 2F r (n -r -1). We distinguish the following cases.

• If n -r ≥ 8, then consider function g(n) = f (n) -h(n). It follows from the definition of h(n) that g(n) = F r+1 ((n -r) 2 -3(n -r) + 3) + 2F r (n -r -1) -n 2 + 3n -3.
The first derivative of g is

n(2F r+1 -2) + 2F r -2rF r+1 -3F r+1 + 3.
It is equal to zero when n = (2r + 3)F r+1 -2F r -3 2F r+1 -2 < r + 3, and strictly positive for larger values of n. Hence, for n ≥ r + 8, g(n) is an increasing function, and since g(r + 8) = 43F r+1 + 14F r -r 2 -13r -43 > 0 for all r > 1, we conclude that g(n) is strictly positive (i.e., h(n) < f (n)) for all n ≥ r + 8.

• If n -r = 7 then f (n) = 29F n-6 + 12F n-7 > n 2 -3n + 3 = h(n) for all n ≥ 9. • If n -r = 6 then f (n) = 18F n-5 + 10F n-6 > n 2 -3n + 3 = h(n) for all n ≥ 8. • If n -r = 5 then f (n) = 11F n-4 + 8F n-5 > n 2 -3n + 3 = h(n) for all n ≥ 8 and f (7) = 30 > 29 = h(7). • If n -r = 4 then f (n) = 7F n-3 + 6F n-4 > n 2 -3n + 3 ≥ h(n) for all n ≥ 7 and f (6) = 20 > 18 = h(6). • If n -r = 3 then f (n) = 5F n-2 + 4F n-3 > n 2 -3n + 3 ≥ h(n) for all n ≥ 5.
Lemma 11. Let n 1 , n 2 be two integers such that n 1 ≥ 3 and n 2 ≥ 3. Then

h(n 1 )h(n 2 ) > h(n 1 + n 2 ).
Proof. We analyze three cases.

• If both n 1 and n 2 belong to {3, 4, 5, 6, 7} then the result can be checked with the values in Table 1. • If n 1 ∈ {3, 4, 5, 6, 7} and n 2 ≥ 8, then the following inequalities (which are easy to verify) prove the result:

n
5(n 2 2 -3n 2 + 3) > (n 2 + 3) 2 -3(n 2 + 3) + 3, 7(n 2 2 -3n 2 + 3) > (n 2 + 4) 2 -3(n 2 + 4) + 3, 11(n 2 2 -3n 2 + 3) > (n 2 + 5) 2 -3(n 2 + 5) + 3, 18(n 2 2 -3n 2 + 3) > (n 2 + 6) 2 -3(n 2 + 6) + 3, 29(n 2 2 -3n 2 + 3) > (n 2 + 7) 2 -3(n 2 + 7) + 3.
• If both n 1 and n 2 are at least equal to 8, then let p(n) = n 2 -3n + 3 and assume without loss of generality that n 2 ≥ n 1 . We then have

h(n 1 )h(n 2 ) -h(n 1 + n 2 ) = p(n 1 )p(n 2 ) -p(n 1 + n 2 ) =(n 2 1 -3n 1 + 3)p(n 2 ) -((n 1 + n 2 ) 2 -3(n 1 + n 2 ) + 3) =n 1 [n 1 (p(n 2 ) -1) -(3p(n 2 ) + 2n 2 -3)] + (3p(n 2 ) -n 2 2 ) + (3n 2 -3).
Observe that the last two terms are strictly positive since 3p(n 2 ) -n 2 2 > 0 for all n 2 ≥ 4 and 3n 2 -3 > 0 for all n 2 ≥ 2. Hence,

h(n 1 )h(n 2 ) -h(n 1 + n 2 ) > n 1 (p(n 2 ) -1) -(3p(n 2 ) + 2n 2 -3) = (n 1 -3)p(n 2 ) -n 1 -2n 2 + 3 ≥ 5p(n 2 ) -3n 2 + 3 = 5n 2 2 -18n 2 + 18. Since 5n 2 2 -18n 2 + 18 > 0 for n 2 ≥ 0, we have h(n 1 )h(n 2 ) -h(n 1 + n 2 ) > 0.
Notation. Let v be a vertex in a graph G. For an integer r ≥ 1, we denote by G r v the graph obtained by adding an edge between v and one of the endpoint of P r in G ∪ P r (see Figure 1). For r = 0, we set

G 0 v = G. w r w 2 w 1 v G Figure 1: The graph G r v .
Lemma 12. Let v be a vertex in a graph G and consider any integer r ≥ 0. Then

B(G r v ) = F r+1 B(G) + F r B(G \ {v}).
Proof. If r = 0 the result is clearly valid since G 0 v = G, F 1 = 1 and F 0 = 0. So assume r ≥ 1. Let {w 1 , . . . , w r } be the vertex set of P r and let {w 1 w 2 , . . . , w r-1 w r } be its edge set so that G r v is obtained from G by adding the edge vw

1 to G ∪ P r . It follows from Equation (2) in Property 1 that B(G r v ) = B(G r v +w r w r-1 )+B(G r v|w r w r-1 )
, where w 0 = v. Observe that G r v + w r w r-1 is isomorphic to the graph obtained by adding a dominating vertex to G r v \ {w r } while G r v|w r w r-1 is isomorphic to the graph obtained by adding a dominating vertex to G r v \ {w r , w r-1 }. By Property 2, we therefore have

B(G r v ) = B(G r v \ {w r }) + B(G r v \ {w r , w r-1 }). If r = 1 then G r v \{w r } = G and G r v \{w r , w r-1 } = G\{v}, which means that B(G r v ) = F r+1 B(G)+F r B(G\{v}) (since F 2 = F 1 = 1)
. So assume r > 1 and suppose the result is valid for smaller values of r. Observe that G r v \{w r } G r-1 v and G r v \{w r , w r-1 } G r-2 v . It then follows from the induction hypothesis that

B(G r v ) = B(G r-1 v ) + B(G r-2 v ) = F r B(G) + F r-1 B(G \ {v}) + F r-1 B(G) + F r-2 B(G \ {v}) = F r+1 B(G) + F r B(G \ {v}).
Definition 1. Consider r disjoint paths P 1 , . . . , P r with endpoints u i , v i (i = 1, . . . , r). We call an eglantine the graph obtained by adding a vertex c to r i=1 P i and linking c to u i and v i (i = 1, . . . , r). Each subgraph induced by c and one of the paths is called a petal, and c is the center of the eglantine.

For illustration, all eglantines of order 9 are shown on Figure 2, with the middle vertex as center. Lemma 13. Let G be the complement of an eglantine of order n with r ≥ 1 petals. If all petals of G have exactly 3 vertices, then

B(G) = 2 r-1 (3r + 2).
Proof. The proof is by induction on r. If r = 1 then G C 3 and B(C 3 ) = 5 = 2 0 (3 + 2). So suppose r ≥ 2 and assume the result is true for smaller values of r. Let c be the center of G, and let u, v be the two other vertices of a petal in G. Also, let H denote the eglantine with r -1 petals obtained by removing u and v from G. 3). Hence, we know from Lemma 12 that B(G + cu) = 2B(H) + B(H \ {c}), where H \ {c} (r -1)K 2 . Moreover, by Corollary 5, we have B(rK 2 ) = 2 r and B((r -1)K 2 ) = 2 r-1 . In summary, it follows from the induction hypothesis that We are now ready to prove the main result.

It follows from Equation (2) in Property 1 that B(G) = B(G + cu) + B(G |cu ). Observe that G |cu is isomorphic to G \ {c} rK 2 while G + cu = G -cu H 2 c (see Figure
B(G) = 2(2 r-2 (3(r -1) + 2)) + 2 r-1 + 2 r = 2 r-1 (3r + 2).
Theorem 14. Let G = (V, E) be a graph of order n ≥ 3 and maximum degree

∆(G) ≤ n -3. Then B(G) ≥ h(n)
with equality if and only if G K n-2 ∪ K 2 when n ≥ 8, and G C n otherwise.

Proof. We can assume G ∈ G n n-3 since B(G) > B(G + uv) for every pair of non-adjacent vertices u, v in G. The proof is by induction on n. The theorem is clearly valid for n = 3 since C 3 is the only graph in G 3 0 and B(C

3 ) = 5 = h(3). So assume n ≥ 4 and suppose G is extremal in G n n-3 . By definition of function h(n), it is sufficient to prove that if G is not isomorphic to K 2 ∪ K n-2 or to C n , then B(G) > h(n). So assume G K 2 ∪ K n-2 and G C n .
Assume G has k ≥ 2 connected components. Let V 1 be the vertex set of one connected component of G and let

V 2 = V \V 1 , with n i = |V i |. Note that n i ≥ 3 because ∆(G) = n -3. Then G is the join of G[V 1 ] and G[V 2 ] and it follows from Property 4 that B(G) = B(G[V 1 ])B(G[V 2 ]). Moreover, G[V 1 ] ∈ G n 1 n 1 -3 and G[V 2 ] ∈ G n 2 n 2 -3
. It then follows from the induction hypothesis and Lemma 11 that B(G) ≥ h(n 1 )h(n 2 ) > h(n).

So assume G is connected, let A be the set of vertices of degree 2 in G and let B = V \ A. Also, let W 1 , . . . , W k be the vertex sets of the connected components of G[A]. Since G is extremal, we know from Property 3 that G[B] is a stable set, which means that in G, every vertex in B has at least three neighbors in A. Moreover, every subgraph G[W i ] of G is a cycle or a path, and if it is a cycle, it is a connected component of G. Hence, since G is connected and G C n , we have k ≥ 2 and every G[W i ] is a path.

If a set W i contains only one vertex v, then let x i and y i be the two neigbhors of v in G. If |W i | > 1, then let u i and v i be the endpoints of the path G[W i ]. In G, u i has one neigbhor in W i and another one in B, which we denote by x i . Also, we denote by y i the neighbor of v i in B (with possibly x i = y i ). We now distinguish two cases. Case 1: x i = y i for i = 1, . . . , k.

In such a case, every W i contains at least two vertices, and because G is connected, it is an eglantine with k ≥ 2 petals and a center c equal to x i and y i for all i = 1, . . . , k. If all petals have 3 vertices, then n = 2k + 1 and we know from Lemma 13 that B(G) = 2 k-1 (3k + 2). Since k ≥ 2 (because we assume n ≥ 4), it is easy to verify that 2

k-1 (3k + 2) > 4k 2 -2k + 1. Hence B(G) > 4k 2 -2k + 1 = (2k + 1) 2 -3(2k + 1) + 3 = n 2 -3n + 3 ≥ h(n).
So assume that at least one petal G[W i ∪ {c}] (1 ≤ i ≤ k) has at least 4 vertices. Let {c, w 1 , . . . , w r } (r ≥ 3) be the vertex set of this petal and let {cw 1 , w 1 w 2 , . . . , w r-1 w r , w r c} be its edge set. Also, let H be the eglantine obtained by removing w 1 , . . . , w r from G.

B(G + cw 1 ) = F r+1 B(H) + F r B(H \ {c})).
Also, G |cw 1 is isomorphic to the graph obtained by adding a dominating vertex to G \ {c, w 1 } P r-1 + H \ {c}. It then follows from Properties 2, 4 and 6 that

B(G |cw 1 ) = B(G \ {c, w 1 }) = F r B(H \ {c}).
In summary, we have

B(G) = F r+1 B(H) + 2F r B(H \ {c}).
Since H has n -r vertices and maximum degree n -r -3, we know by induction that B(H) ≥ h(n -r). Also, because B(H \ {c}) has n -r -1 vertices and maximum degree n -r -3, we know from Theorem 8 that B(H \ {c}) ≥ n -r -1. Hence,

B(G)) ≥ F r+1 h(n -r) + 2F r (n -r -1)
and it follows from Lemma 10 that B(G) > h(n).

Case 2: there is at least one index i such that x i = y i .

Let G[W i ] be a connected component of G[A] of maximum order among those for which x i = y i . Suppose W i contains r ≥ 1 vertices w 1 , . . . , w r so that u i = w 1 , v i = w r and w 1 w 2 , . . . , w r-1 w r are the edges in G[W i ]. Note that n ≥ 5 because |B| ≥ 2 and every vertex in B is linked to at least three vertices of A in G.

It follows from Equation (2) in Property 1 that B(G) = B(G

+ x i w 1 ) + B(G |x i w 1 ). By denoting H = G \ W i , we have G + x i w 1 = G -x i w 1 H r y i (see Figure 4 (b)
), and we know from Lemma 12 that

B(G + x i w 1 ) = F r+1 B(H) + F r B(H \ {y i }).
Moreover, because G |x i w 1 is isomorphic to the graph obtained by adding a dominating vertex to G \ {x i , w 1 }, it follows from Property 2 that B(G |x i w 1 ) = B(G \ {x i , w 1 }). By denoting H = H \ {x i }, we have G \ {x i , w 1 } H r-1 y i (see Figure 4 (b)), and we therefore know from Lemma 12 that B(G \ {x i , v 1 }) = F r B(H ) + F r-1 B(H \ {y i }).

Since H has n -r vertices while its maximum degree is n -r -3, we know by induction that B(H) ≥ h(n -r). Also, because H \ {y i } and H have n -r -1 vertices while their maximum degree is n-r-3, we know from Theorem 8 that both B(H \{y i }) and B(H ) are at least equal to n -r -1. Moreover, we clearly have B(H \ {y i }) ≥ 1. Altogether, we therefore get B(G) = B(G + x i w 1 ) + B(G |x i w 1 ) = F r+1 B(H) + F r B(H \ {y i }) + F r B(H ) + F r-1 B(H \ {y i }) ≥ F r+1 h(n -r) + 2F r (n -r -1) + F r-1 .

Since F r-1 ≥ 0, it follows from Lemma 10 that B(G) > h(n) if r ≥ 2. We can therefore assume r = 1, which means that W i = {w 1 }, H G \ {w 1 } and B(G) ≥ B(G \ {w 1 }) + 2(n -2) ≥ h(n -1) + 2(n -2).

At this point, it might be helpful to remind some of our assumptions: G K 2 ∪ K n-2 and G

C n ; n ≥ 5; G is connected; all connected components G[W 1 ], . . . , G[W k ] of G[A] are paths; if x j = y j , then W j contains only one vertex; in particular, W i = {w 1 }.

• If n ≥ 9, then B(G) ≥ (n-1) 2 -3(n-1)+3+2(n-2) = n 2 -3n+3 = h(n). In order to have B(G) = h(n), we must have B(G\{w 1 }) = h(n-1) = (n-1) 2 -3(n-1)+3, which means (by induction) that G \ {w 1 } is isomorphic to K 2 ∪ K n-3 , and G is therefore isomorphic to K 2 ∪ K n-2 , a contradiction. So B(G) > h(n).

• If n = 8, then it is easy to verify that there are only five graphs G that satisfy all our assumptions. Their complements are represented on Figure 5 and their number B(G) of non-equivalent colorings are respectively equal to 73, 64, 55, 57 and 59. All these numbers are strictly larger than B(K 2 ∪ K 6 ) = 43. n=7 n=8 n=9

  B(G) = B(G + uv) + B(G |uv ) if u and v are non-adjacent vertices.

Property 3 . [ 9 ]

 39 If u and v are two non-adjacent vertices in an extremal graph G, then max{d(u), d(v)} = ∆(G).

  is a stable set, which means that in G, every vertex in B has at least two neighbors in A. • If p ≥ 1, then consider two adjacent vertices u and v in G[A]. Note that G[{u, v}] is a connected component of G and B(G[{u, v}]) = 2. Moreover, n ≥ 4 because G \ {u, v} contains at least one vertex in A and another one in A ∪ B. Since G \ {u, v} has n -2 vertices and maximum degree n -4, we know by induction that B(G \ {u, v}) ≥ n -2. It then follows from Corollary 5 that

Figure 2 :

 2 Figure 2: The eglantines of order 9.

2 GFigure 3 :

 23 Figure 3: Illustration of the proof of Lemma 13.

Figure 4 :

 4 Figure 4: Illustration of the proof of Theorem 14.

  It follows from Equation (2) in Property 1 that B(G) = B(G + cw 1 ) + B(G |cw 1 ). Since G + cw 1 = G -cw 1 H r c (see Figure 4 (a)), we know from Lemma 12 that

Figure 7 :

 7 Figure 7: Extremal graphs G with maximum degree ∆(G) = n -4.

Table 1 :

 1 Values of h(n) for 3 ≤ n ≤ 14.

	3 4 5 6 7 8 9 10 11 12 13 14
	h(n) 5 7 11 18 29 43 57 73 91 111 133 157

Polytechnique

Montréal, funded by two grants from Université Libre de Bruxelles and Fédération Wallonie-Bruxelles.

• If n = 7, then it is easy to verify that there are only two graphs G that satisfy all our assumptions. Their complements are shown on Figure 6 and their number B(G) of non-equivalent colorings are respectively equal to 38 and 45. These two numbers are strictly larger than B(C 7 ) = 29. • If n = 6, then B(G) ≥ h(5) + 2(6 -2) = 19 > 18 = h(6).

• If n = 5, then B(G) ≥ h(4) + 2(5 -2) = 13 > 11 = h(5).

Concluding remarks

We have determined a sharp lower bound on the graphical Bell number B(G) of graphs G of order n and maximum degre n-3. We have also characterized the extremal graphs in G n n-3 . Similar results were obtained in [START_REF] Hertz | Counting the number of non-equivalent vertex colorings of a graph[END_REF] for graphs with maximum degree r = 1, 2, n-2 and n -1. It would be interesting to determine a sharp lower bound on B(G) for graphs in G n r with r in {3, 4, . . . , n -4}. The extremal graphs in this case do not seem to have a simple structure, as was the case for ∆(G) = 1, 2, n -3, n -2, n -1. Indeed, we have determined some of them by exhaustive enumeration, and as an example, we represent in Figure 7 the only graphs G of order n = 7, 8, 9 with minimum value B(G) when ∆(G) = n -4.

Notice also that several graphs with minimum value B(G) are non-connected. This is the case, for example for K 2 ∪ K n-2 in G n n-3 , and for

It would be interesting to determine such extremal graphs with the additional constraint that G must be connected. Also, it could be interesting to characterize the graphs G that minimize or maximize B(G) when the order and the size of G are fixed.


Rogelio Cecena

R&D Benjamin Boniface

Associate Professor, IMT Gilles Dusserre

Rogelio Ceceña

Mines Albi-Institu

Clément Ader

Process optimization for high-temperature aerospace applications

Nowadays, the aircraft parts that are installed in zones subjected to temperatures between 200°C and 400°C are almost entirely metallic. Alloys based on aluminium, titanium and steel dominate these zones while composites are a no-go for this temperature range. This is due to the limitation in operating temperatures of conventional thermosets as well as the lack of economic viability of ceramic-matrix composites at temperatures below 400°C. Simulation of a known process, combined with innovative materials, might overcome the rule and increase the percentage of composites used in the ''warm'' zones of aricraft.


Zhesheng Chen
email: zheshengchen@gmail.com

Zailan Zhang

Johan Biscaras

Abhay Shukla
email: abhay.shukla@upmc.fr.

High performance self-driven photodetector based on graphene/ InSe/MoS 2 vertical heterostructure

Keywords: self-driven, built-in potential, graphene/InSe/MoS 2, heterostructure, photodetector

ABSTRACT: Hybrid Van Der Waals heterostructures comprising ultrathin layers of different materials and offering the possibility of novel properties and unusual charge transport characteristics have become a reality in recent years. Here, we vertically stack graphene, a transition metal dichalcogenide and a III-VI semiconductor together and report a novel selfdriven photodetector based on graphene/InSe/MoS 2 heterostructure. The device shows rectifying and bipolar behavior. In the self-driven mode it exhibits high photoresponsivity (110 mA/W), fast photo-response (less than 1 ms) and high detectivity (over 10 10 Jones). It also shows ambient operational stability over one month of operation and nearly uniform photocurrent distribution because of the efficient electron-hole pair separation arising from the large built-in potential at the interface of MoS 2 and InSe. Our graphene/InSe/MoS 2 heterostructure holds the promise for novel self-driven optoelectronics based on III-VI/ transition metal dichalcogenide heterojuctions.

TOC:

Graphene

Introduction

Vertical Van Der Waals (VDW) heterostructures with different materials have attracted tremendous interest not only for their outstanding properties such as strong interlayer interaction and ultrafast charge transfer process but also in optoelectronic applications such as photodiodes, light emitting diodes (LED) and solar cells with self-driven and rectifying features [1][2][3][4][5][6][7][8] . The ultrathin vertical VDW heterostructure is compatible with gating and the large junction area makes for an efficient and tunable device. The most studied vertical VDW heterostructures are p-n junctions assembled from various transition metal dichalcogenides (TMDs). For example, vertical MoS 2 /WS 2 heterostructures show gate tunable diode like current rectifying behavior and can act as a photovoltaic cell or as a self-driven photodetector with photo-switching ratio exceeding 10 3 . Moreover the collection of photo-excited carriers can be enhanced by adding top and bottom graphene layers 9 . In addition, other vertical VDW heterostructures with p-n junctions such as GaTe/InSe, black phosphorus/MoS 2 and SnSe/MoS 2 have been reported. All of them show diode-like behavior and photovoltaic effect due to the type-II VDW heterojunctions [10][11][12][13] . The self-driven behavior in a photodetector is based on the photovoltaic effect and should in principle be seen in devices based on type-II VDW heterojunctions because of the built-in potential. However the final performance is also determined by the properties of the hybrid device such as interlayer coupling, charge transfer and carrier lifetime. Few devices have been reported where the junction is between materials with the same kind of majority charge carriers like for example WS 2 /MoS 2 with n-n type and Se/polyaniline with p-p type heterostructures 14,15

Results and discussions

In the VDW graphene/InSe/MoS 2 heterostructure fabrication process, few-layer MoS 2 (~ 6 nm) was firstly mechanically exfoliated onto SiO 2 /Si substrate and then multilayer InSe (~ 50 nm) was randomly dry transferred onto MoS 2 , forming InSe/MoS 2 heterostructures. The interface between InSe and MoS 2 is very clean thanks to the dry transfer process and has been confirmed by our previous study 16 . Graphene was prepared on glass using the anodic bonding method 19 and wedging transferred onto the InSe/MoS 2 followed by plasma etching to complete the graphene/InSe/MoS 2 heterostructure. The detailed fabrication process and corresponding optical images can be seen in supplementary information 1 and Figure S1.

The schematic diagram and optical image of graphene/InSe/MoS 2 heterostructure are shown in Figure 1 (a) and (b). The thickness of InSe is around 50 nm so as to maximize light absorption and the bottom MoS 2 layer is only 6 nm so as to optimize electrostatic doping by the back gate.

In Figure 1(c), we characterize the heterostructure using Raman spectroscopy and show that vibration modes corresponding to InSe and MoS 2 are found in the range of 100 cm -1 ~500 cm -1 .

Moreover, the significant modification of Raman intensity (~5 times less from both InSe and property in this two materials is consistent with earlier transport measurements (I ds -V gs ) 20,21 . The electron affinity of InSe (-4.9 eV) 22 is smaller than that of MoS 2 (-4.3 eV) 23 The schematic photodetector based on graphene/InSe/MoS 2 heterostructure is shown in Figure 2(a). It is worth pointing out that a Schottky barrier is known to exist at Au contacts on InSe 22,24 .

An immediate benefit of the graphene layer on the top of the InSe/MoS 2 heterostructure is the ohmic contact between graphene and InSe 25 . The voltage between the source and drain determines the device current, and the voltage between the source and the back gate changes the Fermi level of the MoS 2 layer by electrostatic doping, which further modifies the band bending and the built-in potential at the interface. We first characterize the graphene/InSe/MoS 2 device by photocurrent mapping with 532 nm focused laser illumination without bias voltage or back gate voltage as shown in Figure 2(b-c). The photocurrent is negative because photoelectrons flow from InSe (drain) to MoS 2 (source) under the influence of the built-in potential. It is uniformly -25 nA in the whole device (device area: ~180 μm 2 ) except the area surrounding the Au contact where it is higher. This is in contrast to the InSe/MoS 2 heterostructure without the graphene transport layer where the photocurrent vanishes rapidly away from the Au contact under the same illumination conditions, as shown in Figure S2. The non-uniform photocurrent distribution in the InSe/MoS 2 heterostructure and the nearly uniform distribution of the photocurrent in the graphene/InSe/MoS 2 heterostructure can be understood from the simplified device model shown in Figure 2 (d-e). In general, the total resistance of each device can be split up into several contributions and we will discuss them separately in the following. In a device based on the InSe/MoS 2 heterostructure as show in Figure 2 (d), the total resistance includes the contact resistance Au/InSe (R 1 ) and MoS 2 /Au (R 8 ), the InSe (R 2,3,4 ) and MoS 2 (R 5,6,7 ) sheet resistance and the resistance at the InSe/MoS 2 interface R D (we use a diode symbol here because of the rectifying behavior shown in Figure S3). However, InSe sheet resistance is several orders of magnitude higher than MoS 2 sheet resistance. During photocurrent mapping, the transport paths of photo-excited carriers differ depending on the spot where light is focused resulting in different measured device resistance. When light is focused close to the top Au contact, holes can move through InSe to the top Au contact easily because of the short distance, while electrons can also travel through MoS 2 to the other Au contact easily because of the low resistance. Thus the total resistance for current through the device can be schematized by R 1 + R D + R 5 +R 6 +R 7 + R 8 . When light is focused away from the Au contact, total resistance for the current through the device is much higher due to the difficult hole transport through InSe schematized by R 1 + R 2 + R 3 +R 4 +R D + R 8 . Thus the photocurrent decreases rapidly away from the top Au/MoS 2 contact. In contrast in the device based on graphene/InSe/MoS 2 heterostructure shown in Figure 2 (e), the sheet resistance of InSe is inconsequential because the graphene sheet short-circuits the InSe layer. In this case, photo-excited carriers travel through the InSe layer vertically and then laterally through the graphene layer because of the much higher mobility of graphene reducing carrier recombination and leading to a persisting photocurrent away from the Au contact. This is the second advantage of placing graphene on top of the InSe/MoS 2 heterostructure. The photocurrent is still not totally uniform in the whole graphene/InSe/MoS 2 device, with a somewhat higher current in the region close to the Au/MoS 2 contact because of the difference in the sheet resistance of graphene and MoS 2 . This however could be corrected by fabricating a symmetric graphene/InSe/MoS 2 /graphene heterostructure. We then measure drain-source current-voltage (I ds -V ds ) curves in the absence of back gate voltage (V gs =0 V) under different light illumination power as shown in Figure 3(a). The dark response shows diode-like behavior as seen from the insert of Figure 3(a), and under illumination the curves present the rectifying characteristics of a photodiode. The absolute photocurrent increases rapidly when the bias voltage assumes negative values and starts to saturate at a certain point depending on the illumination power. The photovoltaic effect is also observed in our device corresponding to a short-circuit current (I sc ) at zero V ds or an open-circuit voltage (V oc ) at zero I ds under illumination (Insert panel of Figure 2(a) and Figure S3). In the graphene/InSe/MoS 2 device, I sc and V oc are both higher than the values in the InSe/MoS 2 heterostructure at the same illumination power. For example, I sc and V oc are -4.8 nA and 220 mV

Au Au InSe MoS 2 R 1 R 2 R 3 R 4 R 5 R 6 R 7 R 8 R D Au Au InSe MoS 2 R' 1 R' 2 R' 3 R 5 R 6 R 7 R 8 R D Gr R' 4 R o (d) (e)
in graphene/InSe/MoS 2 device, three times higher than in the device without graphene (-1.5 nA and 70 mV) at the illumination power of 126 mW/cm 2 , as shown in Figure S3. The photovoltaic effect leads to a self-driven device (V ds =0 and V gs =0). Figure 3 (b) shows the photo-response of the graphene/InSe/MoS 2 device in the self-driven mode and 532 nm incident light with 1.26 mW/cm 2 illumination power. Our device shows fast photo-response (τ) with both rise time and fall time less than 1 ms, which is three orders of magnitude less than that in homostructural InSe and MoS 2 photodetectors 26,27 . Furthermore we note that non-zero V g can change the Fermi level in MoS 2 , increasing (V g >0) or decreasing (V g <0) the Schottky interface barrier between InSe and MoS 2 and thus effectively modifying the photocurrent, as shown in Figure 3 (c). Finally, detectivity (D*) is another key parameter that represents the ability of the photodetector to distinguish signal from noise:

, where is the photoresponsivity shown in Figure 4, A is the device area and is the dark current. In the self-driven mode, the D* in our graphene/InSe/MoS 2 heterostructure is as high as 1.08×10 10 Jones (cm Hz 1/2 /W) because of the high light on-off ratio caused by the large built-in potential between InSe and MoS 2 . The ambient stability of the device is an important criterion in applications of 2D materials.

The stability of our graphene encapsulated device was characterized by monitoring the photocurrent and photo-response as a function of time as shown in Figure 3 (d). No change was observed over a month-long period under ambient conditions in strong contrast to the instability of a similar InSe/MoS 2 device (50% decrease in photocurrent over the same period as shown in Figure S4). The difference indicates that graphene plays a key role on the protection of InSe, which is in good agreement with Raman characterization from our previous study 18 . This is the third benefit of placing graphene on top of the InSe/MoS 2 heterostructure. , where and are fitting parameters 18,28 . Under 1.26 mW/cm 2 illumination power of 532 nm light, R and EQE are calculated to be 110 mA/W and 25.7 % in the self-driven mode. Both R and EQE decrease as P increases and moreover, this variation follows the power law very well with the fitting exponents and . This demonstrates that this fitting function can also be applied in self-driven photodetectors. We also compare these values with the performance R and EQE of InSe/MoS 2 photodetectors with different areas in Figure (4).

In these devices, as mentioned earlier, the active area is only around the contacts so that R and EQE, which are normalized to area, decrease as the InSe/MoS 2 device area increase. However in the graphene/InSe/MoS 2 device, the whole detector area is active with a non-zero photocurrent as shown in Figure 2(c), which results in a much higher R and EQE compared to similar sized InSe/MoS 2 photodetectors. As explained earlier this is due to efficient transport of photo-excited carriers through monolayer graphene and significant reduction of electron-hole pair recombination during carrier transport. Compared to the InSe/MoS 2 device, the graphene/InSe/MoS 2 photodetector is active over the whole detector area with a nearly uniform photocurrent distribution and high operational stability under ambient conditions. It also has a high photoresponsivity of 110 mA/W, a fast photoresponse time less than 1 ms and high detectivity resulting from a large built-in potential at the InSe/MoS 2 interface. The performance of our device indicates the promise of graphene encapsulated VDW heterostructures with III-VI and TMD semiconductors and highlights the benefits of a new generation of self-driven photodetectors.

Acknowledgement

MoS 2 )

 2 in our graphene/InSe/MoS 2 heterostructure indicates interference effects and a good interface quality. In order to investigate the electronic band alignment between InSe and MoS 2 , we measured the valence band electronic structure by means of angle resolved photoemission spectroscopy (ARPES) around the Brillouin zone center (Γ point) using the He Iα resonance line as shown in Figure1(d)and (e). From the ARPES spectra of cleaved MoS 2 and InSe, one can deduce that the surface chemical potential is at 1.10 eV (MoS 2 ) and 0.98 eV (InSe) above the valence band maximum (Γ point) indicating intrinsic n-doping of these two semiconductors with band gaps of 1.33 eV (6 nm MoS 2 ) and 1.26 eV (bulk InSe) respectively. The high n-doping

  causing ~0.7 eV difference of Fermi levels between the two materials before band alignment. However, when InSe and MoS 2 are in contact, the Fermi levels must coincide at the interface and a large built-in potential is formed. The corresponding type II band alignment diagram of the InSe/MoS 2 heterostructure in Figure1(f) shows that under illumination photoelectrons flow from InSe to MoS 2 thus causing a "self-driven" photocurrent without any applied bias or gate voltage.

Figure 1 .

 1 Figure 1. (a) Vertical VDW stacking, (b) the corresponding optical image and (c) Raman spectra of graphene/InSe/MoS 2 heterostructure. (d-e) ARPES spectra of MoS 2 and InSe acquired with He Iα resonance line. (f) Band alignment between InSe and MoS 2 at interface.

Figure 2 .

 2 Figure 2. (a) Schematic of graphene/InSe/MoS 2 photodetector. (b-c) Nearly uniform photocurrent distribution in graphene/InSe MoS 2 heterostructure. The red dash square represents the corresponding photocurrent mapping area. (d-e) Schematics of the simplified device models based on InSe/MoS 2 and graphene/InSe/MoS 2 heterostructures.

Figure 3 .

 3 Figure 3. (a) I ds -V ds characterization of graphene/InSe/MoS 2 heterostructure as a function of illumination power. Insert: focused dark current and I sc . (b) Rise (fall) time when light is turned on (off) from dark (illumination) status. (c) Photocurrent variation according to different back gate voltages. (d) Stability of graphene/InSe MoS 2 photodetector under ambient conditions and illumination power of 1.26 mW/cm 2 .

Figure 4

 4 Figure 4 shows photoresponsivity (, where I ph is photocurrent and P is

Figure 4 .

 4 Figure 4. R and EQE as a function of illumination power for graphene/InSe/MoS 2 heterostructure and compared with different sized devices of InSe/MoS 2 heterostructures. Symbols: experiment data, red line: parameterized fit of phenomenological function (see text).

  , in which similar rectifying properties and bipolar behavior can also exist. In our previous study, we combined n-type InSe with n-type MoS 2 and studied light interaction within vertical VDW

	InSe/MoS 2 heterostructures 16 . We found significant Raman intensity modification due to strong
	light coupling in the heterostrucure with a high quality interface and optimized light harvesting
	by controlling the layer thicknesses. In this work we add a monolayer graphene layer on top and
	fabricate a vertical graphene/InSe/MoS 2 heterostructure procuring three benefits: (1) ohmic
	contact between graphene and InSe; (2) efficient charge transport in graphene reducing carrier
	recombination 17 ; (3) encapsulation of InSe by covering graphene layer thus enhancing the
	stability of the device 18 . Our graphene/InSe/MoS 2 device shows current rectifying and bipolar
	behavior and functions as an air-stable self-driven photodetector with nearly uniform
	photocurrent distribution, fast photo-response speed and high detectivity at zero bias and gate
	voltage.

Table 1

 1 The detectivity in the device is two orders of magnitude higher than that of the MoTe 2 /MoS 2 device. Again, we ascribe this to the large built-in potential at the InSe/MoS 2 interface and the high carrier mobility in the graphene layer.

summarizes the key parameters of our graphene/InSe/MoS 2 photodetector, and compares them to some InSe based and MoS 2 based vertical VDW heterostructure photodetectors in the self-driven mode published in literature. The R and τ in our graphene/InSe/MoS 2 device are comparable to that of the previously reported state-of-the-art devices.

Table 1 .

 1 Comparison of key parameters between this work and some InSe based and MoS 2 based vertical VDW heterostructure photodetectors under self-driven mode. Only one other work calculated the detectivity explicitly. In summary, we fabricate a graphene/InSe/MoS 2 VDW heterostructure photodetector by random dry transfer of multilayer InSe on few-layer MoS 2 and wedging transfer of monolayer graphene on top. ARPES measurements show that both InSe and MoS 2 are highly n-doped and a type-II heterostructure is formed. In the self-driven mode, R and EQE are calculated at 110 mA/W and 25.65% which are comparable to other InSe based or MoS 2 based heterostructures.

	Devices	thickness (nm) R (mA/W) D* (Jones)	τ(ms)	light intensity (mW/cm 2 )	Ref
	GaTe/MoS 2	4.5/3.5	1370	-	10	100	29
	MoTe 2 /MoS 2	3.3/7	46	1.06 ×10 8	0.06	~73	30
	GaTe/InSe	18/10	13.8	-	0.02	0.16	10
	MoS 2 /perovskite	0.8/-	60	-	>200	0.7	31
	Graphene/InSe/MoS 2	0.34/40/6	110	1.08 ×10 10	<1	1.26	This work
	4. Conclusions

This work was partly supported by French ANR project (Grant No. ANR-11-BS04-0019), by LABEX MATISSE through a post-doctoral grant, by "Investissement d'Avenir" Labex PLAM (ANR-10-LABX-0039-PALM) and a scholarship from China Scholarship Council (CSC, Grant No. 2011618135). We acknowledge the "Consortium des salles blanches d'Ile de France" for access to clean room facilities. We acknowledge M Rosticher and J Palomo for help with clean room work, and A. Chevy for the precursor materials used in this work.

Supporting Information

Fabrication process of graphene/InSe/MoS 2 heterostructure; photocurrent distribution in InSe/MoS 2 heterostructure; I ds -V ds measurement in InSe/MoS 2 photodetector and stability of InSe/MoS 2 photodetector under ambient condition.