Thomas Ropars
email: thomas.ropars@inria.fr

Amina Guermouche

Bora Uçar
email: bora.ucar@ens-lyon.fr

Esteban Meneses
email: emenese2@illinois.edu

Laxmikant V Kalé
email: kale@illinois.edu

Franck Cappello
email: cappello@illinois.edu

On the Use of Cluster-Based Partial Message Logging to Improve Fault Tolerance for MPI HPC Applications

Fault tolerance is becoming a major concern in HPC systems. The two traditional approaches for message passing applications, coordinated checkpointing and message logging, have severe scalability issues. Coordinated checkpointing protocols make all processes roll back after a failure. Message logging protocols log a huge amount of data and can induce an overhead on communication performance. Hierarchical rollback-recovery protocols based on the combination of coordinated checkpointing and message logging are an alternative. These partial message logging protocols are based on process clustering: only messages between clusters are logged to limit the consequence of a failure to one cluster. These protocols would work efficiently only if one can find clusters of processes in the applications such that the ratio of logged messages is very low. We study the communication patterns of message passing HPC applications to show that partial message logging is suitable in most cases. We propose a partitioning algorithm to find suitable clusters of processes given the communication pattern of an application. Finally, we evaluate the efficiency of partial message logging using two state of the art protocols on a set of representative applications.

Introduction

The generation of HPC systems envisioned for 2018-2020 will reach Exascale from 100s of millions of cores. At such scale, failures cannot be considered as rare anymore and fault tolerance mechanisms are needed to ensure the successful termination of the applications. In this paper, we focus on message passing (MPI) HPC applications. For such applications, fault tolerance is usually provided through rollback-recovery techniques [START_REF] Elnozahy | A Survey of Rollback-Recovery Protocols in Message-Passing Systems[END_REF]: the state of the application processes is saved periodically in a checkpoint on a reliable storage to avoid restarting from the beginning in the event of a failure. In most cases, rollbackrecovery protocols are used flat: the same protocol is executed for all processes.

Flat rollback-recovery protocols have several drawbacks. Coordinated checkpointing requires to restart all processes in the event of a failure leading to a massive waste of resources and energy. Message logging protocols need to log all messages contents and delivery order, which leads to high storage resource occupation, high communication overhead and high energy consumption.

One way to cope with these limitations is to use hierarchical rollback-recovery protocols [START_REF] Guermouche | Uncoordinated Checkpointing Without Domino Effect for Send-Deterministic Message Passing Applications[END_REF][START_REF] Ho | Scalable Group-Based Checkpoint/Restart for Large-Scale Message-Passing Systems[END_REF][START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF][START_REF] Monnet | Hybrid Checkpointing for Parallel Applications in Cluster Federations[END_REF][START_REF] Yang | Trading Off Logging Overhead and Coordinating Overhead to Achieve Efficient Rollback Recovery[END_REF]. Clusters of processes are defined during the execution and different protocols are used inside and between clusters, giving a hierarchical aspect to the protocol. Typically, a partial message logging protocol logs messages between clusters to confine the effects of a process failure to one cluster and a coordinated checkpointing protocol is used within clusters [START_REF] Yang | Trading Off Logging Overhead and Coordinating Overhead to Achieve Efficient Rollback Recovery[END_REF][START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF]. The efficiency of these protocols depends on two conflicting requirements: i) the size of the clusters should be small to limit the rollbacks in the event of a failure; ii) the volume of inter-cluster messages should be low to limit the impact of logging on failure free performance.

This paper provides three main contributions: i) it shows that suitable process clustering can be found in most applications; ii) it proposes a bisection-based partitioning algorithm to find such clusters in an application based on its execution communication pattern; iii) it shows that partial message logging limits the amount of computing resources wasted for failure management compared to flat rollback-recovery protocols.

The paper is organized as follows. Section 2 details the context of this work, and presents the related work on MPI applications communications analysis. Section 3 analyzes a set of execution communication patterns in MPICH2. Section 4 presents our bisection-based partitioning algorithm, designed to address the partitioning problem for a partial message logging protocol. Using this algorithm, we evaluate the performance of two state-of-the-art partial message logging protocols [START_REF] Guermouche | Uncoordinated Checkpointing Without Domino Effect for Send-Deterministic Message Passing Applications[END_REF][START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF] on a set of representative MPI HPC applications. Results are presented in Section 5. Finally, conclusions are detailed in Section 6.

Context

In this section, we first present existing hierarchical rollback-recovery protocols. Then we present the applications studied in this paper. Finally, we detail the related work on analyzing characteristics of MPI HPC applications.

Hierarchical Rollback-Recovery Protocols Rollback-recovery techniques are based on saving information during the execution of an application to avoid restarting it from the beginning in the event of a failure. One of the main concerns is the amount of resources wasted with respect to computing power or energy. Several factors are contributing to this waste of resources: i) the overhead on performance during failure free execution; ii) the amount of storage resources used to save data; iii) the amount of computation rolled back after a failure.

Rollback-recovery protocols are usually divided into two categories: checkpointing-based and logging-based protocols [START_REF] Elnozahy | A Survey of Rollback-Recovery Protocols in Message-Passing Systems[END_REF]. In checkpointing protocols, pro-cesses checkpoints can be either coordinated at checkpoint time, taken independently or induced by the communications. For all these protocols, a single failure implies the rollback of all processes in most cases, which is a big waste of resources. On the other hand, message logging protocols log the content as well as the delivery order (determinant) of the messages exchanged during the execution of the application to be able to restart only the failed processes after a failure. However, logging all messages during a failure free execution can be very wasteful regarding communications and storage resources.

Hierarchical rollback-recovery protocols divide the processes of the applications into clusters and apply different protocols for the communications inside a cluster and for the communications among clusters. Our work focuses on partial message logging protocols [START_REF] Guermouche | Uncoordinated Checkpointing Without Domino Effect for Send-Deterministic Message Passing Applications[END_REF][START_REF] Ho | Scalable Group-Based Checkpoint/Restart for Large-Scale Message-Passing Systems[END_REF][START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF][START_REF] Yang | Trading Off Logging Overhead and Coordinating Overhead to Achieve Efficient Rollback Recovery[END_REF] which apply a checkpointing protocol inside the clusters and a message logging protocol among the clusters. They are attractive at large scale because only one cluster has to rollback in the event of a single failure. These protocols can work efficiently if the volume of inter-cluster messages is very low in which case the cost of message logging is small.

In this paper, we use two of these protocols for evaluations. Meneses et al. [START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF] use a coordinated checkpointing protocol inside clusters. Intra-cluster messages determinants are logged to be able to replay these messages in the same order after the failure and reach a consistent state. Considering a single failure, determinants can be logged in memory. In this study [START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF] is comparable to [START_REF] Ho | Scalable Group-Based Checkpoint/Restart for Large-Scale Message-Passing Systems[END_REF] and [START_REF] Yang | Trading Off Logging Overhead and Coordinating Overhead to Achieve Efficient Rollback Recovery[END_REF]. Guermouche et al. [START_REF] Guermouche | Uncoordinated Checkpointing Without Domino Effect for Send-Deterministic Message Passing Applications[END_REF] propose an uncoordinated checkpointing protocol without domino effect, relying on the send-determinism of MPI HPC applications. Using this protocol, an ordered set of p clusters can be defined. Only messages going from one cluster to a higher cluster are logged, limiting the number of clusters to roll back after a failure to (p + 1)/2 on average. Thanks to send-determinism, this protocol does not require any determinant to be logged and can tolerate multiple concurrent failures.

Applications Studied We study a representative set of MPI HPC applications to see if partial message logging could be used. Thirteen dwarfs have been defined and seven of them represent seven main classes of computation and communication patterns corresponding to numerical methods for high-end simulation in the physical sciences [START_REF] Asanovic | The Landscape of Parallel Computing Research: A View from Berkeley[END_REF]: dense linear algebra, sparse linear algebra, spectral methods, N-body methods, structured grids, unstructured grids, and MapReduce. This paper does not consider MapReduce applications because rollback-recovery is not adapted in this case. Our set of applications includes five of the NAS Parallel Benchmarks (NPB) [START_REF] Bailey | The NAS Parallel Benchmarks 2.0[END_REF] containing BT, LU, CG, FT, and MG; three of NERSC-6 Benchmarks [START_REF] Antypas | NERSC-6 Workload Analysis and Benchmark Selection Process[END_REF] (GTC, MAESTRO, and PARATEC); one of the Sequoia Benchmarks, http://asc.llnl.gov/sequoia/benchmarks/, (LAMMPS); and an Nbody kernel. Table 1 summarizes the dwarfs covered.

HPC Applications Communications Characteristics

The communication patterns of most of the applications considered in this paper have already been published [START_REF] Antypas | NERSC-6 Workload Analysis and Benchmark Selection Process[END_REF][START_REF] Riesen | Communication Patterns[END_REF]. Previous studies highlighted some properties. First, most MPI 1. Dwarfs covered by the studied applications applications make use of collective communications, but in general with a very small payload size that remains invariant with respect to the problem size [START_REF] Vetter | Communication Characteristics of Large-Scale Scientific Applications for Contemporary Cluster Architectures[END_REF]. Second, the communication graphs of many MPI HPC applications have a low degree of connectivity [START_REF] Kamil | Understanding ultra-scale application communication requirements[END_REF], which might indicate that processes can be partitioned into clusters.

Communication Patterns

In this section, we first study the communication patterns of some collective operations in MPICH2. Then, we present communication patterns we collected by running some applications. To get the communication patterns of MPI applications execution, we modified the code of MPICH2 1 to collect data on communications. The applications run on Rennes Grid'5000 cluster over TCP.

We focus on collective communications because they could generate patterns that are hard to cluster (they involve all processes in the application). Figure 1 presents the set of communication patterns used for MPICH2 collective communications, for a power-of-two number of processes (64 processes) and short messages (4 bytes). Details on the implementation of collective communications in MPICH2 can be found in [START_REF] Thakur | Optimization of Collective Communication Operations in MPICH[END_REF]. The recursive doubling algorithm, Fig. 1(a), is used to implement MPI Allgather and MPI Allreduce operations. The recursive halving algorithm, used in MPI Reduce scatter, has the same communication pattern. A binomial tree, Fig. 1(b), is used in MPI Bcast, MPI Reduce, MPI Gather and MPI Scatter. These two patterns are easily clusterizable using for instance clusters of size 16. The last pattern, corresponding to a store-andforward algorithm, is the one used for MPI Alltoall. This pattern is more difficult to cluster, but the use of MPI Alltoall should be limited in applications targeting very large scale. For large messages, many of the collective operations do involve communications between all application processes. Clustering is also difficult to apply to such patterns. However, as mentioned in Section 2, the payload for collective communications is usually small: on the set of applications we tested, we only found this pattern in NPB FT.

Although some of the collective communication patterns define natural clusters (see Fig. 1(a)), some others do not (see Fig. 1(c)). Furthermore, the pointto-point communications can render the patterns more sophisticated, so much so that the size or the number of clusters cannot be known a priori. Figure 2 presents the communication pattern of MAESTRO and GTC executed on 512 processes. As exemplified by these two applications, an automated means of clustering is strictly required to identify clusters of processes in an application. Here, we propose a bisection-based partitioning algorithm to automate clustering for partial message logging protocols. The objectives of the clustering are to reduce the inter-cluster communications, to increase the number of clusters, and to limit the maximum size of a cluster. We first illustrate the limits of existing tools and then present our new method.

Two Possible Approaches

A simplified version of the partitioning problem in which the objectives are to minimize the size of the logged messages and the maximum size of a part corresponds to the NP-complete graph partitioning problem (see the problem ND14 in [START_REF] Garey | Computers and Intractability; A Guide to the Theory of NP-Completeness[END_REF]). This can be easily seen by considering a graph whose vertices represent the processes and whose edges represent the communication between the corresponding two processes. Using heuristics for the graph partitioning problem would require knowing the maximum part size. This can be done but requires an insight into the application and the target machine architecture. A common variant of the above graph partitioning problem specifies the number of parts (in other words, specifies the average size of a part) and requires parts to have similar sizes (thusly reducing the maximum size of a partition). The problem remains NP-complete [START_REF] Bui | Finding good approximate vertex and edge partitions is NP-hard[END_REF]. Tools such as MeTiS [START_REF] Karypis | MeTiS: A Software Package for Partitioning Unstructured Graphs, Partitioning Meshes, and Computing Fill-Reducing Orderings of Sparse Matrices Version 4[END_REF] and Scotch [START_REF] Pellegrini | SCOTCH 5.1 User's Guide[END_REF] can be used to solve this problem. A possible but not an economical way to use those tools in our problem is to partition the processes for different number of parts (say 2, 4, 8, . . .) and try to select the best partition encountered.

Bisection-Based Partitioning

Bisection based algorithms are used recursively in graph and hypergraph partitioning tools, including PaToH [START_REF] Aykanat | PaToH: A multilevel hypergraph partitioning tool[END_REF], MeTiS and Scotch, to partition the input graph or hypergraph into a given number of parts. Simply put, for a given number K of parts, this approach divides the original graph or hypergraph into two almost equally sized parts and then recursively partitions each of the two parts into K/2 parts until the desired number of parts is obtained.

We adapt the bisection based approach and propose a few add-ons to address our partitioning problem. The proposed algorithm is seen in Algorithm 1. The algorithm accepts a set of P processes and a matrix M of size P × P representing the communications between the processes where M (u, v) is the volume of messages sent from the process u to the process v. The algorithm returns the number of parts K and the partition Π = P 1 , P 2 , . . . , P K . In the algorithm, the operation Π ← Π P k ⊕ P k1 ⊕ P k2 removes the part P k from the partition Π and adds P k1 and P k2 as new parts, where P k = P k1 ∪ P k2 . We use M (U, V to represent the messages from processes in the set U to those in the set V .

Algorithm 1 The proposed partitioning algorithm

Input: The set of P processes = {p1, p2, . . . , pP }; the P × P matrix M representing the communications between the processes Output: K : the number of process parts; Π = P1, P2, . . . , P K : a partition of processes 1: The algorithm uses bisection to partition the given number of processes into an unknown number of parts. Initially, all the processes are in a single part. At every step, the largest part is partitioned, if it should be, into two (by the subroutine Bisect), and then if the bisection is acceptable (determined by the subroutine AcceptBisection), that largest part is replaced by the two parts resulting from the bisection. If this new partition is the best one seen so far (tested in subroutine BestSoFar), it is saved as a possible output. Then, the algorithm proceeds to another step to pick the largest part. We now add a few details. Let P k be the largest part, then if P k should be partitioned (tested in the subroutine ShouldPartition), it is bisected into two P k1 and P k2 and tested for acceptance; if not accepted then the bisection is discarded and P k remains as is throughout the algorithm. Notice that when a bisection of P k is accepted, then P k is partitioned for good; if any BestSoFar test after this bisection returns true, then P k will not be in the output Π .

K ← K ← 1; B ← B ← 0 2: Π ← Π ← P1 = {p1, p2, . . . ,
The computational core of the algorithm is the Bisect routine. This routine accepts a set of processes and the communication between them and tries to partition the given set of processes into two almost equally sized parts by using existing tools straightforwardly. One can use MeTiS or Scotch quite effectively if the communications are bidirectional (M (u, v) = 0 =⇒ M (v, u) = 0), in which case M +M T 2 can be used. Alternatively one can use PaToH as each communication (bidirectional or not) can be uniquely represented as an edge.

The routines ShouldPartition, AcceptBisection, and BestSoFar form the essence of the algorithm. The routine BestSoFar requires a cost function to evaluate a partition. The cost function should be defined based on the metric the user wants to optimize, e.g., the performance overhead, and the characteristics of the targeted partial message logging protocol. It is function of the part sizes and of the volume of inter-parts messages. We define a cost function for both of the mentioned rollback-recovery protocols later in Section 5.1.

ShouldPartition is used to stop partitioning very small parts. It returns false if the size of the part in question is smaller than a threshold. Although we mostly used 1 as threshold in our experiments, using a larger threshold will make the algorithm faster. One could select this threshold based on a minimal part size considering the properties of the target machine architecture.

The routine AcceptBisection returns true in two cases. The first one is simply that for Π and Π = Π P k ⊕ P k1 ⊕ P k2 , we have cost(Π ) ≤ cost(Π) according to the cost function. If this does not hold, we may still get better partitions with further bisections. To this end, we have adapted the graph strength formula [START_REF] Cunningham | Optimal attack and reinforcement of a network[END_REF]. For a partition Π = P 1 , P 2 , . . . , P K , we use strength(Π) = B K-1 as its strength, B being the volume of inter-parts data. We accept the bisection if strength(Π ) ≤ strength(Π). If the bisection reduces the strength, then it can be beneficial, as it increases the size of the logged messages only a little with respect to the number of parts.

Evaluation

This section presents our experimental results. For the evaluations, we consider the two partial message logging protocols described in Section 2. We start by defining a cost function for each of them, evaluating the amount of wasted computing resources for failure management. Then we present the results obtained by running our implementation of the proposed partitioning algorithm using these cost functions, on the set of applications executions data we collected. We first show that these results overcome coordinated checkpointing and message logging protocols. Then we validate our partitioning algorithm by comparing our results to the results obtained by the existing graph partitioning tools.

Defining a Cost Function

We define a cost function to evaluate the amount of computing resources wasted by failure management, for the two partial message logging protocols. We would like to stress that the cost functions we define in this section do not give a very precise evaluation of the protocols, because a lot of parameters would have to be taken into account. We use some parameters, that we consider realistic, to provide an insight of the protocols cost. The main goal in this section, is to show that we manage to find clustering configurations with good trade-off between the clusters size and the amount of data logged.

We model the cost of a partial message logging protocol for a given clustering configuration as a function of the total volume of logged messages and the size of the clusters. We use a formula of the form

cost(Π) = α × L + β × R (1)
where L is the ratio of logged data and R is the ratio of processes to restart after a failure. The multipliers α and β are the cost associated with message logging and with restarting processes after a failure, respectively.

To evaluate the overhead of message logging, we use results from [START_REF] Guermouche | Uncoordinated Checkpointing Without Domino Effect for Send-Deterministic Message Passing Applications[END_REF], where message logging impact on communications (latency and bandwidth) on a high performance network results in a 23% performance drop on average. Therefore, α = 23% can be considered as a maximum theoretical overhead induced by message logging for a communication-bounded application. For the sake of simplicity, we do not consider in this study the cost of logging intra-cluster messages reception order in Meneses et al. protocol [START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF].

To evaluate the amount of computing resources wasted during recovery after a failure, we consider the global performance of the system. While a subset of the processes are recovering from a failure, the other processes usually have to wait for them to progress. We assume that the resources of the hanging processes could be temporarily allocated for other computations, until all application processes are ready to resume normal execution. So β includes only the amount of resources wasted by rolling back processes after a failure.

To compute β, we consider an execution scenario with the following parameters [START_REF] Cappello | Fault tolerance in petascale/ exascale systems: Current knowledge, challenges and research opportunities[END_REF]: a Mean Time Between Failure (M T BF ) of 1 day; 30 minutes to checkpoint the whole application (C); 30 minutes to restart the whole application (Rs). The optimum checkpoint interval (I) can be computed using Daly's formula [START_REF] Daly | A model for predicting the optimum checkpoint interval for restart dumps[END_REF]: I = (2 × C × (M T BF + Rs)) = 297min. This formula was originally used in coordinated checkpointing protocols, but we think we can safely apply it to our case. Assuming that failures are evenly distributed a failure occurs at time I 2 on average. The total time lost per MTBF period can be approximated by I 2 + Rs = 179min, that is 12, 4% of the period, giving β = 12.4%. The two other parameters, L and R, are protocol-dependent. In the protocol proposed by Meneses et al. [START_REF] Meneses | Team-based Message Logging: Preliminary Results[END_REF], all inter-cluster messages are logged and only the cluster where the failure occurs roll back. Considering a partition Π = P 1 , P 2 , . . . , P K which entails a volume of B inter-cluster messages over a total volume of messages D,

cost(Π) = 23% × B D + 12.4% × k |P k | 2 P 2 (2)
where

P k |P k | 2 P 2
is the average number of processes restarting after a failure if the failures are evenly distributed among the P application processes. As described in Section 2, the protocol proposed by Guermouche et al. [START_REF] Guermouche | Uncoordinated Checkpointing Without Domino Effect for Send-Deterministic Message Passing Applications[END_REF] only logs half of the inter-cluster messages but requires to roll back K+1 2 clusters on average after a failure:

cost(Π) = 23% × B 2 × D + 12.4% × K + 1 2 × k |P k | 2 P 2 . (3)

Results

We first evaluate the cost of partial message logging. Table 2 presents the results obtained by running our partitioning algorithm with PaToH and cost function [START_REF] Asanovic | The Landscape of Parallel Computing Research: A View from Berkeley[END_REF]. The results show that for all applications except FT and MAESTRO, our tool was able to find a clustering configuration where less than 15% of the processes have to roll back on average after a failure, while logging less than 20% of the data exchanged during the execution.

To get an insight on the quality of the costs obtained, they might be compared to the cost of a coordinated checkpointing and a message logging protocol. With a coordinated checkpointing (L = 0, R = 1), the amount of wasted resources would be 12.4%. With a message logging protocol (L = 1, R = 1 P ), this amount would be 23%. The results show that for all applications except FT, using the partial message logging protocol minimizes the cost. For applications based on dense linear algebra (BT, LU, SP, and PARATEC) or N-body methods (GTC, LAMMPS, and Nbody), the cost obtained is always below 5%.

Table 3 presents the results obtained by running our partitioning algorithm with PaToH and cost function [START_REF] Bailey | The NAS Parallel Benchmarks 2.0[END_REF]. We evaluate it only with the applications having a symmetric communication pattern, because in this protocol, messages logging is not bidirectional and our partitioning tool does not take this in account yet. In all the tests except MAESTRO, we managed to find clusters such that the ratio of rolled back processes is around 55% while logging less than 5.3% of the messages. The cost function of the two partial message logging protocols cannot be compared because the first one is only valid for a single failure case while the second one can handle multiple concurrent failures. Table 3. Partitioning for the protocol of Guermouche et al.

To validate our bisection-based partitioning algorithm, we used PaToH, MeTiS and Scotch as outlined in Section 4.1. Table 4 presents the result of running the three tools on PARATEC with the cost function [START_REF] Asanovic | The Landscape of Parallel Computing Research: A View from Berkeley[END_REF]. We ran the experiment with K = 2, 4, 8, 16, 32, 64, and also with K = 13 which is the result provided by our tool (see table 2). First, it has to be noticed that the number of clusters found by our tool is close to the number of clusters that minimizes the cost function with PaToH and Scotch. Second, only Scotch manages to slightly improve the cost compared to our tool. However, if we use our tool with Scotch instead of PaToH, we obtain exactly the same cost. This is mostly due to the fact that Scotch obtains well balanced partitions for any given K, and that the β in (2) is relatively high. In cases where β is smaller, the partitioner has a higher degree of freedom. Whereas the proposed method automatically exploits this leeway, it is hard to specify the imbalance parameter for the three existing tools we have used (we have not reported these experiments, but this was observed for β around 5). We conclude from these results that the proposed algorithm manages to find good clusters without taking a number of clusters as input. The row "run time" below PaToH contains the running time of PaToH with the given K compared to that of the proposed algorithm (which finds K = 13). As is seen, two runs of PaToH take more time than a single run of the proposed algorithm (despite the overheads associated with repeated calls to the library, including converting the data structure). Partial message logging protocols, combining a checkpointing and a message logging protocol, are an attractive rollback-recovery solution at very large scale because they can provide failure containment by logging only a subset of the application messages during the execution. To work efficiently, such protocols require to form clusters of processes in the application, such that inter-cluster communications are minimized. In this paper, we showed that such clustering can be done in many MPI HPC applications. To do so, we analyzed the communication patterns we gathered from the execution of a representative set of HPC MPI applications. To find clusters, we proposed a bisection-based partitioning algorithm that makes use of a cost function evaluating the efficiency of a partial message logging protocol for a given clustering configuration. Contrary to existing graph partitioning tools, this algorithm does not require the number of clusters as input. We defined a cost function for two state-of-the-art partial message logging protocols and ran tests on our set of execution data. With both protocols, results show that we were able to get a good trade-off between the size of the clusters and the amount of logged messages. Furthermore, with Meneses et al. protocol, percentage of processes to rollback and of message to log is in many cases around 10%, which is an order of magnitude improvement compared to flat rollback-recovery protocols. This result encourages us to continue our work on partial message logging protocols.

1Fig. 1 .Fig. 2 .

 12 Fig. 1. Communication patterns inherent in collective communications of MPICH2

2

 2

Table

	Dense linear	Sparse linear	Spectral	N-body simu-	Structured	Unstructured
	algebra		algebra		methods	lations	grids		grids
	BT,	LU,	CG,	MAE-	FT,	GTC,	MG,	GTC,	MAESTRO
	PARATEC	STRO		PARATEC	LAMMPS,	MAESTRO,
						Nbody	PARATEC

Table 2 .

 2 Partitioning for the protocol of Meneses et al.

		Size	Nb Clusters	Min/Max cluster size	Processes to roll back	Log/Total Amount of data (in GB)	Cost
	NPB LU	1024	16	64/64	53.1%	34/700 (4.8%)	7.7
	MAESTRO 1024	4	250/262	62.5%	27/309 (8.9%)	9.78
	PARATEC 1024	16	61/68	53.1%	1285/23914 (5.3%) 7.8
	LAMMPS	1024	8	127/129	56.2%	0.15/4 (3.8%)	7.9

Table 4 .

 4 Evaluating three tools on PARATEC with the cost function[START_REF] Asanovic | The Landscape of Parallel Computing Research: A View from Berkeley[END_REF].

			4	8	13	16	32	64
	PaToH	7.07	4.52	3.42	3.47	3.36	3.95	5.04
	run time 0.60	0.57	0.60	0.59	0.60	0.61	0.60
	MeTiS	7.47	5.36	5.38	5.47	5.92	3.82	4.96
	Scotch	7.09	4.56	3.44	3.37	3.22	3.67	4.94

Acknowledgments

Experiments presented in this paper were carried out using the Grid'5000 experimental testbed, being developed under the INRIA ALADDIN development action with support from CNRS, RENATER and several Universities as well as other funding bodies (see https://www.grid5000.fr). This work was supported by INRIA-Illinois Joint Laboratory for Petascale Computing and the ANR RES-CUE project.


PhD Fanny Orlhac
email: orlhacf@gmail.com

PhD Frédérique Frouin

PhD Christophe Nioche

PhD Nicholas Ayache

Irene Buvat

PhD Irène Buvat

Validation of a method to compensate multicenter effects affecting CT radiomic features

Keywords: Acrylonitrile Butadiene Styrene CCR, Credence Cartridge Radiomics FBP, Filtered Back Projection Lung reconstruction algorithm

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

Introduction

Since 2012, the concept of Radiomics is expanding [START_REF] Gillies | Radiomics: Images Are More than Pictures, They Are Data[END_REF] in oncology with the objective to characterize tumor heterogeneity from medical images. Radiomics extracts features from medical images that quantify tumor shape, intensity histogram and texture of the lesions more precisely and more accurately than visual assessment by a radiologist, in order to build models involving such features to assist patient management.

In particular, texture analysis from CT images has led to promising results to distinguish between tumor lesions with different histopathological characteristics and to predict treatment response or patients survival [START_REF] Lubner | CT Texture Analysis: Definitions, Applications, Biologic Correlates, and Challenges[END_REF]. However, several studies have highlighted the sensitivity of radiomic features to CT acquisition and reconstruction parameters using phantoms [START_REF] Zhao | Exploring Variability in CT Characterization of Tumors: A Preliminary Phantom Study[END_REF][START_REF] Mackin | Measuring Computed Tomography Scanner Variability of Radiomics Features[END_REF][START_REF] Mackin | Harmonizing the pixel size in retrospective computed tomography radiomics studies[END_REF][START_REF] Shafiq-Ul-Hassan | Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels[END_REF][START_REF] Caramella | Can we trust the calculation of texture indices of CT images? A phantom study[END_REF][START_REF] Berenguer | Radiomics of CT Features May Be Nonreproducible and Redundant: Influence of CT Acquisition Parameters[END_REF] or patient data [START_REF] Kim | Impact of Reconstruction Algorithms on CT Radiomic Features of Pulmonary Tumors: Analysis of Intra-and Inter-Reader Variability and Inter-Reconstruction Algorithm Variability[END_REF][START_REF] Lu | Assessing Agreement between Radiomic Features Computed for Multiple CT Imaging Settings[END_REF][START_REF] Shafiq-Ul-Hassan | Voxel size and gray level normalization of CT radiomic features in lung cancer[END_REF][START_REF] He | Effects of contrast-enhancement, reconstruction slice thickness and convolution kernel on the diagnostic performance of radiomics signature in solitary pulmonary nodule[END_REF]. Indeed, feature values are affected by slice thickness, pixel size, reconstruction kernel, tube voltage, tube current and contrastenhancement. They also differ between different scanners with the same settings [START_REF] Berenguer | Radiomics of CT Features May Be Nonreproducible and Redundant: Influence of CT Acquisition Parameters[END_REF]. Moreover, the impact of imaging protocols varies according to the texture pattern and the radiomic feature [START_REF] Mackin | Measuring Computed Tomography Scanner Variability of Radiomics Features[END_REF].

One of the most widely cited studies in radiomics [START_REF] Aerts | Decoding tumour phenotype by noninvasive imaging using a quantitative radiomics approach[END_REF], which included 1019 patients, used different CT imaging protocols involving different CT scanners, different pixel size and slice thickness, with or without intravenous contrast, without accounting for this variability in the data analysis. To reduce that variability, it has been proposed to resample images with a fixed voxel size, to filter the images [START_REF] Mackin | Harmonizing the pixel size in retrospective computed tomography radiomics studies[END_REF] or to change the definition of features [START_REF] Shafiq-Ul-Hassan | Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels[END_REF][START_REF] Shafiq-Ul-Hassan | Voxel size and gray level normalization of CT radiomic features in lung cancer[END_REF]. These approaches require a modification of the CT images or are not applicable to all radiomic features. The same issue is encountered in PET imaging, where radiomic features are sensitive to the acquisition protocol and reconstruction algorithm [START_REF] Yan | Impact of Image Reconstruction Settings on Texture Features in 18F-FDG PET[END_REF]. A compensation method was initially described in genomics [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF], where the so-called batch effect is the source of variations in measurements caused by the handling of samples by different laboratories, different technicians and on different days. The batch-effect is conceptually similar to variations induced by the scanner or the protocol effects in radiomics. The compensation method identifies a batch-specific transformation to express all data in a common space devoid of batch effects. It has been shown to be effective in PET to realign the radiomic feature distributions between 3 different protocols for healthy liver tissue and breast lesions, without altering the biological information [START_REF] Orlhac | A Postreconstruction Harmonization Method for Multicenter Radiomic Studies in PET[END_REF].

The purpose of this study was therefore to determine whether this compensation method could also correct for the CT protocol effect, using phantom and patient data.

Materials and Methods

All patient data were anonymized and publicly available in supplemental data of ( 9) and [START_REF] Lu | Assessing Agreement between Radiomic Features Computed for Multiple CT Imaging Settings[END_REF]. All authors had control of the data and information submitted for publication.

Phantom experiments

The phantom data used in our study have been produced by Mackin et al [START_REF] Mackin | Measuring Computed Tomography Scanner Variability of Radiomics Features[END_REF] and are publicly available in supplemental data of (4). The Credence Cartridge Radiomics (CCR) phantom consists of 10 layers with different materials corresponding to different texture patterns. This phantom was scanned using 17 different imaging protocols from four medical institutes involving various reconstruction kernels, scan types, slice thickness, pixel spacing, spiral pitch factor and effective milliamperage. Additional information on phantom and acquisition characteristics are provided in Supplemental Tables 1 and2. For each layer, 16 non-overlapping volumes of interest (VOI) with a cubic volume, on average, of 8 cm 3 (range: [7.6 -9 cm 3 ] corresponding to [2708 -14332 voxels] depending on the imaging protocols) are also made available in Dicom-RTstruct format. For each VOI and each imaging protocol, we (FO with 7 years and CN with 20 years of research experience in medical imaging) computed 40 radiomic features using the LIFEx freeware (17) (www.lifexsoft.org, Inserm, Orsay, France, Supplemental Table 3), with a fixed bin size [START_REF] Orlhac | 18F-FDG PET-Derived Textural Indices Reflect Tissue-Specific Uptake Pattern in Non-Small Cell Lung Cancer[END_REF] set to 10 HU between -1 000 UH and 3 000 UH without any spatial resampling. We performed the radiomic analysis for 16 imaging protocols out of 17 due to a reading issue with acquisition CCR1-GE2.

Patients

Publicly available radiomic features from two patient databases (#1 and #2) were used in our study. First set of features was derived from cohort #1 of 42 patients with a lung cancer between September and October 2013 [START_REF] Kim | Impact of Reconstruction Algorithms on CT Radiomic Features of Pulmonary Tumors: Analysis of Intra-and Inter-Reader Variability and Inter-Reconstruction Algorithm Variability[END_REF], including 19 males (mean age, 60.4 years; range, 31-81 years, Table 1). All patients underwent a CT scan with the same machine and protocol (Supplemental Table 4), and CT images were reconstructed using three algorithms: filtered back projection (FBP), Sinogram Affirmed Iterative Reconstruction (Siemens Healthcare, Forchheim, Germany) with a noise reduction strength of level 3 (called S3 thereafter) and 5 (called S5). For each patient, the dominant tumor lesion was segmented manually three times, twice by a radiologist and once by a technologist. For each of the 3 VOIs per patient and each reconstruction, 15 radiomic features were calculated. We (FO) excluded 5 geometric features (volume, diameter, surface, sphericity and compactness) from the analysis as they mostly depend on the segmentation.

Second set of features was obtained from 32 patients of cohort #2 between January and September 2007 (16 males; mean age, 62.1 years; range, 29-82 years, Table 1) with a lung cancer who underwent two CT scans (Supplemental Table 4) within 15 minutes [START_REF] Lu | Assessing Agreement between Radiomic Features Computed for Multiple CT Imaging Settings[END_REF]. This dataset was originally collected in the clinical trial NCT00579852 to evaluate the reproducibility of tumor volume and diameter measurements and is part of the Reference Image Database to Evaluate therapy Response (RIDER) project [START_REF] Armato | The Reference Image Database to Evaluate Response to therapy in lung cancer (RIDER) project: a resource for the development of change-analysis software[END_REF]. The CT images were reconstructed using 6 protocols combining two reconstruction algorithms (lung and standard abbreviated as L and S) and three slice thicknesses (1.25 mm, 2.5 mm and 5 mm) [START_REF] Lu | Assessing Agreement between Radiomic Features Computed for Multiple CT Imaging Settings[END_REF]. For one lesion per patient (29 primary and 3 metastatic lesions), a tumor VOI was obtained from a consensus among the manual segmentations by 3 radiologists. After resampling the VOI voxels to 0.5x0.5x0.5 mm 3 using a trilinear interpolation, 89 radiomic features were calculated for the 6 imaging protocols (2 reconstructions x 3 slice thicknesses) and for each of the 64 scans (32 patients with 2 scans).

Compensation method

To correct for differences in features caused by the various imaging protocols, we (FO) used the ComBat compensation method [START_REF] Johnson | Adjusting batch effects in microarray expression data using empirical Bayes methods[END_REF]. This method has been used for cortical-thickness measurements from MR images [START_REF] Fortin | Harmonization of cortical thickness measurements across scanners and sites[END_REF] and for radiomic features from different PET protocols [START_REF] Orlhac | A Postreconstruction Harmonization Method for Multicenter Radiomic Studies in PET[END_REF]. It is a data-driven method that identifies the protocol effect assuming that the value of each feature y measured in VOI j with imaging protocol i can be written as:

y ij =α+ γ i +δ i ε ij Equation 1
where α is the average value for feature y ij , γ i is an additive protocol effect, and δ i is a multiplicative protocol effect affected by an error term (ε ij ). The compensation consists in estimating the model parameters α, γ i and δ i using a maximum likelihood approach based on the set of available observations y:

y ij ComBat = y ij -α ̂-γ i δi ̂+α ̂ Equation 2
where α ̂, γ i ̂ and δ i ̂ are estimators of α, γ i and δ i .

We used the non-parametric form of the model in which no assumptions are made regarding the laws followed by the parameters. In this setting, ComBat determines a transformation for each feature separately. For each texture pattern of phantom data and of each patient dataset, we used the R (version

Statistical analysis

To determine whether the protocol setting (independent variable i in Equation 1) impacted the distributions of radiomic feature values (dependent variables yij in Equation 1), we (FO, FF with 30 years of experience) performed 2-sided Friedman tests before and after ComBat compensation for each feature as summarized in Supplemental Table 5. Null hypothesis is that there is no difference between the distributions. Benjamini-Hochberg procedure was used to control the false discovery rate [START_REF] Benjamini | Controlling the False Discovery Rate: A Practical and Powerful Approach to Multiple Testing[END_REF]. P-values less than .05 were regarded as statistically significant. As the goal of the compensation is to realign the distributions in terms of mean and standard deviation, a p-value of the Friedman test greater than .05 means that the realignment was successful.

For the phantom data, we also performed a principal component analysis (PCA) of the 2 560 samples (16 VOI x 10 texture patterns x 16 imaging protocols) described by 40 variables (radiomic features). PCA was performed before and after ComBat to visualize the impact of the compensation method on the distinction between patterns. We also studied whether two textural patterns could be distinguished when pooling data from the 3 imaging protocols before and after compensation and for balanced and un-balanced groups.

The statistical analysis was performed using R software.

Results

Patient characteristics are shown in Table 1.

Phantom experiments

In the phantom data, 399/400 p-values of the Friedman tests performed for all features based on 16

imaging protocols and 10 texture patterns were lower than .05 before compensation (Table 2, Supplemental Table 6). Only one p-value for Skewness was larger than .05 for pattern 7 (dense cork, p=.46).

After compensation, all p-values of Friedman tests were higher than .05, demonstrating that the protocol effect was no longer detectable.

These results were confirmed visually using the projection of the data in the space spanned by the first 2 principal components of PCA. Figure 1 shows an overlapped between textural patterns before ComBat, due to the large variability of radiomic feature values computed from 16 different CT protocols. For each textural pattern (each color), several clusters corresponding to different CT protocols could be identified.

After ComBat, textural patterns could be clearly distinguished and were no longer composed of different clusters, demonstrating that the compensation method properly corrected for the scanner effect while retaining the specific characteristics of each texture pattern. Interestingly, the variance explained by the first two components was higher after ComBat (65.6% versus 53.2%), with approximately the same features contributing to the first two principal components before and after compensation (data not shown).

Based on three CT acquisitions (GE1, P2 and S2, see Supplemental Table 2), Figure 2 shows that when data were pooled without realignment, the sensitivity for distinguishing cork from dense cork was 67% (32/48 VOI) with a specificity of 98% (47/48 VOI) using the cutoff maximizing the Youden index. After ComBat, both sensitivity and specificity were 100% (48/48 VOI). For unbalanced groups, Supplemental Figure 1 shows that the compensation method also yielded a perfect distinction between these two patterns.

Patient data

For patient datasets #1 and #2, 100% (10/10) and 98% (87/89) of Friedman tests respectively were lower than .05 between imaging protocols before ComBat (Table 2, Supplemental Tables 78). After ComBat, 30% (3/10) of p-values for dataset #1 and 15% (13/89) of p-values for dataset #2 were lower than .05.

Visual inspection of the radiomic feature value distributions when Friedman tests remained significant after ComBat showed that the residual difference between protocols was always small and that the protocol effect was much reduced (Supplemental Figure 2), demonstrating the effectiveness of the compensation. As illustrated in Supplemental Figure 3 for Homogeneity feature, ComBat corrects the protocol effect with a realignment of feature values among the three protocols for dataset #1 and among the 6 protocols for dataset #2. For instance, before ComBat, the plot shows a shift in distribution with greater Homogeneity values for reconstruction S5 than for S3 and FBP. This was expected since S5 involved higher noise reduction. After compensation, the distributions between the three reconstructions better overlapped. Figure 3 shows three examples of realignment of features between different reconstruction algorithms, reconstruction kernels and slice thicknesses.

Discussion

As widely reported in the literature, radiomic features are sensitive to the acquisition and reconstruction parameters of CT images. Feature values are therefore not directly comparable between different imaging protocols, limiting their use in multicenter studies. Here, we demonstrated that the ComBat method can realign radiomic features computed from different CT imaging protocols. Using phantom data, we showed that ComBat removed the scanner and protocol effect while preserving the differences between texture patterns. The correction for the scanner effect was confirmed using patient images reconstructed with different imaging protocols.

The use of this compensation method should facilitate multicentric radiomic analyses that are absolutely needed to demonstrate the practical usefulness of radiomic features for patient management. Data harmonization is currently a hot topic in the international imaging community with increasing awareness of the need to reduce the variability in image quality between centers and machines [START_REF] Sullivan | Metrology Standards for Quantitative Imaging Biomarkers[END_REF][START_REF] O'connor | Imaging biomarker roadmap for cancer studies[END_REF]. ComBat offers a solution to realign radiomic features with several advantages. ComBat is easily available to all and fast (a function available for free in R software). The transformations are estimated based on the measured feature values, without the need to go back to images or to perform phantom experiments. No learning set is needed. Unlike other options described in the literature, ComBat does not change the feature definitions [START_REF] Shafiq-Ul-Hassan | Intrinsic dependencies of CT radiomic features on voxel size and number of gray levels[END_REF][START_REF] Shafiq-Ul-Hassan | Voxel size and gray level normalization of CT radiomic features in lung cancer[END_REF] and can therefore be used with all software/algorithms and with any radiomic features. This is illustrated in three datasets using three different implementations for the radiomic feature calculation (Supplemental ComBat can model these covariates in the compensation process as illustrated in PET imaging with different proportions of cancer subtypes in different departments [START_REF] Orlhac | A Postreconstruction Harmonization Method for Multicenter Radiomic Studies in PET[END_REF], as long as enough patients with the same characteristics are available.

As demonstrated using the phantom data (Table 2), before compensation, the values of radiomic features were significantly different between imaging protocols for a given pattern. Ignoring this effect in multicentric studies might bias the findings and lower statistical power. A recent study highlighted that the 4 features selected in (13) to build the radiomic signature were highly correlated with tumor volume [START_REF] Vallieres | Dependency of a validated radiomics signature on tumor volume and potential corrections[END_REF], which might explain why the model remained robust on data from different centers. Using ComBat might help determine whether radiomic features reflecting the lesion biological heterogeneity but affected by the center effect more than the lesion volume also have some predictive value.

Our study has some limitations: our findings should be confirmed on other cancer types, for other imaging protocols and scanners and the actual impact on diagnostic performance on clinical data needs to be demonstrated. Independent multi-center validation of radiomic models is also essential for them to become mainstream [START_REF] Gillies | Radiomics: Images Are More than Pictures, They Are Data[END_REF][START_REF] Buvat | Tumor Texture Analysis in PET: Where Do We Stand?[END_REF][START_REF] Reuzé | Radiomics in Nuclear Medicine Applied to Radiation Therapy: Methods, Pitfalls, and Challenges[END_REF].

In summary, ComBat makes it possible to pool radiomic features from different CT protocols. This method appears promising to deal with the center effect in multicenter radiomic studies and to possibly raise the statistical power of such studies. ComBat is data-driven meaning that the transformations identified by

ComBat to set all data in a common space should be estimated for each study involving data from different centers/protocols. Our analysis was based on less than 50 patients for each acquisition protocol demonstrating the efficiency of the method even for small patient cohorts. Using simulations in which we gradually removed patient data (results not shown), we found satisfactory results with as few as 20

patients per imaging protocol. The minimum number of patients required per imaging protocol to successfully apply ComBat remains to be comprehensively investigated. 2). When pooling all radiomic feature values, the optimal cutoff could not perfectly distinguish the patterns while after compensation of scanner effects, a perfect distinction was observed.

Figure 1 :

 1 Figure 1: Phantom data: principal component scores for 2560 samples corresponding to 16 VOI x 10 texture patterns (in colors) x 16 imaging protocols described by 40 radiomic features on the first two principal components before (A) and after ComBat (B). After compensation, each texture pattern was no longer composed of different clusters, demonstrating that the scanner effect has been correctly compensated.

Figure 2 :

 2 Figure 2: Example of ComBat application in phantom experiments: two texture patterns (cork and dense cork) were scanned using three different imaging protocols with 16 volumes of interest in each case (GE1, P2, S2, see Supplemental Table2). When pooling all radiomic feature values, the optimal cutoff could not perfectly distinguish the patterns while after compensation of scanner effects, a perfect distinction was observed.

Figure 3 :Supplemental Figure 1 :Supplemental Figure 2 :Supplemental Figure 3 :Supplemental Table 1 :

 31231 Figure 3: Probability density function(%) of Homogeneity before (left) and after ComBat (right) in patient data using two CT reconstruction algorithms (FBP and S5), two reconstruction kernels (Standard and Lung) and two voxel thicknesses (1.25 mm and 5 mm). CT images are from Kim et al, PLoS ONE 2016;11:e0164924 and Lu et al, PLoS ONE 2016;11:e0166550. Displayed p-values are given for Friedman tests.

  3.4.2, R foundation for Statistical Computing, Vienna, Austria) function called ComBat, available at https://github.com/Jfortin1/ComBatHarmonization to identify the transformation parameters.

Table 5

 5

). ComBat does not require spatial resampling of the CT images to a single pixel size and/or image filtering (5). It is applicable when only radiomic features values are available or when images are not available. ComBat can account for covariates of interest if the patients scanned with different imaging protocols do not have the same characteristics (eg, different age distributions).

Table 1 :

 1 Patient characteristics.

	Patient dataset #1	Patient dataset #2

Table 2 :

 2 number of significant Friedman tests (p<.05) without and with compensation in the phantom and clinical datasets.

	w/o realignment	with realignment

Table 2 :

 2 acquisition characteristics for phantom experiments.

	Designation	Manufacturer/Model	Reconstruc tion kernel	Scan type	Voxel size (mm 3 )	Spiral factor pitch	kVp	mAs
	GE1	GE/Discovery CT750 HD	standard	helical	0.49x0.49x2.5	0.98	120	81
	GE3	GE/Discovery CT750 HD	standard	helical	0.78x0.78x2.5	0.98	120	122
	GE4	GE/Discovery ST	standard	helical	0.98x0.98x2.5	1.35	120	143
	GE5	GE/LightSpeed RT	standard	helical	0.98x0.98x2.5	0.75	120	1102
	GE6	GE/LightSpeed RT16	standard	helical	0.98x0.98x2.5	0.94	120	367
	GE7	GE/LightSpeed VCT	standard	helical	0.74x0.74x2.5	0.98	120	82
	P1	Philips/Brilliance Big Bore	B	helical	0.98x0.98x3	0.94	120	320
	P2	Philips/Brilliance Big Bore	C	helical	0.98x0.98x3	0.94	120	369
	P3	Philips/Brilliance Big Bore	B	helical	1.04x1.04x3	0.81	120	320
	P4	Philips/Brilliance Big Bore	B	helical	1.04x1.04x3	0.81	120	369
	P5	Philips/Brilliance 64	B	helical	0.98x0.98x3	0.67	120	372
	S1	Siemens/Sensation Open	B31s	axial	0.52x0.52x2	1.00	120	26-70
	S2	Siemens/Somatom Definition Flash	170f, 2	helical	0.54x0.54x3	0.60	120	17-28
	T1	Toshiba/Aquilion	FC18	helical	0.63x0.63x3	1.11	120	135
	T2	Toshiba/Aquilion	FC18	helical	0.63x0.63x3	1.11	120	135
	T3	Toshiba/Aquilion One	FC18	helical	0.98x0.98x3	0.99	120

Supplemental


Scale Free Controllability of Large-Scale Networks: an Output Controllability Approach Giacomo Casadei, Carlos Canuda-de-Wit and Sandro Zampieri Abstract-In this paper we consider the problem of controllability and energy consumption for large scale networks. Instead of controlling separately all the nodes of the network we control an output which is defined as some measurement (for instance the average) of the nodes which are not directly controlled. We thus exploit the concept of Output Controllability and the Output Controllability Gramian to analyze the properties of the system. In this context, we show that it is possible to obtain a reduced-order model which makes the Gramian compution and control design much easier. Simulations show that the reduced model is consistent with the original one and for low ratios of controlled nodes, more robust and performing with respect to the original.

I. INTRODUCTION

In recent years, the interest in the topic of networks control have been growing. This interest is motivated by the vast number of applications which nowadays can be framed into the context of networks. Robot swarms ( [START_REF] Olfati-Saber | Flocking for multi-agent dynamic systems: Algorithms and theory[END_REF]), power networks ( [START_REF] Dörfler | Synchronization and transient stability in power networks and non-uniform Kuramoto oscillators[END_REF]) and social networks ( [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: Convergence conjectures and theorems[END_REF]), are just few examples which pushed researchers in studying control problems for networks. Many results about consensus and synchronization can nowadays be found both in the context of linear ( [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF], [START_REF] Wieland | An internal model principle is necessary and sufficient for linear output synchronization[END_REF], [START_REF] Seyboth | On robust synchronization of heterogeneous linear multi-agent systems with static couplings[END_REF], [START_REF] Liu | Robust decentralized output regulation with single or multiple reference signals for uncertain heterogeneous systems[END_REF]) and nonlinear systems ( [START_REF] Arcak | Passivity as a design tool for group coordination[END_REF], [START_REF] Stan | Analysis of interconnected oscillators by dissipativity theory[END_REF], [START_REF] Isidori | Robust Output Synchronization of a Network of Heterogeneous Nonlinear Agents Via Nonlinear Regulation Theory[END_REF], [START_REF] Casadei | About disconnected topologies and synchronization of homogeneous nonlinear agents over switching networks[END_REF], [START_REF] Casadei | Multi-Pattern Output Consensus in Networks of Heterogeneous Nonlinear Agents with Uncertain Leader: a Nonlinear Regression Approach[END_REF]). In recent years, due to the increasing complexity of networks applications, the focus has shifted towards the framework of large scale networks. The analysis of large scale networks has been widely covered in the literature and notable contributions can be found in [START_REF] Newman | The structure and dynamics of networks[END_REF], [START_REF] Newman | The structure and function of complex networks[END_REF], [START_REF] Albert | Statistical mechanics of complex networks[END_REF]. In particular, in [START_REF] Barabási | Emergence of scaling in random networks[END_REF], [START_REF] Barabási | Scale-free networks[END_REF], it has been observed that many real world large networks such as the world wide web present a scaling effect in their growth, namely the connectivity degree follows a power law distribution.

The problem of controlling large scale networks is with many extents still an open problems. A fundamental aspect that has to be considered is the energy necessary to control a system and with this respect the analysis of large networks have been addressed by several authors ( [START_REF] Liu | Controllability of complex networks[END_REF], [START_REF] Cowan | Nodal dynamics, not degree distributions, determine the structural controllability of complex networks[END_REF], [START_REF] Yan | Controlling complex networks: How much energy is needed?[END_REF]). It has been shown that several aspects such as network centrality [START_REF] Bof | On the Role of Network Centrality in the Controllability of Complex Networks[END_REF], [START_REF] Lindmark | A driver node selection strategy for minimizing the control energy in complex networks[END_REF], exact controllability [START_REF] Yuan | Exact controllability of complex networks[END_REF], [START_REF] Lindmark | Minimum energy control for networks of coupled harmonic oscillators[END_REF] and number of control nodes [START_REF] Yan | Spectrum of Controlling and Observing Complex Networks[END_REF] plays a fundamental role in G. Casadei and C. Canudas-de-Wit are with Univ. Grenoble Alpes, CNRS, Inria, GIPSA-lab, F-38000 Grenoble, France (email:

giacomo.casadei@gipsa-lab.fr, carlos.canudas-de-wit@gipsalab.fr). S. Zampieri is with Universitá di Padova, (e-mail: sandro.zampieri@dei.unipd.it): The research leading to these results has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation program (grant agreement N • 694209). assessing controllability. It is however important to distinguish between the classic controllability concept and practical controllability. Classic controllability is defined as a property of the state and input maps while practical controllability is related to the energy necessary to control the system. With this respect, in [START_REF] Pasqualetti | Controllability Metrics, Limitations and Algorithms for Complex Networks[END_REF] the authors have shown that an anlysis of the practical controllability of a large network is really hard to perform.

This work is carried out within the context of the ERC, Scale-FreeBack (see [START_REF]Scale-Free Control for Complex Physical Network systems[END_REF]) which aims to develop holistic scale-free modeling, estimation and control methods for complex network systems. Scale-free control refers to the ability of a control algorithm to easily scale to any network dimension. Several large-scale systems are described by homogeneous networks in which nodes have practically the same degree and are inherently difficult to model and control. These networks can however be abstracted/aggregated/reduced to a Scale-Free graph characterized by heavy-tailed node/link distributions obeying a power law, for which the underlying network structure is dominated by relatively small numbers of nodes (hubs) having a large number of connections. It is expected that controlling scale-free networks by acting only on a limited number of lumped nodes, the control design for large-scale systems will be tractable. This idea boils dawn to control some few outputs (for instance representing the average value of a large group of states), rather than trying to control the network full state. In this context, we can adapt the conventional controllability to the concept of output controllability and consequently the Controllability Gramian to the Output Controllability Gramian. The latter allows to evaluate the controllability and the energy necessary to steer the output to a desired value. By its definition however, the Output Controllability Gramian still requires the computation of the Controllability Gramian and thus presents the problem mentioned previously.

To overcome this issue, in this paper we propose a outputbased model reduction to evaluate the controllability of a large scale network. In particular, we define a reduced model which aggregates the dynamics which contributes to the output. The reduced-model simplifies the computation of the Controllability Gramian and the minimum energy control inputs to steer the ouput of the system to a desired value. Simulation results shows that this reduced system is more robust to noise and ill-conditionement of the Gramian and, for small ratios of controlled nodes, control design based on the reduced model performs better than the design based on the original system. Furthermore, thanks to the simplicity of the reduce system Gramian, a relationship between the control energy, the networks control-degree and inertia can be displayed.

II. PROBLEM FORMULATION Consider a system

ẋ = Lx + Bu (1)
with x ∈ R n , u ∈ R m , L ∈ R n×n a stable matrix and B ∈ R n×m . As far as controllability is concerned, we recall that the Controllability Gramian is defined as

W t = t τ =0 e L(t-τ ) BB T e L T (t-τ ) dτ . (2)
Given a certain A, we could define a family of W which varies according to m. It can be proved that (1) could be generically controllable for any m ≥ 1: on the other hand, it could be practically uncontrollable from an energetic point of view depending on the particular values of m. The energy necessary to steer the system from x(0) = 0 to a desired final value x(t) = x f is defined as

E(x f , T ) = t 0 u o (τ ) 2 dτ = x T f W -1 t x f (3)
where u o (•) is the optimal input which steers the system from the initial state x 0 = 0 to x(t) = x f while minimizing the control energy. As a matter of fact, it is possible to define different controllability metrics which derives from the Controllability Gramian W t :

• worst case energy: simple calculation shows that max

x f =1 E(x f , t) = λ min (W -1 t )
namely the the larger λ min (W K,t ) the more controllable the system; • average control energy: it can be shown that

Tr(W -1 t )
measures the average control energy: in particular the smaller Tr(W -1 t ) the more controllable the system; • quality of complete controllability: the determinant of W t measures the volume of the hyperellipsoid of final states x f reachable from the initial state x 0 = 0 and with control input energy E(x f , t) ≤ 1.

The analysis of the Controllability Gramian plays a fundamental role in assessing if a system is practically controllable or not. However, for large scale systems the computation of the Gramian (2) is computationally costly and typically suffers of noise and ill-conditionement. The computation of the metrics introduced before is thus a difficult task, subject to huge computational inaccuracy. Furthermore, if one seeks to solve an optimal control problem to steer a system to a desired final state x f , solutions which depends on the Controllability Gramian would be really unreliable.

A. Output controllability

In many engineering applications we are interested in controlling an output of the system rather than its full state. A system

Σ x : ß ẋ = Lx + Bu y = Cx (4)
with y ∈ R r , is said to be output controllable if for any y f ∈ R r and t > 0 there exists an input function u(τ ), 0 < τ ≤ t, such that the output goes from y(0) = 0 to y(t) = y f . This property can be checked through the rank of the output controllability matrix.

Theorem 1 System (4) is output controllable if and only if the matrix C O ∈ R r×nm C O = [CB CLB CL 2 B . . . CL n-1 B]
has full row rank.

In general, assuming output controllability is less restrictive than assuming (standard) controllability. In terms of energy, it seems plausible that controlling an output in place of each single state of the system would require less energy. As a matter of fact, the energy necessary to steer the output y(0) = 0 to a desired final value y(t) = y f can be expressed as a function of the Output Controllability Gramian

W o = C W C T . (5)
where for the sake of simplicity in the notation we dropped the subscript K, t.

Lemma 1 Given a system in the form of (4), the minimal energy necessary to steer the output y(0) = 0 to a desired final value y(t) = y f is given by

E(y f , t) = y T f W -1 o y f (6)
where W o is (5).

Proof: The minimal energy necessary to steer the output to the final value can be obtained by solving the optimal problem

min u(•) s.t. y(t) = y f t 0 u(τ ) T u(τ ) dτ (7)
We have

H = u T u + λ T (Lx + Bu) (8)
H x = λ T L (9)
H u = 2u T + λ T B ( 10
)
with constraint

X (x f ) : R n → R m , X (x f ) = Cx (11)
We seek a solution to the two-boundary value problem

ẋo = Lx o + Bu o x o 0 = x 0 (12) λo = -H T x λ o f = X T x ν (13)
Namely we have

λ = -λ T L ⇒ λ(t) = e -L T t λ 0 λ f = e -L T t f λ 0 = C T ν
Thus we have λ 0 = e L T t f C T ν and as a consequence

λ(t) = e L T (t f -t) C T ν (14)
and as a consequence

H u = 0 ⇒ u = - 1 2 B T λ = - 1 2 B T e L T (t f -t) C T ν (15)
In order to find ν, we have

x(t) = e Lt x 0 - 1 2 ñ t 0 e A(t-τ ) BB T e A T (t-τ ) dτ ô C T ν
which for t = t f gives us r constraints

y f = Ce Lt f x 0 - 1 2 C ñ t f 0 e L(t f -τ ) BB T e L T (t f -τ ) dτ ô C T ν we obtain ν = 2 W o Ce At f x 0 -Cx f (16)
Finally, we have

u o (t) = -B T e A T (t f -t) C T W -1 o Ce At f x 0 -Cx f (17)
For the sake of simplicity we set x 0 = 0 and by introducing ( 17) into ( 7) we obtain [START_REF] Wieland | An internal model principle is necessary and sufficient for linear output synchronization[END_REF].

The analysis of the output Controllability Gramian defined in [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF] allows to understand the input-output energy relationship for system [START_REF] Mirtabatabaei | Opinion dynamics in heterogeneous networks: Convergence conjectures and theorems[END_REF].

At this stage, the Output Controllability Gramian W o still requires the computation of the Controllability Gramian W. In case of large scale networks, the problem introduce in previous sections are still present. In order to overcome them, in Section III we introduced a reduced model which allows to simplify the calculation of the Controllability Gramian and thus the Output Controllability Gramian and the minimal energy control input.

III. A REDUCED MODEL NETWORK

In this section we consider a generic network system in the form of [START_REF]Scale-Free Control for Complex Physical Network systems[END_REF]. We suppose that a certain number m of nodes can be directly controlled. Thus system (4) can be expressed as ïz

δ ò = ï L 11 L 12 L 21 L 22 ò ï z δ ò + ï I m 0 ò u y = 0 1×m H T ï z δ ò (18) with L Metzler, z = [x 1 , . . . , x m ] T , δ = [x m+1 , . . . , x n ]
T and H T is a full vector, namely

H T = h 1 . . . h n-m , h i = 0 ∀ i .
The partition of L, B and C follows the partition of the system state into [z, δ] T . In particular, z represents the nodes which can directly be controlled through the control input u while δ represents the uncontrolled nodes in the network. The output is then defined as a measurement of the uncontrolled nodes (for instance their average).

With the goal of controlling the output in mind, in order to simplify the computation of the Output Controllability Gramian and minimal energy control input, we aim to define a reduced model for [START_REF] Barabási | Scale-free networks[END_REF]. We consider the possibility to aggregate the dynamics of the uncontrolled nodes and to express the evolution of the output as a reduced dynamical system.

To this end, we define a matrix P ∈ R (m+1)×n as

P = ï I m 0 0 H T ò . ( 19
)
Accordingly we introduce a new representation for (1)

x = P x .
As a matter of fact, by considering [START_REF] Barabási | Scale-free networks[END_REF] we can write

ẋ = ï L 11 L 12 H T L 21 H T L 22 ò x + ï I m 0 ò u = ï L 11 H T L 21 ò z + ï L 12 H T L 22 ò δ + ï I m 0 ò u (20)
Furthermore, we have that

P + = ï I m 0 0 1 σ H ò , P + P = ï I m 0 0 1 σ HH T ò
where σ = H 2 and P + ∈ R n×(m+1) is such that P P + = I m+1 . As a consequence we can write

x = P + x + (I n -P + P )x = ï I m 0 0 1 σ H ò x + ï 0 0 0 I η -1 σ HH T ò ï z δ ò
where η = n -m and which through simple computations leads to

δ = 0 1 σ H x + (I η - 1 σ HH T )δ
Thus (20) can be rewritten as

ẋ = ñ L 11 1 σ L 12 H H T L 21 1 σ H T L 22 H ô x + ï I m 0 ò u + ñ L 12 (I η -1 σ HH T ) H T L 22 (I η -1 σ HH T ) ô δ (21)
System ( 21) is a reduced model for (1), namely x ∈ R n and x ∈ R m+1 . In case y represents the average of the aggregated nodes, the third term in [START_REF] Yan | Controlling complex networks: How much energy is needed?[END_REF] can be seen as the variance of the measure of the states [x m+1 , . . . , x n ] with respect to y.

To simplify the computations in the following section, we set the control input u as

u = 0 1 σ L 12 H x + υ (22)
and rewrite [START_REF] Yan | Controlling complex networks: How much energy is needed?[END_REF] as

ẋ = ï L 11 0 H T L 21 1 σ H T L 22 H ò x + ï I m 0 ò υ + ñ L 12 (I η -1 σ HH T ) 1 σ H T L 22 (I η -1 σ HH T ) ô δ (23)
In compact form we rewrite [START_REF] Lindmark | A driver node selection strategy for minimizing the control energy in complex networks[END_REF] as

ẋ = Lx + Bυ + Dδ y = Ĉ x (24)
where Ĉ = 0 1×m 1 . It is readily seen that the system is output controllable as long as H T L 21 = 0. Thus, by using Lemma 1, the following holds.

Lemma 2 Given a system in the form of [START_REF] Yuan | Exact controllability of complex networks[END_REF], the minimal energy necessary to steer the output y(0) = 0 to a desired final value y(t) = y f is given by

Ê(y f , t) = y 2 f Ŵ-1 o ( 25
)
where Ŵo is Controllability Gramian of system (24).

Lemma 2 holds for any system [START_REF] Yuan | Exact controllability of complex networks[END_REF]. In particular if D = 0 m+1×η , the energy E(y f , t f ) for [START_REF] Barabási | Scale-free networks[END_REF] and Ê(y f , t) for ( 24) coincides. However the case of D = 0 m+1×η , commonly referred to as exact reduction, holds only for particular structure of the system and particular choices of H (see [START_REF] Aoki | Control of Large-Scale Dynamic Systems by Aggregation[END_REF] for more details on the subject).

We are thus interested in comparing the Output Controllability Gramian W o for the original system [START_REF] Barabási | Scale-free networks[END_REF] with the Ŵo for [START_REF] Yuan | Exact controllability of complex networks[END_REF]. In the next Section we will show that Ŵo displays some interesting properties and simulation results confirms that its properties are closely related to W o .

A. The Output Controllabilty Gramian of the reduced system

We recall that for a matrix  ∈ R (m+1)×(m+1) with structure

 = ï A 0 a T α ò ( 26
)
where A ∈ R m×m , a ∈ R m×1 and α ∈ R, simple computations show that

Ân =    A n 0 a T n-1 i=0 α n-i-1 A i α n   
As a consequence, by using the Taylor expansion

e Ât = ∞ i=0 Âi t i i!
the following holds.

Lemma 3 The exponential of (26) preserves a block triangular structure and in particular can be expressed as

e Ât = ï e At 0 Λ e αt ò ( 27
)
where

Λ = ∞ i=0 a T (αI m -A) -1 (α i I m -A i ) i! t i .
Proof: We define A 1 = blkdiag(A, α) and

A 2 = ï 0 0 a T 0 ò where A 2 is nilpotent. We have that e Ât = ∞ i=0 (A 1 + A 2 ) i t i i!
Thus by considering (A 1 + A 2 ) i we have

(A 1 + A 2 ) 2 = A 2 1 + ï 0 0 αa T ò + ï 0 0 a T A 0 ò (A 1 + A 2 ) 3 = A 3 1 + ï 0 0 α 2 a T 0 ò + ï 0 0 αa T A 0 ò + ï 0 0 a T A 2 0 ò . . . = . . . (A 1 + A 2 ) i = A i 1 + a T (αI m -A) -1 (α i I m -A i )
Then, by definition of A 1 and A 2 , (27) follows straightforward.

For the sake of compactness we rewrite L as

L = ï L 0 T ϕ ò
and with a slight abuse of notation

e Lt = ï e Lt 0 Λ e ϕt ò .
By considering the definition given in ( 5), the Output Controllability Gramian of ( 24) is

Ŵo = Ĉ Ŵ ĈT = Ĉ ñ t 0 e L(t-τ ) B BT e LT (t-τ ) dτ ô ĈT (28)
In general, the (full) Controllability Gramian Ŵ of (24) will have a structure of the form As a consequence, the Output Controllability Gramian reads as

Ŵ = ñ Ŵ11 Ŵ12 Ŵ21 Ŵ22 ô which,
Ŵo = T L ñ e 2 τ Im 2ϕ + (2L) -1 e 2Lτ
-2(ϕI m + L) -1 e (ϕIm+L)τ t τ =0 LT (30) Since L is stable Metzler matrix, ϕ < 0 and L is Hurwitz: thus the first and third terms in the square brackets are vanishing. As a consequence, the second term is vanishing to. Then for t → ∞ the Output Controllability Gramian converges to

Ŵo = -T L ï 1 2φ I m + (2L -1 ) -2(ϕI m + L) -1 ) ò L = -T L (ϕI m -L) 2 2ϕL(ϕI m + L) LT = T Φ (31) where Φ is defined as Φ = - 1 2ϕ (ϕL + L) -1
As far as T and Φ are concerned, we can observe that: i) T express the relationship between the controlled nodes and the output. In networks, this boils down to saying that

T ∝ d z→δ (32)
where d z→δ is the out-degree from controlled nodes z to aggregated nodes δ; ii) Φ can be seen as the network inertia of controlled nodes and the aggregated nodes: it can be bounded as

Φ ≤ 1 -2ϕ 2 λ min (L) -2ϕλ min (L) 2 (33)
where λ min (L) is the smallest eigenvalue of L. From i), we can infer that controllability depends on the number of connections between controlled nodes and the aggregated nodes: the more connections, the less energy required to control the output of the aggregated nodes. From ii), we can conclude that the smaller the time-constant of the aggregated nodes ϕ and the minimum eigenvalue of the controlled nodes λ min (L), the bigger Φ and thus the less the energy required to control the output of the aggregated nodes.

IV. APPLICATION TO AVERAGE-CONTROL OF LARGE SCALE FLOW NETWORKS

In this section we consider the application of the concept presented in this paper to the case of Large Scale flow networks. We consider networks of n nodes with a ratio of controlled nodes K = m n that varies between K = { 1 2 , . . . , 1 40 }. As long as the ratio of controlled nodes is K ≥ 1 5 , full controllability attains the desired performance, however the energy required is much higher (10 5 order) than the case of output control. Once K < 1 5 , the illconditioned Controllability Gramian makes the system diverge from the desired output and the energy grow even more.

The controlled nodes are selected for each network size and ratio as the one with lower betweenness centrality, namely the boundaries of the network. The network dynamics can be written as

Σ x : ß ẋ = Lx + Bu y = Cx ( 34
)
where L is a Hurwitz Metlzer matrix: this choice is consistent with transportation networks where the dynamics of the nodes can be described by a conservation law (Laplacian matrix) with some additional dissipative terms on the boundaries (sinks of the network) [START_REF] Rantzer | Scalable control of positive systems[END_REF]. Following the transportation networks application, the graph has a pseudo-normal distribution centered around 4, representing a Manhattan grid. The Laplacian matrix of the network is normalized and thus |λ i (L)| ≤ 2 [START_REF] Chung | Spectral Graph Theory[END_REF]. By sorting the agents in such a way that the controlled nodes of the graph are separated from the others, we can rewrite Σ x as [START_REF] Barabási | Scale-free networks[END_REF] with

B = ï I m 0 ò C = 0 1×m H T (35)
In particular, we set H T = 1 η 1 T , namely the average of the uncontrolled nodes. This choice is motivated again by flow networks where typically we are allowed to control the boundaries of the graph with the goal to steer some measure of the interior nodes (in this case the average) to a desired value. In case of traffic networks for instance, the controlled nodes would be the gates to the networks and the output the average density of the interior nodes.

A. Full Controllability and Output Controllability

In order to justify the premise to this paper, we first show that in case of large scale networks Output Controllability is computationally much more convenient than Full Controllability. This advantage appears evident already for small networks. We consider a network of n = 100 nodes and we compare the case of Full Controllability with the Output Controllability: we set a target value y f = 1 and a t f = 10 sec and seek for a solution to [START_REF] Seyboth | On robust synchronization of heterogeneous linear multi-agent systems with static couplings[END_REF]. The Full Controllability solution for system (34) is given by

u o (t) = -B T e L T (t f -t) W -1 e Lt f x 0 -1 η y f (36)
while the Output Controllability solution is given by

u o out (t) = -B T e L T (t f -t) C T W -1 o Ce Lt f x 0 -y f . (37)
We define a performance measure as the ratio between the output obtained by setting the two control inputs just defined and the desired output y f .

In Figure 1, it can be seen that as long as K ≥ 1 5 , the full controllability attains the same desired output of the output controllability however with a much higher energy. Once K < 1 5 , the Controllability Gramian of the full system is ill-conditioned to the point of making the system unable to reach the desired value. The energy as well increases rapidly up to an unrealistic value. On the other hand, Output Controllability still reaches the desired value and the energy increases proportionally to the decrements of control nodes.

B. Boundary Control of Normally Distributed Graphs

We are now interested in considering the properties of the reduced system and to evaluate the performance obtained by designing control based on this approximation. In the case presented in this section, the reduced system (23) reads as Σx : By considering the reduced system Σx , the computation of the Controllability Gramian Ŵ becomes much easier and in turn the computation of the minimal energy control input ûo . The solution to [START_REF] Seyboth | On robust synchronization of heterogeneous linear multi-agent systems with static couplings[END_REF] for Σx is given by

               ẋ = ï L 11 L 12 1 1 η 1 T L 21 1 η 1 T L 22 1 ò x + ï I m 0 ò u + ñ L 12 (I η -1 η 11 T ) 1 η 1 T L 22 (I η -1 η 11 T ) ô δ y = 0 1×m 1 x (
ûo out (t) = -BT e LT (t f -t) ĈT Ŵ-1 o î Ĉe Lt f x0 -y f ó . (39)
We are interested in comparing the behavior of the system Σ x in the case in which the control input u is designed according to (37) and in the case in which it is designed based on the reduce system, namely (39). We again set y f = 1 and t f = 10 sec. For the sake of simplicity, we impose x 0 = 0 and consequently x0 = 0.

As shown in Figure 3, the specturm of L overlaps the spectrum of L fro any K: clearly, due to the reduction of size from n to m + 1, the spectrum of L has less eigenvalues but they cover the same interval. This is a crucial aspect since the input design highly depends on the spectrum of L and L: with this respect it it worth pointing out that considering the normalized Laplacian matrix plays an important role in guaranteeing that the original spectrum and the reduced spectrum are covering the same range (a condition that otherwise would hold only for regular graphs). As a consequence, we can expect that (37) and (39) will have a similar behavior.

In order to compare the performance of the system with respect to the two inputs, we consider the ratio between the output of (34) when (37) is used and when (39) is used with the reference value y f = 1. Figure 4 shows the perfromance of these two cases for different values of control nodes m. As long as the K ≥ 1 5 the control input (37) (red plot) outperforms (39) (black plot) of around 20%. However, we find out that around the ratio K = 1 6 a phase-transition appears in the controllability properties of a network.

In particular, when considering input (37) designed over the original system, we observe a critical change in the performance of the system. On the other hand, this phenomena does not appear when we design the input as (39), namely on the reduced system. The performance of the system with (37) degrade abruptly and greatly, when K ≤ 1 6 while the behavior of the network with (39) is much more regular, namely the performance decreases proportionally to a decreasing number of control nodes.

The origin of this behavior seems to lie in the eigenvalue of the Controllability Gramian. The Controllability Gramian computed over (34) is fragile with respect to K and suddenly becomes inaccurate and highly ill-conditioned. On the other hand, the reduced system (38) is more robust and the control input (39) preserves acceptable performance for smaller K. As reported in Figure 5, this behavior is scale-free with respect to the size of the network. For input (37), we can still observe the phase-transition appearing around K = 1 6 and an abrupt degradation of the performances for smaller number of control nodes. With input (39), the system display a much more robust behavior for smaller K.

In a bigger picture, we can say that the design (39) is also scale-free also with respect to K, namely the degradation of performance caused by a reducing number of controlled nodes is almost neglectable.

C. The out-degree of controlled nodes, the Inertia and the Control Energy

The the performance observed in previous simulations is consistent also with the evolutions of the Output Controllability Gramian W o of the original system and Ŵo of the reduce system. For a network of n = 500 nodes, Figure 6 shows the ratio between between Ŵo and W o . For any K, underestimates W o around 15%-20%. Regarding the Output Controllability Gramian, at the end of Section III-A we showed that Ŵo is closely related to the out-degree of the controlled nodes z towards the aggregated nodes δ (32) and to the network inertia (33). These two aspects, which can be easily inferred through the analysis of the reduced order system, are indeed confirmed by simulation results.

In Figure 7, we considered again a network of n = 500 nodes with a fixed number of controllers m = 100. We plot the performance of the system and the control energy with the two inputs (37) and (39) respectively, for increasing connectivity degrees d z→δ . It can be seen that the performance with input (39) improve for higher d z→δ . Most importantly, the control energy decreases proportionally to an increasing d z→δ . Similarly, in Figure 8, we plot performance and energy of the network with respect to Φ in (33). We can clearly see that for an increasing Φ , the energy required to achieve the desired goal decreases greatly.

These simulations allows to say that the reduced system (38) is a good approximation of the original system (34) from an energetic perspective. As a matter of fact, the relationship between the control energy and the network structure (obtained in (32) and (33)) is shown to have a real impact.

V. CONCLUSION

In this paper we considered the problem of controlling a large scale network. Rather than controlling the full state, we consider the case in which we are interested in controlling an output which depends on the non-controlled nodes. Thanks to this choice, the network is practically controllable independently from the size.

In order to simplify the computation of the Controllability Gramian and minimum energy control inputs, we considered a model-reduced approximation of the system, base on the aggregation of the non-controlled nodes. We showed that this new systems is reliable both for estimating the energy necessary to control the output and also to design optimal control. The reduction of complexity due to the aggregation simplifies a lot the computation and ultimately shows to be more reliable than performing computation over the large original system.

0 e 0 e

 00 given definition of Ĉ = 0 1 , leads to the Output Controllability GramianŴo = Ŵ22As a matter of fact, considering the structure of B we have I m -e Lτ )(e ϕτ I m -e L T τ )dτ LT[START_REF] Rantzer | Scalable control of positive systems[END_REF] where we have defined L = (ϕI m -L) -1 . The integral term reads as t 0 e 2ϕτ I m -e ϕτ Im e Lτ -e ϕτ Im e L T τ + e Lτ e L T τ dτ = = t 2ϕτ I m -e (ϕIm+L)τ -e (ϕIm+L T )τ + e Lτ e L T τ dτ which, in case L = L T (the case of undirected graphs), becomes t 2ϕτ Im + e 2Lτ -2e (ϕIm+L)τ dτ .

Fig. 1 :

 1 Fig.1: Full Controllability and Output Controllability. As long as the ratio of controlled nodes is K ≥1 5 , full controllability attains the desired performance, however the energy required is much higher (10 5 order) than the case of output control. Once K <1 5 , the illconditioned Controllability Gramian makes the system diverge from the desired output and the energy grow even more.

Fig. 2 :

 2 Fig. 2: Connection Distribution for the original system and of the reduced order system: the green histograms shows the original network, namely a network with pseudo-Normal Distribution for n = 500. The blue histograms represents the distribution of the network obtained after the model reduction is preformed for K = { 1 10 , 1 20 , 1 40 }.

Fig. 3 :

 3 Fig. 3: Eigenvalue Distribution of L and L: for a network of n = 500 nodes and K = { 1 10 , 1 20 , 1 40 }. The spectrum of L covers the same spectrum of L for any K.

  38)For a graph of n = 500 nodes, the connection distribution before and after the reduction is shown in Figure2: as already mentioned the original network as a Manhattan structure, with degree centered around 4 connections. The distribution for (38) depends on the number of controlled nodes and here we display the case of K = { 1 10 , 1 20 , 1 40 }. We can clearly see that the aggregated nodes become a unique node which possesses many connections towards the boundary nodes.

Fig. 4 :

 4 Fig. 4: Output of the system with two different control inputs u o and ûo , respectively calculated on the original system and on the reduced system for different size of the network. Control nodes are sitting on the boundary of the network.

Fig. 5 :Fig. 6 :

 56 Fig. 6: Ratio between the Outptu Controllability Gramian of Σx and Σx calculated over 100 different graphs.

Fig. 7 :

 7 Fig. 7: Performance and Energy with respect to the degree of controlled nodes towards aggregated nodes: performances both with u o out and ûo out increases. Most importantly, the control energy decreases with an increasing d z→δ .

Fig. 8 :

 8 Fig. 8: Performance and Energy with respect to inertia of the network Φ : performances both with u o out and ûo out increases. Most importantly, the control energy decreases for an increasing Φ .