M Colonna

Ph G F Burgio 1 '-

M Di Chomaz^

J Toro~

Randrup

M Colonna 1 '

G F Burgio

-'

Ph Chomaz

M Di Toro'

J Randrup

Simulating the Langevin force by simple noise in nuclear one-body dynamics *

For the purpose of addressing catastrophic phenomena in nuclear dynamics, we explore the possibility of simulating the stochastic part of the collision integral in the Boltzmann-Langevin model by the numerical noise associated with the finite number of test particles in the ordinary BUU treatment. Considering idealized twodimensional matter, for which it is practical to simulate the Boltzmann-Langevin equation directly, we demonstrate that the number of test-particles per nucleon can be adjusted so that the corresponding BUU calculation yields a good reproduction of the spontaneous clusterization occurring inside the spinodal region. This approximate method may therefore provide a relatively easy way to introduce meaningful fluctuations in simulations of unstable nuclear dynamics.

Introduction

Intense experimental and theoretical efforts are currently being directed at understanding the mechanisms responsible for the observed production of complex fragments in heavy-ion collisions. In recent years, intermediate-energy heavy-ion collisions have been studied extensively with BUV-type models, in which the collisionless Vlasov equation for the reduced one-body phase-space density /(r,p) is augmented by a Pauli-blocked Boltzmann collision term of the Uehling-Uhlenbeck form [1]. The resulting BUU equation describes the average one-body density and is suited for the description of one-body observables, such as inclusive particle spectra. Models of this class have been quite successful in describing the structure of the equilibrated primary sources formed during the first stage of strongly damped and (incomplete) fusion reactions [2].

However this mean-trajectory approach cannot provide a description of phenomena exhibiting large fluctuations or instabilities, such as multifragmentation processes. In such cases, a fully dynamical picture is essential for predicting particular space geometries in fragment production which should be ultimately related to fundamental properties of nuclear interactions in the nuclear medium. Thus there is a clear need for extending the usual mean-trajectory transport models so as to include fluctuations.

Generally speaking, in order to describe the formation of many fragments one should have to solve the full many-body problem, which is impractical. However, alternative avenues considering ensembles of one-body densities have been investigated [3,4,5,6]. These approaches are equivalent to the so-called Boltzmann-Langevin model (BL) in which the stochastic part of the two-body interaction is added to the mean BUU evolution. In particular, using an exact lattice simulation of the BL equation, it was shown [7] that in presence of instabilities, the fluctuations corning from the stochastic part of the collision integral are propagated and amplified by the mean field leading finally to the formation of clusters. Up to now this approach has been applied only to an idealized two-dimensional system [7], and it appears to be unfeasible to treat a three-dimensional scenario, such as a nuclear collision, because of the large amount of computational effort required.

We have therefore investigated the possibility of replacing the (usually complicated) physical fluctuations by a simpler noise. The main idea is that the dynamics of very unstable systems may be less sensitive to the particular nature of the fluctuations, because the most unstable modes will quickly dominate. Specifically, we study the possibility of replacing the physical fluctuations arising from the collision integral by the numerical noise present when the ordinary BUU model is solved by the test-particle method, as is most commonly done. From a practical point of view, such a simplified approach would allow us to perform approximative calculations for realistic three-dimensional processes.

v \ t*
if.

The Boltzmann-Langevin model

The Boltzmann-Langevin equation can be written on compact form as follows,

(1)

UL

where f(r,p,t) is the one-body phase-space density. On the left-hand side, #[/](r,p) is the self-consistent effective single-particle Hamiltonian governing the collisionless Vlasov evolution. One the right-hand side, /[/] denotes the effect of the residual two-body collisions, which can be decomposed into an average effect, /[/], and the fluctuating part SI[f]. The standard B U U treatment retains only the average term /[/] and therefore produces a single dynamical trajectory for /(r,p). The fluctuating collision term 6I[f] acts like a random "kick"on the one-body density and may therefore produce branchings of the dynamics, resulting in a bundle of different histories.

The physical origin of the fluctuations can be understood as follows [8]: the actual number Af of elementary collisions has a Poisson distribution characterized by its mean value N and so has the variance CTY = N. The average part of the collision integral, /[/], expresses the evolution corresponding to putting N = N for all such elementary processes, while the fluctuating part accounts for the stochastic remainder SN = N -N.

In order to gain insight into the growth of fluctuations we start from linear response theory, which is valid for small fluctuations and therefore is well suited for describing the early evolution of the system. Let us denote by v a particular unstable mode of the BUUdynamics and let its associated imaginary energy be E u = ih/t u . Looking at the linear response to the Boltzmann-Langevin evolution, it is possible to demonstrate [9] that early on the average density fluctuation (rf, in the unstable mode v satisfies an equation of motion of the form,

d 37 at 2 -0V t., (2)
where D v is the source term arising from the stochastic nature of the collision integral. It can be expressed as follows [9],

drdp r dr dp where a(r,p;r',p') = d/dt -< Sf(r,p)8f(r',p') v, =u is the early growth rate of the correlation function and F^(r,p) is the eigenmode associated with the unstable mode.

It follows that the fluctuations of the mode v exhibit the following behaviour in time, *l(t) * D u t v (e 2 '/'--1) + cr 2 (0) e 21 /'-, (4) where 0^(0) is the initial fluctuation. This result shows that the terms D u t v and 0^(0) play similar roles in the evolution. Therefore, it may be possible to replace the physical source term D v by a suitable initial noise 0^(0) determined by the relation

^(0) = ZU, . (5)
This is the idea what we will now test in a simplified scenario.

t v \
3 Test in an idealized system

In order to make a first quantitative test of the idea descibed above, we have considered a gas of Fermions situated on a two-dimensional torus with side lengths equal to I,* = 63 fm and L y = 21 fm. For the effective one-body field we employ a simplified Skyrme interaction, pa pa

with A --100.3 MeV and B = 48 MeV. Moreover, p u is the saturation density and p(x) is the average of the density p(x,y) with respect to the transverse direction y and folded in the x direction with a Gaussian of width a = 0.87 fm. Since the effective field U then depends on x only, the same holds for those modes that can be excited, and this simplifies the analysis considerably. The folding simulates the finite range of the interaction which introduces a natural cutoff scale for the unstable modes. In determining the parameters of J7, we have sought to mimic standard threedimensional matter. Thus we have required a Fermi momentum of Pp = 260 MeV/c, and a binding energy of 16 MeV per nucleon, and that a density doubling leads to approximately zero binding (corresponding to a compressibility modulus of K = 300 MeV for a calculation in three dimensions). The in-medium "cross-section" was taken as 2.4 fm, corresponding to a interaction radius of 1.2 fm. In the present idealized two-dimensional system, the saturation density is p lt = 0.55 fm~". In order to ensure that the system is initially situated inside the spinodal region, we prepare the system to have a uniform density equal to half the normal density, p(x,y,t = 0) = i/? 0 , and having a temperature of T -3 MeV. This physical scenario is similar to what was considered in ref. [7].

Lattice simulation

Recently a lattice simulation of the BL evolution has been developed [7, 8, 10 j. In this method the stochastic part of the collision number, SN(12; 1'2') (see above), is simulated directly on a lattice in phase space. The size of the phase-space cells must be of the order h 2 in order to simulate the extension of an elementary quantum state and we have used Aa; = 3 fm and Ap -120 MeV/c. The mean-field evolution is treated by means of a standard matrix technique. However, in order to achieve sufficient accuracy, a second finer grid is required and therefore smaller cells having 6x = I fm and 8p = 40 MeV/c are employed. The transformations between these two scales are described in ref. [7]. We first focus on the evolution of the spatial density of the system, which is displayed in fig. la. Initially the system has a uniform density, but soon the fluctuations break this translational symmetry. Subsequently the fluctuations are rapidly amplified by the action of the effective one-body field, thus leading towards fragment formation. (Of course real fragments do not form, because the size of the torus is kept fixed so the system cannot expand.) It is instructive to perform a Fourier analysis of the density for the considered system, p(x). Indeed, in the present simplified scenario, the eigenmodes of the density are plane waves, characterized by a wave number k. The amplitude associated with a given k is given by dx -ikx 6p(x,t), (7) where 6p(x) = p(x) -p 0 represents the fluctuating part of the density and L x is the size of the system in the x direction. The Fourier transform p(k,t) pertains to a particular dynamical history. In reality we consider an entire ensemble of such evolutions (generated by running the code many times) and it is therefore of interest to consider the following ensemble average,

= J / dxdx' Sp(x,t) Sp(x',t) ( 8
)
This quantity is also recognized as the Fourier transform of the spatial correlation function <r(x -x') =-< 6p(x,t) 8p(x',t) > [6].

Invoking linear response theory, we expect that the evolution of tr% is dominated by the associated unstable mode i/&, and so it should exhibit approximatly the behaviour described by eq. ( 4). In the present case, the eigenmodes factorize so F Vlc (x,p) is a plane wave with respect to the position x times a function of the momentum p.

In fig. 2a is shown the time evolution of af. for the fastest growing mode. We note that while <rf. is close to zero during the earliest time steps, the exponential increase predicted by eq. ( 4) is clearly observed at later times. For large times, t > tk, we then expect that <r| % Dktkexp(2t/tk) and we may then extract the key quantities Dk and tk by simple extrapolation. This procedure gives Dk ~ 2 • 10~"' fin""* and tk % 50 fm/c. These values are in perfect agreement with those obtained from eq. ( 3) for Dk and the dispersion relation (see ref. [9]) for tk-Finally, in fig. 3a we show the fluctuation <r\ as a function of the wave number k. It can be seen that certain modes are amplified more rapidly than others, in accordance with their respective characteristic times <*, and the final Fourier spectrum is therefore dominated by the most unstable modes.

The test particle method

In order to test the possibility of replacing the stochastic collision term by a simple noise, we solve the BUU equation with the test particle method, using a collisional procedure based on the concept of nucleon mean free path [11]. Since each nucleon is represented as a collection of J\f test particles, the stochasticity of the collision integral is reduced correspondingly by the factor I/A/". Therefore, since Af is always finite, the test-particle method retains some degree of fluctuation in the dynamics. (This is the feature which we wish to turn to our advantage, since it provides a natural source of initial noise.) When the spinodal region of the nuclear matter phase diagram is not explored small fluctuations of the distribution function are not important because they are not amplified (but destroyed), and consequently a good description of the mean dynamics is obtained with the test-particle method [2]. By contrast, when the system is unstable, all kinds of perturbations of the density / may be important, since they are amplified, and therefore the numerics is no longer under control. This is what we observe in fig. Ib where the initial noise is amplified and leads towards clusterization. In this situation our idea is to carefully adjust the numerical fluctuations (i.e. the number of test particles .V) so as to simulate the physical fluctuations, which are well accounted in the lattice calculations (see ref. [8] for a detailed discussion).

Since the perturbation of the density introduced by using a finite number of test particles is of statistical nature, it may be calculated in a simple manner. Thus the number of test particles n present in a volume V fluctuates with a characteristic dispersion a^t -n, where n = pVAf is the mean number of test particles in V, which is fixed by the specified value of the initial density. Consequently the relative fluctuation of the density p is <r p /p = l/-\/fi, or

_ _ = " n The corresponding Fourier coefficient (8) is then independent of k, r k - (9) (10)
Thus, the numerical noise <r~k scales as I/.A/", and so it is possible to choose the number of test particles A/" so as to reproduce the required magnitude of the physical source term Dk, for a particular value of the wave number k. Indeed, since in the BL lattice calculation we have tr%(0) = 0, it follows from eq. ( 4) that asymptotically (i.e. for t ^> tk) the magnitude of the fluctuations is determined by Dktki whereas it is given by Dktk/Af + 0"i(0) in the BUU test-particle simulation. Therefore, the matching of those two asymptotic evolutions yields the condition

-l)D k t k = (11)
which determines the number of test particles jV. Inserting the values determined above for D k and tjt, together with L x , L y , and /?, we find Af ~ 83 for the fastest mode (the main uncertainty coming from the relatively inaccurate determination of the diffusion coefficient Dk).

We have then performed BUU calculations using a range of values for Af and compared the results with the corresponding lattice calculation. The time evolution of the fastest mode is shown in fig. 2, in addition to the corresponding BL result. The results are seen to bear out our expectation that the time evolutions follow the form given in eq. ( 4). Moreover, we observe that the lattice calculation is indeed well reproduced by the BUU simulation when the value fif -90 is employed, which is in excellent agreement with the above result J\f % 83.

Figure Ib shows the evolution of one particular density distribution (calculated with the optimal value . A/* = 90), and it is seen that the density irregularities quickly attain the same magnitude as those of the lattice calculation.

In fig. 3 we show the quantity <r\ as a function of the wave number k for different times. We observe that not only are the time evolutions nearly the same, when comparing the BL lattice calculation and the BUU test-particle method, but also the spectral shape of the fluctuations are remarkably similar. This feature is v»ry if, if.

important for the viability of the method. The similarity of the spectral shapes arises from the fact that the product D^tk depends only relatively weakly on k, in the neighborhood of the fastest mode, so that the spectral profile depends primarily on the amplification times i/t which are the same in the two treatments.

t 4 Discussion

We have shown that when one is interested in the dynamics of very unstable systems the specific character of the fluctuations may not play a decisive role in the catastrophic evolution, provided the overall magnitude of the fluctuations is suitably adjusted. This is because the various unstable modes, once they have been agitated, are amplified exponentially so that the most unstable ones will tend to quickly dominate. In the present context, the magnitude of the imposed noise is set by the numerical parameter Af, the number of test particles per nucleon, which is determined by the relation (11). In the idealized two-dimensional case considered here, this number is around J\f % 90, which is close to the value J\f -83 pertaining to a full two-dimensional situation. We are presently in the process determining the value pertaining to three dimensions.

It should be noted that the source term D^ depends on the phase-space density /(r, p). For example, at relatively moderate excitations it is proportional to the square of the temperature. Consequently, the optimal value of JV", as given by eq. ( 11) will depend on the the process considered and should therefore be carefully adjusted for the particular scenario under study. This feature is particularly relevant in connection with the very recent calculations showing various exotic structures in central symmetric collisions [12,13,15,14,[START_REF] Batko | LBL-32722 (Sep'92)[END_REF], since these geometries depend sensitively on the time scale for the cluster formation, and hence the results will depend qualitatively on the particular value of Af employed.

The advantage of the proposed test-particle method is that it makes it practically possible to perform meaningful three-dimensional simulations of nuclear collisions, a task that is presently beyond feasibility for the exact Boltzmann-Langevin model. However, we wish to emphasize that the simple method explored here should only be regarded as an approximate method for obtaining physically reasonable simulation of nuclear catastrophic processes using the standard BU[/model.

Finally, it should be mentioned that the two-body collision integral is not the only source of fluctuations in nuclear dynamics. This fact is perhaps most evident at lower energies where the Pauli blocking renders the collision integral ineffectual, and in fact the fluctuations observed for damped reactions have been largely understood as associated with the random character of the individual nucleon transfers.[17] A suggestion is to introduce a fluctuation in the effective field, SH, which can be done already at the level of the Vlasov equation [START_REF] Bromley | [END_REF]. Hopefully, by confronting careful calculations based on well-controlled approximations with experimental data, it will ultimately be possible to ascertain the relative importance of the various possible mechanisms. CD O 0 0 0 q q q q q o d d d d 6

1 fFigure 1 :

 11 Figure 1: Density profiles. The density profile p(x,t) associated with one particular trajectory at four different times t, for the BL lattice calculation (part a) and for the BUUtest-particle simulation (part 6).

  s

Figure 2 :

 2 Figure 2: Growth of the fluctuations. The time evolution of the Fourier transform of the correlation function, a\, for the most unstable mode (which has node number K = kL^/lw = 6), for the BL lattice calculation (part o) and the BUU test-particle simulation with three different values of A/", the number of test-particles per nucleon (part 6). This quantity is the variance of fluctuations having wave number equal to k.

Figure 3 :

 3 Figure 3: Variance spectrum. The variance <TJ! as a function of the mode number K (-kL x /2w), shown at a number of times t for the BL lattice calculation (part a) and for the BUU test-particle simulation with Af = 90 (part 6).

"This work was supported in part by the Director, Office of Energy Research, Office of High Energy and Nuclear Physics, Nuclear Physics Division of the U.S. Department of Energy under Contract No. DE-AC03-76SF00098.


Pierre Hosteins
email: pierre.hosteins@ifsttar.fr

A Compact Mixed Integer Linear Formulation for Safe Set Problems

Keywords: Graph partitioning, network control, Safe Set, Connected Safe Set, Mixed Integer Linear Programming, symmetry-breaking, variable reduction

The Safe Set Problem is a type of partitioning problem with particular constraints on the weight of adjacent partition components. Several theoretical results exist in the literature along with a mixed integer linear formulation. We propose a new, more compact, Mixed Integer Linear Program for the (Weighted) Safe Set Problem and Connected Safe Set Problem with a polynomial number of variables and constraints, based on a flow from an additional fictitious node. The model is enriched by symmetrybreaking considerations and a variable reduction subproblem. We test the model on a benchmark of small instances which underline the difficulty of solving the Safe Set Problem through mathematical programming approaches. We also compare it with the only existing formulation in the literature and find that our approach is usually more efficient on a set of benchmark instances.

Introduction

Consider an undirected and connected graph G = (V, E) with node weights w i ≥ 0 for i ∈ V . We call N i ⊂ V the set of neighbours of node i ∈ V and define a non-empty subset S ⊂ V which induces a subgraph G[S]. The function w : 2 V → R + is defined as w(S) = i∈S w i , i.e. the total weight of a subset of nodes S ⊆ V . S is said to be a safe set if S = ∅ and for any maximally connected component C in G[S] and C in G[V \ S] we have w(C) ≥ w(C ) whenever an edge exists between C and C . The Weighted Safe Set problem (WSSP) consists in finding the safe set with minimum total weight, while the simpler Safe Set Problem (SSP) considers only unit node weights. The Connected Safe Set Problem (CSSP) is a Safe Set Problem where the safe set is connected (i.e. G[S] has only one component). An optimal solution of any Safe Set variant is denoted as S * . We observe that if G is not connected, the problem decomposes completely on each of the connected components of G. Thus, we assume without loss of generality that G is connected. An example of a feasible safe set on a small graph is illustrated in Figure 1.

The SSP was first defined in [START_REF] Bapat | Network majority on tree topological network[END_REF] as a way to determine important influential communities in a social network. Another application is invoked in [START_REF] Fujita | Safe set problem on graphs[END_REF] where one looks for a subspace inside a building that can shelter the people from nearby spaces in a case of a catastrophe, i.e. a building design problem. From a more general point of view, the SSP is designed to determine the smallest subset of a network which we need to control (or equip) in order to control the behaviour of (or assist) the whole network. One can see from its definition above that it can be classified as a type of graph partitioning problem, though only node weights are considered. This differentiates greatly the SSP from partitioning problems based on edge cuts. Other nodebased partitioning problems exist, such as the Vertex Separator Problem [START_REF] Balas | The vertex separator problem: a polyhedral investigation[END_REF] and the k-separator problem [START_REF] Ben-Ameur | The k-separator problem: polyhedra, complexity and approximation results[END_REF] where a subset of nodes of minimum cardinality (or weight) is removed from the graph in order to obtain (possibly balanced) disconnected components of the original graph. Another class of interdiction problem, called the Critical Node Problem [START_REF] Arulselvan | Detecting critical nodes in sparse graphs[END_REF][START_REF] Lalou | The critical node detection problem in networks: A survey[END_REF], is also linked to node-based graph partitioning as some of its versions aim at fragmenting the graph into disconnected subgraphs by deleting a subset S of nodes which is not necessarily a minimal node cut-set. Contrary to such node-cut based graph partitioning problems, no subset of nodes is removed from the graph in a Safe Set solution. Other social network problems study the influence of certain parts of the network on the whole network. A prominent example is the positive influence dominating set problem [START_REF] Dinh | On the approximability of positive influence dominating set in social networks[END_REF][START_REF] Wang | Positive influence dominating set in online social networks[END_REF], where we seek a minimal set of nodes S such that all other nodes in the network have at least a fraction ρ > 0 of their neighbours in S. Such models try to identify a minimal set of participants in the network to solicit so that we can guarantee a positive influence over the rest of the network. Such models, however, are defined based on local properties, such as a subset of neighbours for each node belonging to a certain set. The SSP is defined using more global properties based on whole connected components of a graph partition, which makes it less straightforward to cast as a Mixed Integer Linear Program (MILP) to solve it. The problem is quite recent and has not been studied extensively in the literature. Moreover mainly theoretical contributions on the subject can be found. For example, the weighted problem over trees has been shown to be NP-hard in [START_REF] Bapat | Network majority on tree topological network[END_REF][START_REF] Bapat | Safe sets, network majority on weighted trees[END_REF] while the unweighted problem is shown to be solvable in polynomial time on graphs with bounded treewidth and interval graphs in [START_REF] Águeda | Safe sets in graphs: Graph classes and structural parameters[END_REF]. Bounds on the size of safe sets and connected safe sets are also studied in [START_REF] Fujita | Safe set problem on graphs[END_REF]. While the present manuscript was being revised, the work of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] was published, which introduces a MILP formulation and a greedy heuristic to solve the WSSP on general graphs. The MILP has an exponential number of constraints, which are introduced at each node of the branching tree using a separation algorithm, and can be easily adapted to the connected version of the problem. The MILP and the heuristic are tested on a benchmark of instances with number of nodes ranging from 10 to 30 and different edge densities.

We propose in the following a relatively compact MILP with a polynomial number of variables and constraints to tackle the problem on general graphs. The model is quite different in its approach from the one in [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF]. In Section 2 we detail our MILP formulation for the WSSP and its connected version. In Section 3 we propose a way to break some problem symmetries and solve a simpler MILP in a preprocessing step to reduce the number of variables of our Safe Set model. In Section 4 we present the results of numerical experiments made on a benchmark of small random graphs and in Section 5 we briefly conclude. Additional results on small world graphs are provided in Appendix A. In Appendix Appendix B, we compare the results of our approach with that of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF].

A compact MILP formulation

In the following we define the necessary variables and constraints to provide a MILP for the (Connected) WSSP. We define a model that identifies the connected components of a given solution by using a flow circulation coming from a fictitious node 0, which will be partially absorbed by each node in the graph. Recall that a feasible safe set is a subset S ⊂ V , and we define n := |V |, m := |E|≥ n -1, and W := i∈V w i . We also define E to be the set of directed edges corresponding to E, i.e. E = {(i, j) : i ∈ V, j ∈ V and {i, j} ∈ E}. Since for a general graph, it is impossible to know in advance the number of maximally connected components in the optimal partition, we have to define variables with indices ranging up to the maximal theoretical number of connected components n-11 . As a result, the variables related to potential components which are unused in the optimal solution will be 0. In order to formulate our model, we will define the following variables, corresponding to a feasible safe set S:

• x i : binary variable equal to 1 if node i ∈ S (i.e. is in the safe set), 0 otherwise.

• y ij : binary variable equal to 1 for edge {i, j} ∈ E if i ∈ S and j ∈ S, 0 otherwise.

• y ij : binary variable equal to 1 for edge {i, j} ∈ E if i ∈ V \ S and j ∈ V \ S, 0 otherwise. • a ic : binary variable equal to 1 if node i ∈ S is part of component c ∈ {1, . . . , n -1}.
• a ic : binary variable equal to 1 if node i ∈ V \ S is part of component c ∈ {1, . . . , n -1}.

• t ic : binary variable equal to 1 if node i ∈ S is the unique node in c ∈ {1, . . . , n -1} which receives a non-zero flow from node 0.

• t ic : binary variable equal to 1 if node i ∈ V \ S is the unique node in c ∈ {1, . . . , n -1} which receives a non-zero flow from node 0.

• f ij : continuous flow variables from i ∈ V ∪ {0} to j ∈ V defined for couples (i, j) such that i = 0 or (i, j) ∈ E .

• ω i : total weight of the component in G[S] to which i belongs.

• ω i : total weight of the component in G[V \ S] to which i belongs.

• ν c : total weight of the component c ∈ {1, . . . , n -1}.

• ν c : total weight of the component c ∈ {1, . . . , n -1}.

Variables with a prime refer to quantities linked to the set V \ S. The approach advocated is inspired by a MILP formulation proposed for a Critical Node Problem in [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF] which aims at deleting a subset of nodes in the graph in order to maximise the number of maximally connected components of the resulting subgraph. The rationale behind this strategy is to make sure that the structure of the connected components inside the partitioned graph is correctly taken into account by the MILP. In the modified graph, node 0 is the source of n units of flow while all the other nodes in V are sinks which absorb one unit of net flow each. Moreover, each node is assigned to a connected component inside G[S] or G[V \ S] and only one node per connected component can receive a non-zero amount of flow from node 0. Finally, the flow received (and partially absorbed) by such a node is dispatched to the other nodes in the component until it is completely absorbed. Since a flow can pass only between two nodes connected by an edge {i, j} for which y ij = 1 or y ij = 1, the flow conservation conditions guarantee that we cannot assign two nodes to the same connected component if they are not indeed part of a unique connected component inside

G[S] or G[V \ S].
We illustrate this process in Figure 2 on a small graph where the optimal safe set is made up of the two grey nodes 1 and 2 (with weight w 1 = w 2 = 2) while the nodes in V \ S (with unit weights) are white. A total flow of 8 is created at node 0 and flows along the arcs displayed as dashed lines towards one node in the unique component of the safe set and the two nodes 3 and 6 belonging each to a component of G[V \ S] (the other arcs from node 0 to the remaining nodes in V are omitted for reasons of clarity). The arcs where a non-zero flow circulates are displayed as directed arcs in order to represent the direction of the flow.

We stress that instead of using a multi-commodity flow model as in [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF], we only use a single type of flow. Indeed, the model of [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF] is a bi-level MILP which needs an integrality property for the inner optimisation problem in order to dualise it and obtain a single level reformulation (see Proposition 1 in [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF]). Restricting to a single commodity in the framework of [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF] would nullify this property, while it is perfectly admissible in our case. This certainly introduces a less tight formulation on general grounds. However, it greatly reduces the number of flow variables and accelerates notably the resolution of the linear relaxation of the MILP presented below. Moreover, our tests with a multi-commodity formulation tend to hint at the fact that such a formulation does not in practice improve the lower bound of the model detailed below. Armed with the above variables, we propose the following formulation of the WSSP, called M SS :

(M SS ) minimise i∈V w i x i (1) subject to y ij ≥ x i + x j -1 {i, j} ∈ E, (2)
y ij ≥ 1 -x i -x j {i, j} ∈ E, (3)
y ij ≤ x i {i, j} ∈ E, (4)
y ij ≤ x j {i, j} ∈ E, (5)
y ij ≤ 1 -x i {i, j} ∈ E, (6)
y ij ≤ 1 -x j {i, j} ∈ E, (7)
j∈N i f ij - j∈N i ∪{0} f ji = -1 i ∈ V, (8)
i∈V

f 0i = n (9)
f 0i ≤ (n -1) n-1 c=1 (t ic + t ic ) i ∈ V, (10)
f ij ≤ (n -1)(y ij + y ij ) (i, j) ∈ E , (11)
n-1 c=1

a ic = x i i ∈ V, (12)
n-1 c=1

a ic = 1 -x i i ∈ V, (13)
a ic ≥ a jc + y ij -1 {i, j} ∈ E, c = 1, . . . , n -1, (14)
a jc ≥ a ic + y ij -1 {i, j} ∈ E, c = 1, . . . , n -1, (15)
a ic ≥ a jc + y ij -1 {i, j} ∈ E, c = 1, . . . , n -1, (16)
a jc ≥ a ic + y ij -1 {i, j} ∈ E, c = 1, . . . , n -1, (17)
t ic ≤ a ic i ∈ V, c = 1, . . . , n -1, (18)
t ic ≤ a ic i ∈ V, c = 1, . . . , n -1, ( 19
) i∈V t ic ≤ 1 c = 1, . . . , n -1, ( 20
) i∈V t ic ≤ 1 c = 1, . . . , n -1, (21)
ν c = i∈V w i a ic c = 1, . . . , n -1, (22)
ν c = i∈V w i a ic c = 1, . . . , n -1, (23)
ω i ≥ ν c -W (1 -a ic ) i ∈ V, c = 1, . . . , n -1, (24)
ω i ≤ ν c + W (1 -a ic ) i ∈ V, c = 1, . . . , n -1, (25)
ω i ≤ W x i i ∈ V, (26)
ω i ≥ ν c -W (1 -a ic ) i ∈ V, c = 1, . . . , n -1, (27)
ω i ≤ ν c + W (1 -a ic ) i ∈ V, c = 1, . . . , n -1, (28)
ω i ≤ W (1 -x i ) i ∈ V, (29)
ω i ≥ ω j -W y ij (i, j) ∈ E , (30)
i∈V x i ≥ 1 (31) y ij , y ij ∈ {0, 1} : {i, j} ∈ E, (32)
a ic , a ic ∈ {0, 1} : i ∈ V, c = 1, . . . , n -1, ( 33
)
ω i , ω i ≥ 0 : i ∈ V, (34)
ν c , ν c ≥ 0 : c = 1, . . . , n -1, (35)
f ij ∈ [0, n -1] : i ∈ V ∪ {0}, j ∈ V, (36)
x i ∈ {0, 1} : i ∈ V, (37)
t ic , t ic ∈ {0, 1} : i ∈ V, c = 1, . . . , n -1. ( 38
)
Should we want to extend the problem to the connected version, we must define the variables t ic , a ic and ν c only for c = 1, or else add the following constraint:

i∈V n-1 c=2 a ic = 0. (39)
Let us now explain briefly each constraint family used to model the WSSP. Constraints (2) to ( 7) guarantee the consistency of y variables with respect to the decision variables x. They are a linearisation of constraints [START_REF] Bapat | Safe sets, network majority on weighted trees[END_REF] guarantees that a unit of flow coming from node 0 is absorbed by each node in the graph. Constraint [START_REF] Ben-Ameur | The k-separator problem: polyhedra, complexity and approximation results[END_REF] guarantees that exactly the right amount of flow is provided by fictitious node 0. Constraints [START_REF] Codato | Combinatorial Benders' cuts for mixed-integer linear programming[END_REF] ensure that flow can only go from node 0 to another node in the graph if the associated t or t variable is non-zero. Constraints [START_REF] Cornaz | The vertex k-cut problem[END_REF] then make sure that the flow can only circulate on edges inside a given component of 12) and ( 13) assign a node i ∈ V to only one component, be it in

y ij = x i x j and y ij = (1 -x i )(1 -x j ) for {i, j} ∈ E. Constraint
G[S] or G[V \ S]. Constraints (
G[S] or G[V \ S]
. Constraints ( 14) to [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF] then ensure that two nodes inside the same component are indeed assigned to the same component number c. Indeed, when y ij = 1, the two inequalities ( 14) and ( 15) will transform into the equality a ic = a jc for all values of c (the same is true for ( 16) and ( 17) when y ij = 1). We then define the consistency relation of the t variables with respect to variables a in constraints ( 18) and ( 19) and also that only one node inside component c of

G[S] or G[V \ S]
can have a non-zero t value. This node will represent the bridge through which some flow will enter component c and then be dispatched to all nodes of c, each of which needs to absorb one unit of net flow. The combination of all previous constraints impose that for each partition c of G[S] (respectively G[V \ S]), the necessary quantity of flow will be dispatched to c through a single bridge node, and will then be absorbed by all the nodes in the component. Then, the link between t and a variables ensures that the nodes in this component will be assigned to the corresponding component number.

Constraints [START_REF] Watts | Collective dynamics of small-world networks[END_REF] and ( 23) are used to define the ν variables according to which node is assigned to each component label c and constraints (24) to (29) ensure that the ω variables are defined correctly with respect to the ν variables. Recalling that ω i is the total weight of the component that node i ∈ S belongs to, constraints (30) are the Safe Set constraints which guarantee that ω i ≥ ω j if {i, j} ∈ E with i ∈ S and j ∈ V \ S. Finally, since constraint (30) allows for the model to choose to assign every node to V \ S (and therefore S = ∅), we force the model to have at least one node in S with Equation (31).

Proposition 1. Model M SS provides the optimal solution for the (Connected) WSSP.

Proof. Consider first an optimal solution S * to the WSSP. We can construct a feasible solution to model M SS as follows. First define x * i = 1 for all nodes i ∈ S * and x * i = 0 for i ∈ V \ S * . It follows that y * ij = 1 for i ∈ S * and j ∈ S * with {i, j} ∈ E, while y * ij = 0 if either i / ∈ S * or j / ∈ S * (and conversely y * ij = 1 if i ∈ V \S * and j ∈ V \ S * with {i, j} ∈ E and 0 otherwise), so that constraints (2) to ( 7 12) to [START_REF] Wang | Positive influence dominating set in online social networks[END_REF]. For each index c ∈ {1, . . . , n -1} attributed to a connected component C ∈ C S , we set variable ν * c to the total weight of the nodes of that component and for each node i ∈ S we set ω * i to the total weight of the component that node i belongs to (and we perform the same actions for the components of C S and the variables ν * c and ω * i ). Finally, we can set the flow variables as follows. For each node i ∈ V and c ∈ {1, . . . , n -1} such that t ( ) * ic = 1, we set f * 0i to be equal to the number of nodes in the connected component of node i. After which, we increase f * jk by one unit over each edge {j, k} on a path from node i to each node of its connected component. It is easy to verify that constraints [START_REF] Bapat | Safe sets, network majority on weighted trees[END_REF] to [START_REF] Cornaz | The vertex k-cut problem[END_REF] are satisfied by such assignments. We therefore have constructed a solution to model M SS with an objective function value which is equal to w(S * ) and therefore is not larger than the one of optimal solution S * to the WSSP.

Consider now the optimal solution to model M SS , with binary optimal values for x * i for i ∈ V . We define the safe set S * = {i ∈ V : x * i = 1}. By construction, w(S * ) is the same as the optimal objective function value of M SS so we have constructed a solution to the WSSP with an objective function value which is not larger than the one of the optimal solution to model M SS , which completes the proof.

Proposition 2. The y, y , a and a variables can be relaxed as continuous and still achieve integrality when the x and t are integers.

Proof. Suppose the variables y ( ) and a ( ) were defined as continuous variables bounded by 0 and 1. Consider a feasible solution to M SS with integer x and t ( ) variables. We have a partition of V into S = {i ∈ V : x i = 1} and its complement V \ S = {i ∈ V : x i = 0}. Consequently, the integrality of y and y is straightforward to check. Indeed, y ij = 0 when either x i = 0 or x j = 0 because of constraints ( 4) and [START_REF] Arulselvan | Detecting critical nodes in sparse graphs[END_REF]. In contrast, by [START_REF] Aringhieri | A Genetic Algorithm for a class of Critical Node Problems[END_REF], y ij = 1 when x i = x j = 1 (while the same reasoning applies to the y ).

We now prove the integrality of the a variables; the integrality of the a variables follows by a similar argument. For two nodes i ∈ S and j ∈ S with {i, j} ∈ E, constraints ( 14) and [START_REF] Lalou | The critical node detection problem in networks: A survey[END_REF] will ensure that a ic = a jc for all c = 1, . . . , n -1 and by propagation the same equalities will hold between all pairs of nodes inside the same connected component. By constraint [START_REF] Ben-Ameur | The k-separator problem: polyhedra, complexity and approximation results[END_REF] and flow conservation constraints [START_REF] Bapat | Safe sets, network majority on weighted trees[END_REF], a total flow of n units originates from node 0 and each node has to absorb one unit, so that some non-zero flow f 0i > 0 is required for a subset of nodes i ∈ V . Therefore, by constraints [START_REF] Codato | Combinatorial Benders' cuts for mixed-integer linear programming[END_REF] and [START_REF] Veremyev | An integer programming framework for critical elements detection in graphs[END_REF], for each component index c with i∈V a ic > 0 there will be one (and only one) node i with t ic = 1 and f 0i > 0. For that node, we see from [START_REF] Simonetti | The minimum connected dominating set problem: Formulation, valid inequalities and a branch-and-cut algorithm[END_REF] that a ic = 1 and a id = 0 for all d = c, so that i is assigned to component label c. By the reasoning applied on the a variables earlier, we can extend this result to all the nodes inside the same connected component as i, which guarantees that the a variables are binary whenever the x and t variables are. This completes the proof of integrality for y ( ) and a ( ) variables.

When implementing the model for numerical experiments, the y and y variables are kept continuous but we noticed that the model is solved faster if the a and a variables are kept binary. This unexpected behaviour might be due to the important limitations imposed by setting an a variable to 0 or 1, through a propagation effect of constraints ( 14) to [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF]. In such a case, branching early on a and a might allow for a larger increase of the optimal objective value of the linear relaxation when M SS is solved by a linear solver and speed up the overall convergence of the branching procedure. Moreover, we observe that CPLEX is able to eliminate some a ( ) variables during its presolving step, which might also explain why the reduced model is then easier to solve. In the above MILP we have in total O(n 2 ) binary variables and O(m) continuous variables (or O(n 2 ) continuous variables if the a and a variables are treated as continuous) and we have a maximum of O(mn) constraints, which overall makes the model relatively compact for sparse graphs.

MILP improvements

In this section we provide some improvements which will help solving the M SS faster when using a numerical MILP solver.

Symmetry breaking

The model defined in Section 2 suffers from several symmetries which might hamper its performance when we try to solve it with a MILP solver. For example, any node inside a given component c can be used equivalently to dispatch the flow from node 0 to the rest of the component. We can break such symmetry with the following constraints:

j>i t jc ≤ t ic + 1 -a ic , i ∈ V, c = 1, . . . , n, (40a)
j>i

t jc ≤ t ic + 1 -a ic , i ∈ V, c = 1, . . . , n. (40b)
Such constraints ensure that the node with the smallest index will be the one chosen to receive a non-zero flow from node 0.

We may also observe that for a given solution, a component in the graph can be given a priori any label c = 1, ..., n -1. We can try to break (partially) that symmetry through the following constraints:

ν c ≤ ν c-1 , 2 ≤ c ≤ n -1, (41a)
ν c ≤ ν c-1 , 2 ≤ c ≤ n -1. ( 41b
)

Reducing the number of variables

Given the large initial range for the c index of a and a , M SS has a large number of binary variables, most of which are superfluous in many instances. Moreover, for the linear relaxation of M SS , the initial large number of c indices for a allows for many small fractional values of these variables, suggesting that the nodes are partially "spread" inside many components and can satisfy easily the Safe Set constraint (30).

Limiting the range of available c indices for the nodes in a generic safe set or its complement could help increase the lower bound of the linear relaxation. In the following, we introduce a model whose optimal solution allows us to find an upper bound on the c index for all the variables of M SS .

In order to do so, we consider the "worst case", for any given safe set S, in terms of the number of connected components of G[V \ S]. We prove below that an upper bound on the number of connected components in both G[S] and G[V \ S] can be obtained by finding the subset Σ * ⊂ V of smallest cardinality |Σ * |= k * with at least one edge towards each node of V \ Σ * , which is the definition of the dominating set of minimum cardinality (which is an NP-hard problem [START_REF] Garey | Computers and Intractability: A Guide to the Theory of NP-Completeness[END_REF]). A classic MILP formulation for the Minimum Dominating Set Problem (see e.g., [START_REF] Simonetti | The minimum connected dominating set problem: Formulation, valid inequalities and a branch-and-cut algorithm[END_REF]) uses the following variables:

• x i = 1 if node i ∈ Σ and 0 otherwise.
We also define the set Γ i = N i ∪ {i} which contains node i ∈ V and its neighbours. k * is the optimal objective value of the following optimisation problem, called ∆ SS : Proof. Consider the graph G[V \S] and its connected components. If one such component has more than one node, it is possible to move some of these nodes to the safe set and obtain more (in the worst case, not less) connected components in G[V \ S]. By induction, a solution with the largest possible number of components in G[V \ S] is obtained when each such component has cardinality one. We obtain an upper bound on the number of connected components by ignoring the edges inside G[V \ S], so that each node in V \ S defines an independent connected component, and by relaxing the safe set constraints (30), supposing it is satisfied by any solution we select. However, since each node in G[V \S] has to define a possible independent component, it has to be connected by an edge to a node in the safe set S. Consequently, the problem is akin to finding a dominating set Σ * of minimum cardinality. Moreover, since k * represents the number of nodes attributed to the safe set and each node in

(∆ SS ) minimise k (42a) subject to i∈V x i = k (42b) j∈Γ i x j ≥ 1 i ∈ V, (42c)
k ∈ N, x i ∈ {0, 1} : i ∈ V. (42d
G[V \ Σ * ] is supposed to form a unique component, n -k * is the number of components in G[V \ Σ * ].
We can apply the same reasoning symmetrically to the components of G[S], which completes the proof.

Although the Minimum Dominating Set Problem was proved to be NP-hard, it is solved very quickly for the benchmark instances presented in Section 4 and Appendix A, which makes it an acceptable pre-processing step for the SSP. Using the optimal solution k * of ∆ SS allows us to define variables a Alternative ways to obtain an upper bound on the number of components of the safe set or its complement can be devised. For example, the Critical Node Problem based on removing a subset of nodes of maximum cardinality in order to maximise the number of connected components in the induced subgraph, see e.g. [START_REF] Shen | Exact interdiction models and algorithms for disconnecting networks via node deletions[END_REF][START_REF] Veremyev | An integer programming framework for critical elements detection in graphs[END_REF], can provide the maximum number of connected components in any induced subgraph if one removes the cardinality constraint on the set of deleted nodes. By extension, and since we are interested in the number of components and not their cardinality, we can extend the solution of such a modified Critical Node Problem by deleting nodes further until all remaining components have cardinality one, in which case we would obtain a vertex cover of the graph. Ergo, the solution to the Minimum Vertex Cover Problem would also provide a valid bound on the number of connected components of the safe set and its complement. Such problems, however, tend to be slower to solve than the Minimum Dominating Set Problem and therefore subtract more computational time to model M SS , since the total time budget of our experiments is used for the preprocessing step and the resolution of model M SS together. Preliminary results on select instances, using a relaxed version of the Critical Node Problem, did not indicate a significant improvement with respect to using the Minimum Dominating Set Problem.

Numerical experiments

We now propose to test the model M SS with the improvements of Sections 3.1 and 3.2 on a set of small instances in order to test its viability.

Instance generation and experimental conditions

In order to test our MILP, we generate a first set of benchmark instances, with n ∈ {20, 25, 30, 35, 40, 50} with different possible values for the edge density δ e = 2m/(n(n -1)) ∈ {0.1, 0.2, 0.3, 0.4}. For each value of n and δ e , we generate five instances randomly as follows. First we randomly add edges to create a tree, to ensure that the graph is connected, then we add edges randomly inside the graph until we reach the desired edge density. The graphs share similar properties with the random graphs of the Erdös-Renyi type, except that the constraint of building a connected graph can influence the graph characteristics for small sizes, e.g. the graphs with 10 nodes and density 0.1 have to be trees. The average characteristics of the graph subsets for each value of n and δ e are displayed in Table 1. They are respectively the clustering coefficient C, the average length of all shortest paths l, the diameter D, the average degree d, the dispersion of the nodes degree σ d and finally the ratio of the last two quantities which gives an idea whether the nodes have similar degree or whether there is a large degree dispersion in the graph. The measure of dispersion used the standard deviation. These same instances can be considered either as unweighted or with node weights randomly generated as integers between 1 and 10. Additional results on instances of the small world type are provided and discussed in Appendix A, while a comparison of M SS with the branch-and-cut approach of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] on a different set of benchmark instances is provided in Appendix B.

We add the symmetry-breaking constraints (40) to model M SS . The symmetry-breaking constraints (41) make the linear relaxation much slower and do not allow for any speedup when solving M SS , therefore we do not include them in our numerical experiments. The indices for the a and a variables are reduced using the variable reduction technique introduced in Section 3.2 by solving ∆ SS as a preprocessing step, which is always very fast. For the connected version of the problem, we impose that variables a and t are defined only for index c = 1. We stress that in such a case, the variable reduction technique can still be applied to the complement of the safe set.

n δ e C l D d σ d σ d / d n δ e C l D d σ d σ d / d 20 0
The code is implemented using the C++ library of CPLEX 12.9. Results are obtained on a computer with an Intel Core i7-6600U CPU at 2.60GHz (4 threads) and 32GB of memory. We set CPLEX's absolute gap to 0 and impose a time limit of 3600 seconds (minus the time needed to solve the variable reduction problem, which in practice is almost negligible). We also disable CPLEX's cuts, because they tend to slow down the solver and do not substantially improve the lower bound on the objective value.

Numerical results

We start this section with a numerical study of the effect of the MILP improvements of Section 3. The results are displayed in Table 2 for the instances of the first benchmark where each line represents an instance size and we limit ourselves to the smallest and highest density values δ e = 0.1 and δ e = 0.4 to assess the efficiency of the improvements both for our least and most dense instances. For each instance size we present the average time (in seconds) to solve the problem over the five corresponding instances, the average relative gap (in percent) normalised by the upper bound (i.e. comprised between 0 and 100%) and the number of instances (out of five) that were solved to optimality within the time limit. We include the instances which reached the time limit in the average solution time. The results were computed using the unweighted version of the SSP, respectively by solving M SS alone, by solving M SS with the symmetry constraints (40) and finally by solving M SS with the variable reduction technique of model ∆ SS . From the table, for δ e = 0.1, it appears that each improvement method brings by itself a reduction in the average time needed to solve an instance of the SSP, however the symmetry-breaking technique manages to close more instances than the variable reduction preprocessing. The results when using both the symmetry breaking constraints and the variable reduction technique are instead displayed in the upper left part of Table 3: the added benefit of using both techniques simultaneously is a further increase in the number of solved instances. For denser instances with δ e = 0.4 however, the benefit of the improvements is less clear. More precisely, the addition of the symmetry-breaking constraints seems to be beneficial but the variable reduction mechanism does not provide a clear help. As a consequence, the addition of both MILP improvements does not improve as clearly over the results of M SS without any improvement, though the number of closed instances is larger. This is not necessarily surprising since model M SS already shows better performances on denser instances without the improvements. Average results for the unweighted SSP over sets of five instances for edge density δ e = 0.1 or 0.4 and graph size n. We display results without improvements (M SS ), with the help of symmetry-breaking constraints (M SS +(40)) and of the variable reduction technique (M SS + ∆ SS ). For each graph size, we display the average running time of the model, the average relative gap and the number of instances out of five that are solved within the one hour time limit. The time is displayed in seconds.

M SS M SS +(40) M SS + ∆ SS n δ e time (
We display the results for the unweighted and weighted version of the SSP, with the help of both MILP improvements, in Table 3. The quantities reported in the table are the same as in Table 2. The model is generally able to find the optimal solution for instances of up to 30 nodes. Additionally, instances with a higher density of edges seem to be easier to solve as the average time is generally smaller for each instance size and a larger number of instances are solved within the time limit. This is in sharp contrast to some other types of node partitioning problems, see e.g. [START_REF] Aringhieri | Local Search Metaheuristics for the Critical Node Problem[END_REF]. Therefore, it seems to be a particular feature of our model to be fit for graphs with a somewhat large density, which are generally found in social networks applications (see e.g. some examples of real social networks intances and their characteristics in [START_REF] Aringhieri | A general Evolutionary Framework for different classes of Critical Node Problems[END_REF]). We observe that it is hard to conclude from our results whether it is harder to solve the weighted version of the problem, although it generally provides a larger average relative gap for larger instances, which is not a surprise in the light of complexity results over trees [START_REF] Bapat | Network majority on tree topological network[END_REF][START_REF] Bapat | Safe sets, network majority on weighted trees[END_REF]. Though our model can hope to solve instances of up to 40 nodes depending on the graph density, it is clear that solving instances with 50 nodes is out of reach within a limited computational time of one hour.

The results for the Connected Safe Set are displayed in Table 4 and generally confirm our previous findings. The dependance of the solve time on the edge density is harder to establish though, even if we still solve more instances with δ e = 0.4 than in other cases. It is generally faster to solve the connected version of the problem, which is probably due to the reduction of a and t variables in the model. We can solve more instances than the unconnected case within the one hour time limit. Nevertheless no instance of 50 nodes can be solved either for the C(W)SSP.

Conclusion

We have provided a compact MILP for the (Connected) Weighted Safe Set problem, a somewhat new partitioning problem in the literature with interesting applications in network majority control. The model is relatively compact for sparse graphs. We also provided symmetry-breaking constraints and a variable reduction technique which greatly help the convergence of a numerical solver on small instances. Our results show that instances with up to 35 nodes can be solved within one hour of computational time (depending on the density), which makes it faster than the branch-and-cut algorithm from [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF].

A first possible research direction to tackle the SSP through Mathematical Programming would be to study 3.

alternative improvements to our model, such as tightening further the number of variables or reducing further the symmetries. Given the number of continuous variables introduced to define correctly the connected components of a given solution, it would be interesting to pursue the exploration of our model through Benders decomposition and check whether it can help increasing the size of instances which can be solved. Since the continuous variables are easy to determine once the integer variables are fixed, it could be possible to solve the dual of the slave problem analytically to speed up the convergence of the Benders approach. Since several inequalities in M SS involve big-M coefficients, e.g. equations (24) to (30), alternative decomposition methods tailored for such cases, as in [START_REF] Codato | Combinatorial Benders' cuts for mixed-integer linear programming[END_REF], are also interesting future research directions. Another type of decomposition known to be successful for some partitioning problem is column generation. An approach along the lines of [START_REF] Cornaz | The vertex k-cut problem[END_REF] could provide interesting results on instances larger than those tackled in this work. Given the difficulty of solving even small instances of the Safe Set problem, designing tailored efficient heuristic approaches would be a very welcome avenue of research. For example, carefully designed metaheuristics have had great success for solving difficult node partitioning problems in the literature where exact approaches have failed to obtain good quality solutions [START_REF] Aringhieri | A Genetic Algorithm for a class of Critical Node Problems[END_REF][START_REF] Aringhieri | A general Evolutionary Framework for different classes of Critical Node Problems[END_REF][START_REF] Aringhieri | Local Search Metaheuristics for the Critical Node Problem[END_REF][START_REF] Zhou | Memetic search for identifying critical nodes in sparse graphs[END_REF]. As for the results in Section 4, we first display the effect of the MILP improvements for small world graphs in Table 6, for unweighted instances with parameter k = 4 and k = 10. The results are similar to those of Section 4, in the sense that instances which are less dense tend to profit from the MILP improvements while denser instances are not solved faster on average, while there even are instance sizes where the improvements slightly worsen the performance of the model. This underlines the dependency of the performance of said improvements on the graph topology. 40) and the variable reduction preprocessing for the Connected SSP on small world graphs, both unweighted and weighted. The quantities reported are the same as in Table 7.

n k δe C l D d σ d σ d / d n k δe C l D d σ d σ d / d 20
M SS M SS +(40) M SS + ∆ SS n k time (
however the LB was not reported.

The results for the model of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] are those reported in the tables of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF]. Even though those results were obtained on a different computer and using a somewhat different framework (e.g., we make use of the heuristics of CPLEX, while [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] does not), we believe the important quantitative differences in the running times of both models is enough to assert the relative computational efficiency of our approach. Indeed, in most instances with 20 nodes or more, model M SS closes the instances at least ten times faster. On the smallest instances, and contrary to the numerical analysis in Section 4 and Appendix A, M SS does not always solve denser instances faster for a given graph size. The same pattern, however, can be seen again on the largest unweighted instances. Table 9: UB and running time for the benchmark instances of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] for unweighted instances with densities 0.3, 0.5 and 0.7, for model M SS with improvements and the branch-and-cut algorithm of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF].

Figure 1 :

 1 Figure 1: Example of a feasible safe set for a small instance of the SSP. The nodes in the safe set are displayed in black.

Figure 2 :

 2 Figure 2: Example of a feasible flow configuration for a small instance of the WSSP. The nodes in the safe set are displayed in grey.

  ) are satisfied. Consider the sets C S and C S defined respectively as the connected components of graphs G[S * ] and G[V \ S * ]. For each component C ∈ C S , pick one node i ∈ C and a yet unused index c ∈ {1, . . . , n -1} and set a * ic = 1, as well as t * ic = 1. For nodes j of the same component set a * jc = 1 for the same index c and t * jc = 0, then set all other variables a * , a * , t * and t * related to the nodes of C at value 0. The same thing can be done for the nodes of each component C ∈ C S by changing the roles of a and t variables with a and t variables. Given the values of the y * ij and y * ij variables for arcs {i, j} inside each connected component, the values assigned clearly respect the constraints (

) Proposition 3 .

 3 The optimal value k * for the model (42) is such that n -k * defines an upper bound on the number of maximally connected components in G[S] and G[V \ S] respectively for any solution set S of the SSP.

c

  for the restricted set of indices c = 1, . . . , n -k * .

Table 1 :

 1 Average characteristics (explained in the text) for the sets of five instances, for each edge density δ e and graph size n.

	.1	0.00 7.00 19.00	1.90 1.09 0.58	20 0.2	0.15 2.24 4.40	3.80	2.48 0.65
	25 0.1	0.08 4.53 11.00	2.40 1.49 0.62	25 0.2	0.14 2.10 4.00	4.80	2.98 0.62
	30 0.1	0.05 3.25	6.80	2.87 1.88 0.66	30 0.2	0.19 2.06 3.80	5.80	3.71 0.64
	35 0.1	0.06 3.07	6.60	3.37 2.24 0.66	35 0.2	0.18 1.99 3.60	6.80	4.25 0.62
	40 0.1	0.10 2.79	5.40	3.90 2.54 0.65	40 0.2	0.20 1.98 3.60	7.80	4.89 0.63
	50 0.1	0.08 2.59	4.80	4.88 3.22 0.66	50 0.2	0.19 1.91 3.00	9.80	6.10 0.62
	20 0.3	0.27 1.81	3.00	5.70 3.58 0.63	20 0.4	0.39 1.62 3.00	7.60	4.69 0.62
	25 0.3	0.29 1.77	3.00	7.20 4.53 0.63	25 0.4	0.39 1.60 2.60	9.60	5.73 0.60
	30 0.3	0.29 1.74	3.00	8.67 5.28 0.61	30 0.4	0.39 1.60 2.60 11.60	6.95 0.60
	35 0.3	0.30 1.73	3.00 10.17 6.23 0.61	35 0.4	0.39 1.60 2.40 13.60	8.00 0.59
	40 0.3	0.30 1.72	3.00 11.70 7.15 0.61	40 0.4	0.40 1.60 2.40 15.60	9.28 0.60
	50 0.3	0.30 1.71	3.00 14.68 8.94 0.61	50 0.4	0.40 1.60 2.20 19.60 11.65 0.59

Table 2 :

 2

		s) gap (%) # solved	time (s) gap (%) # solved	time (s) gap (%) # solved
	20 0.1	453.6	0.00	5	273.6	0.00	5	198.2	0.00	5
	25 0.1	3347.0	40.69	1	2291.8	24.18	3	2598.8	25.78	2
	30 0.1	3117.2	30.11	2	1932.2	11.60	4	3037.2	33.98	2
	35 0.1	3600.0	74.57	0	3600.0	48.55	0	3600.0	52.14	0
	40 0.1	3600.0	80.29	0	3600.0	72.76	0	3600.0	76.40	0
	50 0.1	3600.0	92.94	0	3600.0	91.81	0	3600.0	90.84	0
	20 0.4	15.4	0.00	5	18.0	0.00	5	14.0	0.00	5
	25 0.4	103.6	0.00	5	95.6	0.00	5	65.6	0.00	5
	30 0.4	327.2	0.00	5	334.2	0.00	5	237.6	0.00	5
	35 0.4	2782.0	10.16	2	1632.4	0.00	5	2546.8	6.72	2
	40 0.4	3322.4	8.33	3	2672.0	0.00	5	3380.8	16.12	3
	50 0.4	3600.0	95.97	0	3600.0	95.15	0	3600.0	93.77	0

Table 3 :

 3 Average results of model M SS with symmetry breaking constraints (40) and the variable reduction preprocessing for the SSP, over sets of five instances for each edge density δ e and graph size n. For both unweighted and weighted instances, we display the average running time of the model, the average relative gap and the number of instances out of five that are solved within the one hour time limit. The time is displayed in seconds.

			Unweighted			Weighted
	n δ e	time (s) gap (%) # solved	time (s) gap (%) # solved
	20 0.1	94.2	0.00	5	72.2	0.00	5
	25 0.1	532.2	0.00	5	746.4	0.00	5
	30 0.1	1743.8	4.44	4	2554.6	16.65	3
	35 0.1	3538.0	33.84	1	3600.0	39.53	0
	40 0.1	3600.0	67.73	0	3600.0	74.89	0
	50 0.1	3600.0	87.97	0	3600.0	94.55	0
	20 0.2	15.4	0.00	5	28.0	0.00	5
	25 0.2	234.4	0.00	5	140.2	0.00	5
	30 0.2	856.8	0.00	5	1472.0	7.52	4
	35 0.2	3437.0	18.01	2	3522.6	34.49	1
	40 0.2	3600.0	71.13	0	3600.0	72.21	0
	50 0.2	3600.0	93.11	0	3600.0	93.72	0
	20 0.3	14.4	0.00	5	11.8	0.00	5
	25 0.3	77.6	0.00	5	52.0	0.00	5
	30 0.3	362.6	0.00	5	562.2	0.00	5
	35 0.3	2798.2	10.64	3	1956.8	4.31	4
	40 0.3	3600.0	38.39	0	3384.0	43.19	1
	50 0.3	3600.0	91.05	0	3600.0	96.36	0
	20 0.4	14.2	0.00	5	12.8	0.00	5
	25 0.4	98.8	0.00	5	66.4	0.00	5
	30 0.4	267.8	0.00	5	375.0	0.00	5
	35 0.4	1709.4	2.40	4	905.4	0.00	5
	40 0.4	2678.4	11.58	4	2788.0	16.05	3
	50 0.4	3600.0	95.27	0	3600.0	99.41	0

Table 4 :

 4 Results of model M SS with symmetry breaking constraints (40) and the variable reduction preprocessing for the Connected SSP, on both unweighted and weighted instances. The quantities reported are the same as in Table

			Unweighted			Weighted
	n δ e	time (s) gap (%) # solved	time (s) gap (%) # solved
	20 0.1	4.4	0.00	5	4.6	0.00	5
	25 0.1	41.4	0.00	5	55.0	0.00	5
	30 0.1	258.4	0.00	5	284.4	0.00	5
	35 0.1	2423.0	29.30	2	1614.2	11.05	4
	40 0.1	3492.8	38.60	1	3529.6	41.95	1
	50 0.1	3600.0	66.30	0	3600.0	74.76	0
	20 0.2	7.8	0.00	5	6.8	0.00	5
	25 0.2	47.2	0.00	5	35.4	0.00	5
	30 0.2	232.2	0.00	5	193.8	0.00	5
	35 0.2	3206.2	17.22	2	1687.4	5.62	4
	40 0.2	3600.0	42.36	0	3600.0	41.44	0
	50 0.2	3600.0	67.66	0	3600.0	76.88	0
	20 0.3	7.2	0.00	5	6.4	0.00	5
	25 0.3	48.6	0.00	5	51.2	0.00	5
	30 0.3	256.4	0.00	5	217.6	0.00	5
	35 0.3	2265.2	0.00	5	1300.6	0.00	5
	40 0.3	3600.0	45.32	0	3600.0	45.44	0
	50 0.3	3600.0	92.72	0	3600.0	83.52	0
	20 0.4	9.6	0.00	5	9.8	0.00	5
	25 0.4	71.8	0.00	5	74.8	0.00	5
	30 0.4	236.0	0.00	5	302.4	0.00	5
	35 0.4	2485.0	2.57	4	1002.4	0.00	5
	40 0.4	3297.0	18.94	2	2915.0	0.00	5
	50 0.4	3600.0	92.60	0	3600.0	95.21	0

Table 5 :

 5 Average characteristics for the sets of five instances, for the small world graphs. The columns represent the same quantities as in Table1except for the additional column δe which represents the average edge density over the five instances.

	4	0.21 0.46 2.72 5.00 4.00 0.37 0.09	20 6	0.32 0.53 1.98 3.80	6.00 0.48 0.08
	25 4	0.17 0.45 3.15 6.00 4.00 0.35 0.09	25 6	0.25 0.55 2.33 4.20	6.00 0.52 0.09
	30 4	0.14 0.45 3.54 7.60 4.00 0.39 0.10	30 6	0.21 0.54 2.63 5.00	6.00 0.44 0.07
	35 4	0.12 0.44 3.77 7.80 4.00 0.42 0.11	35 6	0.18 0.54 2.85 5.80	6.00 0.52 0.09
	40 4	0.10 0.44 4.25 8.60 4.00 0.38 0.09	40 6	0.15 0.51 2.82 5.40	6.00 0.59 0.10
	20 8	0.42 0.60 1.68 3.00 8.00 0.56 0.07	20 10	0.53 0.62 1.47 2.00 10.00 0.66 0.07
	25 8	0.33 0.59 1.91 3.20 8.00 0.54 0.07	25 10	0.42 0.60 1.67 3.00 10.00 0.64 0.06
	30 8	0.28 0.57 2.11 4.00 8.00 0.53 0.07	30 10	0.34 0.59 1.84 3.00 10.00 0.70 0.07
	35 8	0.24 0.56 2.27 4.40 8.00 0.61 0.08	35 10	0.29 0.60 2.01 3.60 10.00 0.62 0.06
	40 8	0.21 0.57 2.41 4.00 8.00 0.60 0.08	40 10	0.26 0.59 2.13 4.00 10.00 0.66 0.07

Table 6 :

 6 Average results for the unweighted SSP over sets of five instances of small world graphs for parameter k = 4 and k = 10 and number of nodes n. We display results without improvements (M SS ), with the help of symmetry-breaking constraints (M SS +(40)) and with the help of the variable reduction technique (M SS + ∆ SS ). For each graph size, we display the average running time of the model, the average relative gap and the number of instances out of five that are solved within the one hour time limit. The time is displayed in seconds.We display in Tables 7 and 8 the results of model M SS with symmetry breaking constraints (40) and the variable reduction preprocessing on small world graphs, similar to Section 4. The SSP is still very hard to solve for such instances, even for very small graphs. As for random graphs, we can see that instances with a larger k parameter (and hence with a larger edge density) are usually easier to solve. Contrary to Section 4, however, no instance with 40 nodes is ever solved to optimality since the density for such instances do not even reach density 0.3 or 0.4, even with k = 10.

		s) gap (%) # solved	time (s) gap (%) # solved	time (s) gap (%) # solved
	20 4	239.2	0.00	5	140.2	0.00	5	103.6	0.00	5
	25 4	2479.6	18.04	3	2153.2	19.24	3	1936.2	0.00	5
	30 4	3600.0	40.96	0	3600.0	32.73	0	3600.0	44.01	0
	35 4	3600.0	79.16	0	3600.0	62.24	0	3600.0	79.46	0
	40 4	3600.0	83.82	0	3600.0	79.15	0	3600.0	71.42	0
	20 10	27.4	0.00	5	24.2	0.00	5	22.2	0.00	5
	25 10	138.2	0.00	5	110.4	0.00	5	88.2	0.00	5
	30 10	461.2	0.00	5	838.0	0.00	5	373.4	0.00	5
	35 10	2763.8	20.99	2	3264.2	12.03	2	2402.4	11.34	3
	40 10	3600.0	57.59	0	3600.0	68.88	0	3600.0	55.99	0

Table 7 :

 7 Average results of model M SS with symmetry breaking constraints (40) and the variable reduction preprocessing for the SSP, over sets of five instances for small world graphs, for each average number of neighbours k and graph size n. For both unweighted and weighted instances, we display the average running time of the model, the average relative gap and the number of instances out of five that are solved within the one hour time limit. The time is displayed in seconds.

			Unweighted			Weighted
	n k	time (s) gap (%) # solved	time (s) gap (%) # solved
	20 4	44.4	0.00	5	29.4	0.00	5
	25 4	457.6	0.00	5	410.2	0.00	5
	30 4	952.2	0.00	5	1060.0	0.00	5
	35 4	3526.6	21.70	2	3552.0	38.68	1
	40 4	3600.0	47.06	0	3600.0	39.16	0
	20 6	11.8	0.00	5	9.2	0.00	5
	25 6	72.8	0.00	5	68.2	0.00	5
	30 6	1785.8	6.00	4	1397.0	10.49	4
	35 6	3600.0	40.59	0	3180.6	15.97	3
	40 6	3600.0	49.11	0	3600.0	53.85	0
	20 8	10.2	0.00	5	10.0	0.00	5
	25 8	60.6	0.00	5	61.4	0.00	5
	30 8	283.2	0.00	5	301.2	0.00	5
	35 8	2121.6	13.57	3	1953.8	9.29	3
	40 8	3600.0	46.95	0	3600.0	39.27	0
	20 10	16.8	0.00	5	12.0	0.00	5
	25 10	76.4	0.00	5	74.8	0.00	5
	30 10	237.6	0.00	5	203.4	0.00	5
	35 10	787.0	0.00	5	2102.2	6.65	3
	40 10	3600.0	39.81	0	2922.8	12.01	3

Table 8 :

 8 Results of model M SS with symmetry breaking constraints (

  The results confirm that model M SS (with MILP improvements) is able to solve instances of up to 30 nodes without problem.

	n		δ e = 0.3			δ e = 0.5			δ e = 0.7
		M SS	MILP of [16]	M SS	MILP of [16]	M SS	MILP of [16]
		UB time UB	time	UB time UB	time	UB time UB	time
	20	8	11	8	232	9	15	9	223	10	22	10	194
	21	8	16	8	298	10	27	10	409	10	37	10	190
	22	8	18	8	379	10	15	10	715	11	27	11	750
	23	9	34	9	614	11	41	11	467	11	94	11	546
	24	10	38	10	381	11	58	11	658	12	110	12	2434
	25	11	59	11	1091	12	44	12	908	12	123	12	1768
	26	11	78	11	1561	13	97	13	3307	13	148	13	2389
	27	11	150	11	1854	13	134	13	2462	13	168	13	3899
	28	12	284	12	4557	13	248	13	3891	14	194	14	2709
	29	12	182	12	4277	14	130	14	4089	15	178	15	4680
	30	13	472	13	4963	15	359	15	7200	15	332	15	7200

Such a value for the number of connected components of G[V \ S * ] can be obtained when G is a star graph.

Acknowledgements

I would like to thank Rosario Scatamacchia for sharing C++ code useful for generating the benchmark instances and for fruitful discussions, as well as the constructive and detailed comments of three anonymous referees which greatly helped enhance the manuscript.

Appendix A Results over small world networks

In this appendix we present detailed results over small networks generated through the mechanism described in [START_REF] Watts | Collective dynamics of small-world networks[END_REF]. These networks typically have a small diameter and represent social networks fairly well. They are known to be notoriously difficult to solve for some partitioning problems, such as [START_REF] Aringhieri | Local Search Metaheuristics for the Critical Node Problem[END_REF][START_REF] Zhou | Memetic search for identifying critical nodes in sparse graphs[END_REF]. We generated instances with a number of nodes in the set n ∈ {20, 25, 30, 35, 40} as in Section 4. The parameters needed to generate small world graphs are the average number of neighbours k and a rewiring probability of edges which we choose as p = 0.05, similar to [START_REF] Ventresca | Global search algorithms using a combinatorial unranking-based problem representation for the critical node detection problem[END_REF]. We vary parameter k ∈ {4, 6, 8, 10} to obtain instances with different characteristics and generate five instances for each value of n and k. The different average coefficients over those sets of five instances are displayed in Table 5 in a manner similar to Table 1, except for the additional column δe which represents the average density of the instances. We can clearly see that the clustering coefficient is usually much larger for small world graphs while the ratio of the dispersion (standard deviation) in the nodes degree over the average degree (last column) is much smaller, denoting the fact that two different nodes tend to have similar degree in a small world graph.

Appendix B Comparison with the model of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] We compare in the following results of our model M SS with those of the MILP approach in [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] on a set of benchmark instances they introduced. As discussed in Section 1, the algorithm advocated in [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] is a branch-and-cut algorithm based on a model with an exponential number of constraints, implemented as lazy constraints. The results of both models are reported in Table 9 for unweighted instances of SSP and in Table 10 for weighted instances. The instances range from 10 to 30 nodes with densities 0.3, 0.5 and 0.7. We report the upper bound (UB) and running time of both models. For the model of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF], when the running time is 7200 seconds, the algorithm reached its time limit and the instance is not solved at optimality, n δ e = 0.3 δ e = 0. 10: UB and running time for the benchmark instances of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF] for weighted instances with densities 0.3, 0.5 and 0.7, for model M SS with improvements and the branch-and-cut algorithm of [START_REF] Macambira | A new formulation for the safe set problem on graphs[END_REF].