During these three years of thesis, two papers related to graphons were published, which partially originated in the internship corresponding to the Master Degree in Systems, Control and Information Technologies at the Université Grenoble Alpes. Their objective is to study the properties of graphons that can be used for the analysis of multi-agent systems where the network topology is sampled from a graphon. Although some results are used in this thesis, specially in Section 2.

In this thesis we address several problems formulated in an open multi-agent system (OMAS) scenario where the set of agents can change in time, independently of the evolution of the dynamics associated with the system. We analyze the case of OMAS formulated in a fixed finite size network and we use two approaches for the analysis of the systems. In the first approach we consider scenarios characterized by activation/deactivation of agents such that at each time instant a different subset of active agents can interact in the system. In the second approach we study scenarios characterized by replacements of agents where at a specific time instant, an agent can be replaced while the rest of the agents remain the same. In this case, all the agents are able to interact at all time.

Three different problems are considered in this thesis: randomized consensus, resource allocation problem and epidemics. First, we analyze the problem of randomized consensus subject to additive noise where different subset of agents exchange information at each iteration. We define a noise index based on the expected mean squared error and we derive upper bounds. Then, we consider the resource allocation problem where agents can be replaced during the implementation of an optimization algorithm. For this optimization problem, we analyze two different algorithms: weighted gradient descent and random coordinate descent. For the weighted gradient descent, we evaluate the performance of the algorithm in an OMAS subject to packet losses by defining appropriate performance metrics. For the random coordinate descent algorithm, we study the convergence to the minimizer in an OMAS and we provide an alternative analysis using tools inspired from online optimization. Finally, we study a SIS epidemic in continuous time subject to replacements of agents during the evolution of the disease. We perform the analysis using an aggregate function to condensate the information of the system and we derive upper bounds for its asymptotic behavior. i

Aperçu de la thèse

Systèmes multi-agents ouverts (OMAS)

Les systèmes multi-agents sont composés d'entités individuelles interagissant entre elles dans un environnement commun pour accomplir une tâche collective ou atteindre un objectif commun. Les interactions sont généralement définies à l'aide d'un graphe composé d'un ensemble de noeuds représentant les agents et d'un ensemble d'arêtes modélisant les échanges d'information entre les agents. En général, les chercheurs considèrent un ensemble fixe d'agents pour la dérivation des résultats, tandis que l'ensemble d'arêtes peut varier dans le temps. Cependant, dans de nombreux scénarios, l'ensemble d'agents peut également changer puisque des agents peuvent rejoindre ou quitter le système ou simplement être remplacés. Un système multi-agents soumis à des ajouts, des suppressions ou des remplacements d'agents est appelé un système multi-agents ouvert (OMAS). Si ces changements potentiels dans l'ensemble des agents ne sont pas fréquents, l'étude des systèmes multi-agents peut être effectuée sous l'hypothèse que les agents restent les mêmes dans le temps, même si les connexions peuvent varier. Cependant, si des ajouts, des suppressions ou des remplacements se produisent à une échelle de temps similaire à celle du système dynamique sur le réseau, il est nécessaire de considérer des OMAS pour l'analyser. Parmi les applications où l'usage des OMAS peut être envisagée, on trouve par exemple les réseaux sociaux, les épidémies, les convois de véhicules, les smart grids, la robotique en essaim, etc.

Les systèmes multi-agents soumis à des remplacements sont un cas particulier d'OMAS, où le départ d'un agent est immédiatement suivi de l'arrivée d'un nouvel agent de sorte que la taille du système ne change pas. Lorsque l'intervalle de temps entre un départ et une arrivée est négligeable, l'occurrence de ces deux événements peut être considérée comme simultanée, et l'analyse peut être effectuée en ne considérant que les remplacements. Ce comportement particulier peut également être vu comme une approximation de l'OMAS où les taux d'arrivées et de départs sont similaires.

iii

Principaux défis dans les OMAS

Les changements possibles dans la cardinalité de l'ensemble des agents augmentent considérablement la complexité d'OMAS et posent de nombreux défis pour l'analyse de ces systèmes.

Espaces d'états et trajectoires : lorsque l'ensemble des agents change de cardinalité, les agents peuvent appartenir à différents espaces d'états et la définition des trajectoires ne peut pas être facilement formulée.

Modèles d'événements : au cours de l'évolution du système, les OMAS se caractérisent par la survenance de trois types d'événements : les arrivées, les départs et les remplacements. Une identification parfaite des instants temporels auxquels ces événements se produisent ne semble pas plausible et des hypothèses appropriées doivent être considérées pour modéliser des OMAS et garder l'analyse traitable.

Topologies de graphes : il existe quelques cas où les topologies de graphes suivent un modèle spécifique qui peut être complètement étendu pour n'importe quel nombre d'agents. Cependant, il est peu probable que les réseaux réels suivent des topologies de graphes restrictives et conservent également les mêmes connexions à tout moment.

Formulation en temps : la plupart des travaux antérieurs sur les OMAS considèrent une formulation en temps discret, où les agents n'interagissent qu'à des instants de temps précis. Cependant, de nombreux phénomènes complexes, notamment associés aux systèmes chaotiques, ne peuvent pas être reproduits par les systèmes en temps discret et des formulations en temps continu sont nécessaires.

Approches potentielles de l'OMAS

Après avoir identifié les enjeux des OMAS, nous passons maintenant en revue les principales approches qui pourraient être utilisées dans leur étude.

Sur-ensemble fini invariant dans le temps : des OMAS peuvent être sousapproximés en considérant un sur-ensemble fini d'agents où l'ouverture du système est caractérisée par la activation et la désactivation des agents, de sorte qu'à chaque instant un nombre différent d'agents peut interagir dans le système. Cette approche sera considerée dans le Chapitre 2.

Espace d'état fixe : lorsqu'un OMAS est soumis uniquement à des remplacements, la dimension du système ne change pas avec le temps, et lorsqu'un agent quitte le réseau, un nouvel agent rejoint immédiatement le système en gardant le même label. Cette approche sera poursuivie dans les Chapitres 3, 4 et 5.

Systèmes multidimensionnels multimodes : étant donné que la taille du système peut changer en raison de l'occurrence d'événements liés aux OMAS, des systèmes de commutation entre des espaces de différentes dimensions peuvent être utilisés pour l'analyse.

Systèmes dynamiques de dimension infinie : le système dynamique de dimension infinie peut être utilisé comme sur-approximation de OMAS. L'objectif est de résoudre le problème en dimension infinie puis d'utiliser des techniques d'approximation pour obtenir la solution au cas d'OMAS avec une certaine marge d'approximation.

Principales contributions de cette thèse

L'intérêt principal de cette thèse est l'analyse d'OMAS en utilisant deux approches différentes : le sur-ensemble fini invariant dans le temps et l'espace d'état fixe. Nous considérons trois problèmes différents : l'influence du bruit dans le consensus randomisé, le problème d'allocation des ressources et la propagation d'une maladie. Nous formulons ces problèmes en utilisant une approche stochastique et nous donnons des bornes pour des quantités scalaires qui sont définies dans le but d'évaluer les impacts des activations/désactivations ou des remplacements d'agents dans le système.

Dans le Chapitre 2 nous analysons le problème du consensus randomisé perturbé par un bruit additif. Nous considérons la formulation en temps discret d'un algorithme de consensus randomisé linéaire où les matrices d'interaction sont symétriques et tirées d'une distribution commune. Ensuite, nous analysons le cas où le bruit additif modifie les états des agents et nous dérivons une expression sous forme fermée de l'erreur quadratique moyenne induite par le bruit, ainsi que des bornes supérieure et inférieure plus simples à évaluer. Ensuite, nous considérons le cas d'OMAS où les matrices à chaque instant correspondent au sous-ensemble d'agents actifs dans le système. L'activation des agents est déterminée par des variables aléatoires de Bernoulli indépendantes associées à chaque agent potentiel du système. Pour les OMAS, nous exprimons les bornes en utilisant les valeurs propres de la matrice laplacienne du graphe sous-jacent ou la résistance effective moyenne du graphe, prouvant ainsi leur qualité. Enfin, nous dérivons des expressions pour les bornes sur quelques exemples de graphes et les évaluons numériquement. v Le Chapitre 3 se concentre sur l'analyse de l'algorithme Weighted Gradient Descent pour résoudre le problème d'allocation de ressources (RA) sous pertes de paquets. Nous considérons le cas particulier des fonctions de coût quadratiques par morceaux et définissons deux métriques de performance qui mesurent, respectivement, l'écart à la contrainte et l'erreur sur la fonction de coût attendue. Nous dérivons des bornes supérieures sur les deux métriques qui sont proportionnelles à la différence entre la fonction de coût initiale et la fonction de coût évaluée au minimiseur. Ensuite, nous étendons l'analyse de la violation de contrainte aux OMAS où les agents peuvent être remplacés, ce qui implique un changement sur la fonction de coût associée. Nous montrons que la combinaison des remplacements et des pertes fait diverger l'erreur de violation de contrainte avec le temps.

Le Chapitre 4 est consacré à l'analyse de l'algorithme Random Coordinate Descent (RCD) pour résoudre le problème RA dans des OMAS sous l'hypothèse que les fonctions locales sont lisses fortement convexes et ont leurs minimiseurs situés dans une boule donnée. Nous établissons la convergence linéaire de l'algorithme dans les systèmes fermés, en termes d'estimation vers le minimiseur, pour des graphes généraux et des pas appropriés et nous estimons le changement de la solution optimale après un remplacement pour évaluer son effet sur la distance entre les estimation et le minimiseur. A partir de ces deux éléments, nous dérivons des conditions de stabilité dans des OMAS où les agents peuvent être remplacés et établissons la convergence linéaire de l'algorithme vers une erreur en espérance en régime permanent. De plus, en considérant des paramètres plus simples, nous analysons l'algorithme RCD dans des OMAS à l'aide d'outils similaires à ceux couramment utilisés dans l'optimisation en ligne. En particulier, nous comparons la solution issue de l'algorithme RCD avec la solution optimale et la stratégie égoïste. Ce chapitre a été développé en collaboration avec Charles Monnoyer de Galland et Julien Hendrickx de l'Université Catholique de Louvain (UCLouvain).

Dans le Chapitre 5 nous fournissons une analyse préliminaire de la propagation d'une maladie en temps continu dans des OMAS soumis à des remplacements. Nous considérons le modèle épidémique SIS sur un graphe et définissons une fonction agrégée pour l'analyse du système. Premièrement, nous considérons le cas des remplacements sur un graphe donné fixe où l'occurrence des remplacements est déterminée par un processus de Poisson. Ensuite, nous étudions le cas de graphes échantillonnés à partir de graphons où nous dérivons des bornes supérieures qui dépendent des propriétés du graphon. Enfin, nous fournissons une analyse préliminaire de scénarios où les connexions changent également lors des remplacements selon un graphon, et nous formulons une conjecture pour le comportement asymptotique de la fonction d'agrégation. vi First, I would like to thank God for giving me the strength to keep going.

I would like to express my sincere gratitude to my supervisors Elena Panteley and Paolo Frasca for guiding me along these three years of research and for their genuine dedication at all levels. I learned a lot from our interactions and thanks to you I always found the motivation to continue with the research despite the difficulties. Thank you for your advice, help and support.

I would like to express my appreciation to Profs. Giacomo Como, Giuseppe Notarstefano and Franck Iutzeler for your reviews that helped me improving this manuscript. Thanks also to Antoine Girard and Federica Garin for accepting being members of the jury for my defense.

My thanks also go to Prof. Julien Hendrickx for giving me the opportunity of visiting his research lab in the Université Catholique de Louvain. It was a wonderful experience to spend those weeks with him and his fantastic research team.

I would like to thank to my colleague and friend Charles Monnoyer de Galland, with whom I had the opportunity of working together for a long period of time. I really enjoyed our exchanges on math, research and life, together with Iona. Thank you for your guidance during my stay in Louvain-la-Neuve.

I thank all my fellows labmates and friends in GIPSA-lab, whom I had the chance to meet during the years I stayed in the lab. I really enjoyed our discussions and jokes during the lunch and the different activities and games that we had together. Thank you for making the lab a pleasant environment to work in, especially during the epidemic. I wish I had been with you until the end of my doctorate.

Let me also express my thanks to my colleagues and friends in L2S, with whom I spent the last months of my doctorate having pleasant conversations during the lunch.

Last but not least, I would like to thank to my family. My parents, Roberto and Dolores, who always gave me their love and support to overcome all the difficulties in my life. I also would like to thank my cousin Pablo and his wife Margarita, my family in Europe, who always welcomed me warmly during my visits to Madrid List of symbols and acronyms In the last years, an increasing number of real-world problems are being understood as constantly evolving towards a networked environment where single entities interact with each other to accomplish a collective task or achieve a common objective. Systems of this type are called multi-agent systems and have been extensively studied in the scientific community considering different types of dynamics, coupling, goals and time formulation [START_REF] Bullo | Lectures on Network Systems[END_REF].

Regarding the dynamics, each agent is in itself a dynamical system that can be modeled by ordinary differential equations (ODE). At the most basic level, each agent i is commonly considered as a single integrator of the form:

ẋi (t) = u i (t),
where u i (t) is the input of the agent. However, this type of models does not reflect the behavior of more complex systems. For this reason, researchers have considered the use of more complex dynamics to model each single agent. In this line, models such as double integrators, Lagrangian systems, etc. [START_REF] Ren | Distributed Consensus in Multi-vehicle Cooperative Control[END_REF][START_REF] Ren | Distributed Coordination of Multi-agent Networks: Emergent Problems, Models, and Issues[END_REF], have been studied with the objective of extending well-known results to more complex dynamics.

With respect to the goals, researchers have considered several problems associated with multi-agent systems including: stability, consensus, synchronization and formation control. Stability corresponds to the analysis of the equilibrium points of the system. Consensus is a special case of coordination of multi-agent systems where the objective is to reach a common value for all the agents which must be static [START_REF] Olfati-Saber | Consensus and cooperation in networked multi-agent systems[END_REF]. As a generalization of consensus, many researchers also consider the problem of synchronization where agents must reach a common value, but this value is time-varying. The synchronization problem is often considered in networks of oscillators, where the trajectories of agents are more chaotic [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF]. In the case of formation control, agents must reach specific locations determined by desired geometric patterns [START_REF] Oh | A survey of multi-agent formation control[END_REF].

Furthermore, in the last years, the community has paid attention to the extension of optimization problems to the case of multi-agent system. In addition to its states, an agent has a local cost function depending only on the state of the agent, and the objective is the minimization of the sum of the cost functions of all the agents while agents exchange information only with the neighbors1 . Additional constraints can be added depending on the nature of the problem. Generally, most of the algorithms make use of the well known subgradients of the cost fuctions to solve the optimization problem [START_REF] Nedic | Distributed subgradient methods for multi-agent optimization[END_REF]. However, other methods like primal duals have been adapted to multi-agent systems [START_REF] Notarstefano | Distributed Optimization for Smart Cyber-Physical Networks[END_REF].

While the links among the agents can be modified by nonlinear functions as in the case of oscillators [START_REF] Dörfler | Synchronization in complex networks of phase oscillators: A survey[END_REF], the connections in most of the multi-agent systems are characterized by diffusive couplings. For this particular case, the network can be represented by a matrix and tools of Linear Algebra can be used for the analysis of multi-agent systems. Furthermore, due to the link between specific type of matrices and graphs, Spectral Graph Theory [START_REF] Cvetković | An Introduction to the Theory of Graph Spectra[END_REF] plays an important role in the study of interactions of agents over a network .

In many scenarios, the links in multi-agent systems do not remain static and change with time. For this reason, researchers focused on the analysis of classical problems applied to time-varying graphs, where under several assumptions on the connectivity of the possible graph topologies in a time interval, many results were easily extended [START_REF] Moreau | Stability of continuous-time distributed consensus algorithms[END_REF][START_REF]Stability of multiagent systems with time-dependent communication links[END_REF][START_REF] Lin | State agreement for continuous-time coupled nonlinear systems[END_REF]. Furthermore, in many cases, a complete knowledge of a sequence of graph topologies is not available and the attention was focused on randomized interactions, where the graph topologies are sampled from a common distribution and tools from Probability Theory can be used [START_REF] Hatano | Agreement over random networks[END_REF][START_REF] Fagnani | Randomized consensus algorithms over large scale networks[END_REF][START_REF] Tahbaz-Salehi | A necessary and sufficient condition for consensus over random networks[END_REF][START_REF] Frasca | On the mean square error of randomized averaging algorithms[END_REF]. In this way, the complexity of multi-agent systems regarding the network topology were explored considering only possible changes in the current set of edges while the set of agents was assumed to be always static.

Open multi-agent systems (OMAS)

Even if models based on time-varying switching topologies provide a good approach for studying the behavior of systems where communications among the agents do not remain static, they do not encompass possible changes in the set of agents that could significantly affect the evolution of several systems. In the real world, agents can join and leave the network, and therefore the size of many multi-agent systems is changing with time. A multi-agent system subject to additions, removals or replacements of agents is called an open multi-agent system (OMAS) [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF][START_REF] Franceschelli | Proportional dynamic consensus in open multi-agent systems[END_REF]. If these potential changes in the set of agents are not frequent, the study of the multi-agent system can be performed under the assumption that agents remain the same in time, even if the connections may vary. However, if additions, removals or replacements occur in a time-scale similar to that of the dynamical system on the network, it is necessary to consider OMAS for the analysis of the dynamics. Among the potential applications where OMAS can be used we have: Social networks: users enter and leave communities, discussions, forums, etc., modifying the group behavior.

Epidemics: infected individuals can enter or leave towns, villages, cities, etc., affecting the spread of diseases.

Vehicle platoons: during transportation, vehicles may join or leave the platoon according to their particular destinations.

Smart-grids: smart appliances are constantly switched on and off in a network of electrical devices, modifying the number of active devices.

Swarm robotics: the size of a group of robots can change depending on the environmental conditions and goals of the team.

Historically, computer science is one of the first research fields where the concept of OMAS has been considered in the analysis of trust and security, however, no dynamics are associated with the individual agents [START_REF] Huynh | An integrated trust and reputation model for open multi-agent systems[END_REF][START_REF] Pinyol | Computational trust and reputation models for open multi-agent systems: a review[END_REF]. Also, in the field of game theory, researchers have analyzed the problem of price of anarchy bounds in dynamic populations where players can exit and be replaced by new participants [START_REF] Lykouris | Learning and efficiency in games with dynamic population[END_REF]. Regarding dynamical systems, several works have indirectly considered potential cases of addition/removal of agents as an adverse scenario to analyze the robustness of algorithms [START_REF] Zhu | Discrete-time dynamic average consensus[END_REF][START_REF] Patel | Dynamic partitioning and coverage control with asynchronous one-to-base-station communication[END_REF].

Preliminary works focused on consensus in OMAS consider the analysis of the pairwise gossiping algorithm in a complete graph [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF][START_REF] Hendrickx | Open multi-agent systems: Gossiping with random arrivals and departures[END_REF], where a discrete time formulation was used to analyze replacements and arrivals of agents determined by deterministic and random time instants. Furthermore, in [START_REF] Monnoyer De Galland | Modelling gossip interactions in open multi-agent systems[END_REF], authors extended the results to the case of arrivals/departures when the time instants are determined by Poisson processes. In these works, the variance of the system was the object of study, which condensates the information of the system in a single value like an aggregate function. Following this line of research, in [START_REF] Monnoyer De Galland | Fundamental performance limitations for average consensus in open multi-agent systems[END_REF], the authors provided lower bounds on the mean squared error for any averaging algorithm applied to an OMAS. Another work based on gossip interactions is [START_REF] Abdelrahim | Max-consensus in open multi-agent systems with gossip interactions[END_REF], where the MAXconsensus in an OMAS over a complete graph is analyzed.

Regarding deterministic formulations of OMAS, one of the preliminary works is [START_REF] Franceschelli | Stability of open multiagent systems and applications to dynamic consensus[END_REF], where authors study the dynamic consensus in a contractive OMAS under the assumption that the arrival process is bounded. Other works related to the analysis of the dynamic consensus in OMAS are [START_REF] Sanai Dashti | Dynamic consensus on the median value in open multi-agent systems[END_REF] and [START_REF] Deplano | Dynamic min and max consensus and size estimation of anonymous multi-agent networks[END_REF] where the median value and the max value are the object of the study respectively and the assumption of a dwell time is used for the derivation of the results. In [START_REF] Xue | Stability of multi-dimensional switched systems with an application to open multi-agent systems[END_REF], the authors analyze the stability of multi-dimensional switched systems, where restrictions are imposed on the switching signal to guarantee practical consensus.

Other works include the analysis of OMAS under stochastic interactions in a finite superset of agents where the activation of agents are determined by Bernoulli random variables [START_REF] Varma | Open multi-agent systems with discrete states and stochastic interactions[END_REF]. Using a similar approach, an algorithm based on dual averaging was proposed to minimize a global cost function that depends on a time-varying set of active agents in a fixed size network [START_REF] Hsieh | Optimization in open networks via dual averaging[END_REF]. The study of the Friedkin-Johnsen models where agents exchange roles at different time instant may also be encompassed in OMAS since the set of agents can be seen as time-varying due to the changes in roles [START_REF] Wang | Achieving consensus in spite of stubbornness: time-varying concatenated Friedkin-Johnsen models[END_REF].

The case of replacements

Multi-agent systems subject to replacements is a special case of OMAS, where the departure of an agent is immediately followed by the arrival of a new agent such that the size of the system does not change. When the time interval between a departure and an arrival is negligible, the occurrence of these two events can be considered as simultaneous, and the analysis can be performed considering only replacements. This particular behavior can also be seen as an approximation of OMAS where the rate of arrivals and departures are similar. The particular case of replacements has been studied in opinion dynamics, where researchers analyzed scenarios when the states of agents are reset [START_REF] Grauwin | Opinion group formation and dynamics: Structures that last from nonlasting entities[END_REF]. In the analysis of consensus in OMAS, replacements have been studied in [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF][START_REF] Hendrickx | Open multi-agent systems: Gossiping with random arrivals and departures[END_REF], where the authors considered the variance for the analysis of agents interacting in a complete graph. Regarding replacements in optimization, in [START_REF] Hendrickx | Stability of decentralized gradient descent in open multi-agent systems[END_REF], the authors analyzed the gradient descent algorithm in a network subject to replacements of the cost functions of the agents.

Main challenges in OMAS

Possible changes in the cardinality of the set of agents increase considerably the complexity of OMAS and pose many challenges for the analysis of these systems.

State spaces and trajectories: the concept of state spaces and trajectories are fundamental in the study of classical multi-agent systems where it is assumed that the agents always remain in the same state-space and they have an associated trajectory for all time [START_REF] Bullo | Lectures on Network Systems[END_REF]. However, when the set of agents is changing in cardinality, agents can belong to different state-spaces and the definition of trajectories cannot be easily formulated since the states of the agents may not be well-defined for some time intervals. For instance, during the removal of an agent, its associated trajectory is completely removed while when an agent joins the system, a new trajectory is created. Even if only replacements occur in the system, trajectories experiment jumps, which can drastically change the behavior of the current state.

Models of events: during the evolution of the system, OMAS are characterized by the occurrence of three types of events: arrivals, departures and replacements [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF][START_REF] Franceschelli | Stability of open multiagent systems and applications to dynamic consensus[END_REF]. A perfect identification of the time instants at which these events take place do not seem plausible in a real application, and appropriate assumptions must be considered to model OMAS and keep the analysis tractable. In this line, while the use of dwell times can guarantee some lower bounds on the time interval between the occurrence of two distinct events, the use of stochastic processes could model richer behaviors of OMAS that may be closer to real applications.

Network topologies: there are few specific cases where the graph topologies follow a specific pattern that can be completely extended for any number of agents such as: line, cycle, grids, complete or bipartite graphs [START_REF] Fagnani | Introduction to Averaging Dynamics over Networks[END_REF]. However, it is unlikely that real networks will follow graph topologies that are restrictive, and also keep the same connections at all time. For this reason, new graph topologies must be explored, which could include random graphs [START_REF] Bollobás | Random Graphs[END_REF] and graphons [START_REF] Lovász | Large Networks and Graph Limits[END_REF].

Time formulation: unlike few works in continuous time such as [START_REF] Sanai Dashti | Dynamic consensus on the median value in open multi-agent systems[END_REF][START_REF] Xue | Stability of multi-dimensional switched systems with an application to open multi-agent systems[END_REF], most of the previous works in OMAS consider a discrete time formulation, where agents only interact at discrete time instants. Even if a discrete time formulation can in many cases reproduce the behaviors of classical dynamical systems formulated as ordinary differential equations (ODE), many complex phenomena, specially associated with chaotic systems, cannot be modeled with discrete time systems [START_REF] Krivine | Discrete-time and continuous-time modelling: some bridges and gaps[END_REF]. For this reason, a continuous time formulation of OMAS is necessary to study other types of dynamics, which implies the study of more complex processes, where the state of the system also changes between arrivals, departures and replacements.

Potential approaches to OMAS: pros and cons

After having identified the challenges in OMAS, we now proceed to review the main approaches that could be used in their study.

Time-invariant finite superset: OMAS can be under-approximated by considering a finite superset where it is assumed that we have the information about all the potential agents [START_REF] Varma | Open multi-agent systems with discrete states and stochastic interactions[END_REF][START_REF] Hsieh | Optimization in open networks via dual averaging[END_REF]. The openness of the system is characterized by the activation and deactivation of agents, such that at each time instant different number of agents may interact in the system. Since the dimension of the global system is constant, standard tools from dynamical systems can be used for the analysis. However, this approach is restricted to systems where a full knowledge of all the possible agents is available such that unexpected agents cannot join the network. This approach will be pursued in Chapter 2. Fixed state-space: when an OMAS is subject only to replacements, the dimension of the system does not change with time, and the state of the system always remains in the same state-space. Unlike the time-invariant finite superset, all the agents interact at the same time instants, but when an agent leaves the network, a new agent will join the system immediately keeping the same label. Trajectories are well-defined but they may experience jumps during the replacements. Since the state-space does not change, standard concepts of dynamical systems can be applied. However, this approach is restricted to replacements and cannot be used in the case of decoupled arrivals and departures. This approach will be pursued in Chapters 3, 4 and 5.

Multi-mode multi-dimensional systems: since the size of the system may change due to the occurrence of arrivals or departures, switching systems between spaces of different dimensions may be used for the analysis. Preliminary works in this field correspond to [START_REF] Verriest | Multi-mode multi-dimensional systems[END_REF][START_REF]Pseudo-continuous multi-dimensional multi-mode systems[END_REF], where the author studied the switching between a fixed number of spaces of different dimensions. The family of possible vector spaces to which the state belongs is defined as a discrete vector bundle and concepts of Topology are required for the analysis. For this framework, the concept of trajectories is meaningful and a global equilibrium for the system can be defined, allowing the use of some tools of hybrid systems [START_REF] Goebel | Hybrid dynamical systems[END_REF]. However, the classical concept of stability cannot be applied and the approach can be rigid for many applications since it is required to know all the possible vector spaces to define the discrete vector bundle. Infinite dimensional dynamical systems: infinite dimensional dynamical system may be used as an over-approximation of OMAS. In this case, the problem can be formulated as a dynamical system in a space of sequences when we consider a infinite network [START_REF] Noroozi | Small-gain theorem for stability, cooperative control and distributed observation of infinite networks[END_REF] or in a space of mesurable functions by using a graphon [START_REF] Gao | Graphon control of large-scale networks of linear systems[END_REF]. The objective is to solve the problem in infinite dimensions and then use approximation techniques to obtain the solution to the case of OMAS with a certain margin of approximation. The main advantages of this approach is that all the possible agents can be considered in the analysis of the system and advanced tools from infinite dimensional systems can be used. However, the definition of appropriate mappings from the finite network to the infinite system and vice versa could make it difficult to use this approach.

Main contributions of this thesis

The main interest of this work is the analysis of OMAS using two different approaches: time-invariant finite superset and fixed state-space. We consider three different problems: influence of noise in randomized consensus, resource allocation problem and the spread of a disease. We formulate the problems using a stochastic approach and we provide bounds for scalar quantities that are defined with the objective of evaluating the impacts of activations/deactivations or replacements of agents in the system.

In Chapter 2 we analyze the problem of randomized consensus perturbed by additive noise. We consider the formulation in discrete time of a linear randomized consensus algorithm where the interaction matrices are symmetric and drawn from a common distribution. Then, we analyze the case where additive noise modifies the states of the agents and we derive a closed form expression for the mean squared error induced by the noise, together with upper and lower bounds that are simpler to evaluate. Next, we consider the case of OMAS where the interaction matrices at each time instant are generated by the subset of active agents in the system. The activation of agents are determined by independent Bernoulli random variables associated to each potential agent in the system. For OMAS, we express the bounds by using the eigenvalues of the Laplacian matrix of the underlying graph and the graph's average effective resistance, thereby proving their tightness. Finally, we derive expressions for the bounds for particular graph topologies and numerically evaluate them.

Chapter 3 focuses on the analysis of the weighted gradient descent algorithm to solve the resource allocation (RA) problem under packet losses. We consider the particular case of piecewise quadratic cost functions and define two performance metrics that measure, respectively, the deviation from the constraint and the error on the expected cost function. We derive upper bounds on both metrics which are proportional to the difference between the initial cost function and the cost function evaluated at the minimizer. Then, we extend the analysis of the constraint violation to OMAS where agents can be replaced, which implies a change on the associated cost function. We show that the combination of replacements and losses makes the constraint violation error diverge with time.

Chapter 4 is devoted to the analysis of the Random Coordinate Descent (RCD) algorithm to solve the RA problem in OMAS under the assumptions that local functions are smooth strongly convex and have their minimizers located in a given ball. We establish the linear convergence of the algorithm in closed systems, in terms of the estimate towards the minimizer, for general graphs and appropriate step-size and we estimate the change of the optimal solution after a replacement to evaluate its effect on the distance between the estimate and the minimizer. From these two elements, we derive stability conditions in OMAS where agents can be replaced and establish the linear convergence of the algorithm towards a steady state expected error. Additionally, by considering a simple setting, we analyze the RCD algorithm in OMAS using tools similar to those commonly used in online optimization. In particular, we study the accumulated errors that compare solutions issued from the RCD algorithm and the optimal solution or the non-collaborating selfish strategy and we derive some bounds in expectation for these accumulated errors. This chapter was developed in collaboration with Charles Monnoyer de Galland and Julien Hendrickx of the Université Catholique de Louvain (UCLouvain).

In Chapter 5 we provide a preliminary analysis of the spread of a disease in continuous time in OMAS subject to replacements. We study the SIS epidemic model over a network and inspired by the seminal work in [START_REF] Hendrickx | Open multi-agent systems: Gossiping with deterministic arrivals and departures[END_REF], we define an aggregate function for the analysis of the system based on a macroscopic description. First, we consider the case of replacements over a fixed given graph where the occurrence of replacements follows a homogeneous Poisson process and the new values of the replaced agents are determined by a random variable. We analyze the expectation and variance of the aggregate function and we provide some upper bounds. Then, we study the case of graphs sampled from graphons where we derive upper bounds which depend on properties of the graphon. Finally, we provide a preliminary analysis of scenarios where connections also change during the replacements according to a graphon, and we formulate a conjecture for the asymptotic behavior of the aggregate function.

Thesis outline. Chapter 2 analyzes the influence of noise in randomized consensus with an application to OMAS considering a finite superset of agents. Chapters 3 and 4 focus on the resource allocation problem in OMAS where agents can be replaced at each time instant. While Chapter 3 analyzes the weighted gradient descent algorithm under packet losses, Chapter 4 studies the stability of the RCD algorithm and considers the analysis in a particular setting using tools inspired by online optimization. Chapter 5 provides a preliminary analysis of a SIS epidemic in continuous time in an OMAS subject to replacements determined by a Poisson process. In Chapter 6 we summarize the results of the thesis and propose directions for future research. Appendix A introduces some mathematical concepts and results related to graphs, graphons and Poisson processes.

List of publications

The following is an exhaustive list of publications written during these three years of thesis, that are either published, under review or still in preparation. It also contains the publications related to graphons, whose objective is to develop tools for the analysis of OMAS using an infinite dimensional approach.

Journal papers

J1. R.

Chapter 2

Randomized consensus with activation/deactivation of agents subject to noise disturbance

When the set of all the possible agents in the systems is available, OMAS can be analyzed using a superset. The openness of the system is then characterized by a continuous activation and deactivation of agents [START_REF] Varma | Open multi-agent systems with discrete states and stochastic interactions[END_REF], with a network topology determined by the induced subgraph. Moreover, if the time-scale is correctly chosen such that the dependence of the activation/deactivation of agents between iterations is negligible, the interaction matrices at different time instants can be considered as independent events and the behavior of the system can be analyzed using tools associated with randomized algorithms. This corresponds for example to the case of smart devices connected to a large electrical network in a town or city, where plugging/unplugging a device can be considered as an action independent of the state of other devices on the network. Randomized interactions have been extensively studied in consensus systems [START_REF] Tahbaz-Salehi | Consensus over ergodic stationary graph processes[END_REF][START_REF] Ravazzi | Ergodic randomized algorithms and dynamics over networks[END_REF], with a wide range of applications including social networks [START_REF] Acemoğlu | Opinion fluctuations and disagreement in social networks[END_REF], sensor networks [START_REF] Kar | Sensor networks with random links: Topology design for distributed consensus[END_REF], and clock synchronization [START_REF] Bolognani | A randomized linear algorithm for clock synchronization in multi-agent systems[END_REF]. Even if consensus in random networks has been analyzed for a long time, researchers have mostly considered ideal, noiseless interactions [START_REF] Fagnani | Randomized consensus algorithms over large scale networks[END_REF]. Instead, more realistic dynamical models should at least include noise. The influence of noise has indeed been investigated in several works [START_REF] Xiao | Distributed average consensus with least-mean-square deviation[END_REF][START_REF] Jadbabaie | Scaling laws for consensus protocols subject to noise[END_REF][START_REF] Yazıcıoğlu | Structural robustness to noise in consensus networks: Impact of average degrees and average distances[END_REF] regarding deterministic consensus, both in discrete-time [START_REF] Garin | A survey on distributed estimation and control applications using linear consensus algorithms[END_REF] and in continuous-time [START_REF] Patterson | Leader selection for optimal network coherence[END_REF][START_REF] Bamieh | Coherence in largescale networks: Dimension-dependent limitations of local feedback[END_REF] but, relatively little is known about the effects of noise on randomized consensus. For instance, the authors in [START_REF] Wang | Distributed averaging under constraints on information exchange: emergence of Lévy flights[END_REF] studied additive noise in random consensus over directed graphs, showing that, due to asymmetric loss of links, uncommon phenomena such as Lévy flights can arise in the evolution of the system. However, the emergence of this phenomenon is not possible for undirected graphs, provided the links are deactivated in a symmetric way.

13

In the case of OMAS, in addition to the variation of the number of active agents interacting at each time instant, the effect of noise adds an additional perturbation to the behavior of the system. It is clear that noise prevents the states of the nodes from reaching consensus. For this reason, it is important to quantify the effects of noise by a noise index, which is simply the steady-state normalized mean squared error between the states of the nodes and their average. It is well known that for symmetric matrices, this type of indexes are closely related to the spectrum of the interaction matrices, such that closed forms and bounds are generally expressed as a function of the eigenvalues [START_REF] Garin | A survey on distributed estimation and control applications using linear consensus algorithms[END_REF].

In this chapter we study the steady-state normalized mean squared error in a randomized consensus perturbed by additive noise considering agents with scalar states. However, the majority of the results can be extended to agents with higher dimensional states. In Section 2.1 we define a noise index and obtain an explicit expression for it. Since this closed form involves the calculation of an n 2 -dimensional matrix, we derive an upper bound and a lower bound that depend on the eigenvalues of relevant n-dimensional matrices. These results generalize well-known results about deterministic consensus with noise [START_REF] Garin | A survey on distributed estimation and control applications using linear consensus algorithms[END_REF][START_REF] Bamieh | Coherence in largescale networks: Dimension-dependent limitations of local feedback[END_REF]. Section 2.2 presents the analysis of randomized consensus with noise in OMAS characterized by activations and deactivations of agents. We refine our analysis for a specific class of random update matrices related to OMAS, which we refer to as Randomly Induced Discretized Laplacians (RIDL). Given an underlying large graph, these matrices are generated by sampling a subset of active nodes and considering the subgraph induced by the active nodes. When the update matrices are RIDLs, our bounds can be expressed as functions of the Laplacian eigenvalues of the underlying graph. Further rewriting of the bounds as functions of the graph's average effective resistance reveals that they are asymptotically tight. In Section 2.3 we illustrate the behavior of the upper and lower bounds for several sparse and dense graphs. Finally, Section 2.4 presents the conclusions and future work of this chapter.

Randomized consensus with additive noise

We consider a system composed by n agents holding scalar values x i ∈ R, i = 1, . . . , n, that interact through a randomized averaging algorithm in discrete time given by:

x(k + 1) = P (k)x(k), for all k ∈ Z ≥0 (2.1)
where

x(k) = [x 1 (k), . . . , x n (k)]
T is a vector with the states of the nodes of the network and P (k) is a stochastic matrix for all k. One of the most important results on dynamics (2.1) is the determination of sufficient conditions to achieve probabilistic consensus with independent and identically distributed (i.i.d.) matrices P (k).

Lemma 2.1 (Corollary 3.2 [START_REF] Fagnani | Randomized consensus algorithms over large scale networks[END_REF]) The dynamics (2.1) with i.i.d. matrices P (k) reaches consensus with probability 1 if:

P ii (k) > 0 with probability 1 for all i ∈ V, (2.2a)
G P has a globally reachable node,

where V = {1, . . . , n} and G P is the graph associated1 with the expected matrix

P = E[P (k)].
Condition (2.2a) is easy to satisfy in most cases since each agent usually has access to its own state. For undirected graphs, condition (2.2b) is equivalent to having a connected network. Moreover, for undirected graphs, the matrices P (k) of dynamics (2.1) are usually symmetric and, therefore, doubly stochastic: if matrices P (k) are doubly stochastic, then dynamics (2.1) reaches average consensus with probability 1.

Even if under ideal conditions it is possible to reach average consensus for undirected graphs, noise can affect the performance of the system. In order to study its effects, additive noise can be included in the dynamics [START_REF] Fagnani | Introduction to Averaging Dynamics over Networks[END_REF][START_REF] Garin | A survey on distributed estimation and control applications using linear consensus algorithms[END_REF], by defining:

x(k + 1) = P (k)x(k) + w(k), (2.3)
where w(k) is a vector of noise. It is natural to define a performance index to measure the deviation from the consensus point due to the perturbations as:

J noise := 1 n lim k→+∞ E x(k) - 1 n 11 T x(k) 2 , (2.4)
where 1 is a vector of ones. The purpose is studying this index, under the following standing assumption.

Assumption 2.1 The stochastic matrices P (k) are i.i.d., symmetric, and satisfy (2.2a)-(2.2b) for all k.

Our first contribution in this direction is showing that, under suitable conditions on the noise vector w(t), index J noise can be expressed in closed form.

Theorem 2.1 (Mean squared error index) Consider system (2.3) under Assumption 2.1 with a noise vector w(k) uncorrelated w.r.t. both i and k, with zero mean and variance σ 2 . Then, the noise index (2.4) can be expressed as:

J noise = σ 2 n vec T (I n )(I n 2 -K) -1 vec(I n ) -1 , (2.5)
where

K = E[P (k)Ω ⊗ P (k)], Ω = I n -1 n 11
T and I n is the identity matrix of size n.

Proof: For every k we have:

x(k) = Q(0, k -1)x(0) + k-1 s=0 Q(s + 1, k -1)w(s), where Q(ℓ, m) = P (m)P (m -1) • • • P (ℓ + 1)P (ℓ). Let us denote H = 1 n 11 T and δ(k) = x(k) -1 n 11 T x(k) = x(k) -Hx(k). Then, δ(k) = (Q(0, k -1) -H)x(0) + k-1 s=0 (Q(s + 1, k -1) -H)w(s).
The expectation of its squared norm can be expressed as:

E[∥δ(k)∥ 2 ] = E ∥(Q(0, k -1) -H)x(0)∥ 2 (2.6) + 2 k-1 s=0 x T (0)(Q(0, k -1) -H) T (Q(s + 1, k -1) -H)E[w(s)] + k-1 s,r=0 E [(Q(s + 1, k -1) -H)w(s)] T [(Q(r + 1, k -1) -H)w(r)] .
First, we analyze the third term of (2.6) which we denote by U . By applying the trace on the scalar U , we obtain

U = k-1 s,r=0 tr E (Q T (s + 1, k -1)Q(r + 1, k -1) -H) E[w(r)w T (s)] .
Due to the characteristics of the noise vector, we have:

E[w(r)w T (s)] = σ 2 I n , r = s 0, r ̸ = s ,
which yields:

U = σ 2 k-1 s=0 tr E Q T (s + 1, k -1)Q(s + 1, k -1) -H = σ 2 k-1 s=0 tr E Q T (s + 1, k -1)ΩQ(s + 1, k -1) . (2.7)
Then, we express the trace using the vectorization operator (tr(A T B) = vec T (A)vec(B)) as:

U = σ 2 k-1 s=0 vec T (I n )vec E Q T (s + 1, k -1)ΩQ(s + 1, k -1) .
Due to the idempotence of Ω and commutativity of the product P (s)Ω, we can write the matrix Q T (s + 1, k -1)ΩQ(s + 1, k -1) as:

Q T (s + 1, k -1)ΩQ(s + 1, k -1) = P (k -1) • • • P (s + 1)ΩP (s + 1) • • • P (k -1) = P (k -1) • • • P (s + 1)Ω 2 P (s + 1) • • • P (k -1) = P (k -1) • • • P (s + 1)ΩP (s + 1)Ω • • • P (k -1) . . . = P (k -1) • • • P (s + 1)P (s + 1)Ω • • • P (k -1)Ω. (2.8)
Then, we can apply the property vec(ABC) = (C T ⊗ A)vec(B) and we get:

U = σ 2 vec T (I n ) (I n ⊗ Ω)vec(I n ) + k-2 s=0 E (P (k -1)Ω ⊗ P (k -1))vec(P (k -2) . . . P 2 (s+1)Ω . . . P (k -2)Ω) .
Following an iterative procedure we obtain:

U = σ 2 vec T (I n ) (I n ⊗ Ω) + k-2 s=0
E [(P (k -1)Ω ⊗ P (k -1)) . . . (P (s+1)Ω ⊗ P (s+1))] vec(I n ).

Since the matrices P (s) are i.i.d., we get:

U = σ 2 vec T (I n ) k-1 s=1 E[P Ω ⊗ P ] s + (I n ⊗ Ω) vec(I n ) = σ 2 vec T (I n ) k-1 s=0 E[P Ω ⊗ P ] s + (I n ⊗ Ω) -I n 2 vec(I n ),
where we have avoided including the time dependence of matrices P (k). This abuse of notation will be made multiple times in the rest of this chapter. Notice that

(I n ⊗ Ω) -I n 2 = -(I n ⊗ 1 n 11 T ) and vec T (I n )(I n ⊗ 1 n 11 T )vec(I n ) = 1.
Let us recall K = E[P Ω ⊗ P ] and consider the limit k → +∞. By applying Jensen inequality to the spectral norm ∥ • ∥ 2 , we can see that

∥E[P Ω ⊗ P ]∥ 2 2 ≤ E[∥P Ω ⊗ P ∥ 2 2 ] = E[∥P Ω∥ 2 2 ∥P ∥ 2 2 ] < 1.
Therefore, ∞ s=0 K s = (I n 2 -K) -1 , and we have:

lim k→+∞ U = σ 2 vec T (I n )(I n 2 -K) -1 vec(I n ) -1 .
The second term of (2.6) is zero since the noise process has zero mean. The first term of (2.6), denoted by

F = E [∥(Q(0, k -1) -H)x(0)∥ 2 ]
, can be expressed as:

F = E [(Q(0, k -1) -H)x(0)] T [(Q(0, k -1) -H)x(0)] = E x T (0)(Q T (0, k -1)Q(0, k -1) -H)x(0) . Being Q T (0, k -1)Q(0, k -1)
the product of doubly stochastic matrices, per Theorem 4 in [START_REF] Tahbaz-Salehi | Consensus over ergodic stationary graph processes[END_REF] it converges to H with probability 1, so that lim n→∞ F = 0 by applying the Lebesgue Dominated Convergence Theorem. Finally, J noise is given by (2.5). ■ Remark 2.1 (Alternate expressions for J noise ) Since Ω is idempotent and commutes with any stochastic matrix, expression (2.8) could be arranged in two other ways, so that instead of (2.5) one would obtain the analogous expression where the symmetric matrix

K = E[P (k)Ω ⊗ P (k)] is replaced by E[P (k)Ω ⊗ P (k)Ω] or by E[P (k) ⊗ P (k)Ω].

Upper and lower bounds on the noise index

The calculation of (2.5) requires the determination of the n 2 ×n 2 matrix K which is impractical either for theoretical or numerical analysis. For this reason, we derive a lower bound J LB and an upper bound J UB for the noise index, which depend on the spectrum of symmetric n-dimensional matrices.

Theorem 2.2 (Upper and lower bounds) For the noise index (2.4) we have:

σ 2 n n i=2 1 1 -λ 2 i ( P ) ≤ J noise ≤ σ 2 n n i=2 1 1 -λ i ( P ) , (2.9)
where P = E[P (k)] and P = E[P 2 (k)].

Proof: In order to obtain the lower bound for the noise index, we can express the term (2.7) as:

U = σ 2 k-1 s=0 tr E Q T (s + 1, k -1)ΩQ(s + 1, k -1) = σ 2 k-1 s=0 E ∥Q(s + 1, k -1)Ω∥ 2 F ,
where ∥•∥ F is the Frobenius norm. Then we can apply Jensen inequality obtaining:

U ≥ σ 2 k-1 s=0 ∥E [Q(s + 1, k -1)Ω] ∥ 2 F .
Since the matrices P (s) are i.i.d., we can see that E[Q(s + 1, k -1)Ω] = P s Ω such that:

U ≥ σ 2 k-1 s=1 ∥ P s Ω∥ 2 F + ∥Ω∥ 2 F .
For a symmetric matrix P we have

∥P ∥ 2 F = ∥P 2 ∥ S 1 , where ∥A∥ S 1 = n i=1 |λ i (A)| is the Schatten 1-norm. Thus: U ≥ σ 2 k-1 s=1 ∥( P 2 Ω) s ∥ S 1 + ∥Ω∥ S 1 ≥ σ 2 k-1 s=1 ( P 2 Ω) s + Ω S 1 .
When k → +∞, we have:

lim k→+∞ U ≥ σ 2 ∞ s=0 ( P 2 Ω) s + Ω -I n S 1 = σ 2 (I n -P 2 Ω) -1 - 1 n 11 T S 1
.

By using the definition of the Schatten 1-norm we get:

lim k→+∞ U ≥ σ 2 n i=2 1 1 -λ 2 i ( P )
, which is the desired lower bound.

For the upper bound, we consider the summation element in (2.7), denoted by S:

S = tr E Q T (s + 1, k -1)(I n -H)Q(s + 1, k -1) = tr E Q T (s + 1, k -1)Q(s + 1, k -1) -tr (H) .
At this point, we can use the inequality

tr(E[Q T (s + 1, k -1)Q(s + 1, k -1)]) ≤ tr(E[P 2 ] k-s-1 ),
proved in [START_REF] Del Favero | A majorization inequality and its application to distributed Kalman filtering[END_REF]Cor. 5], which yields

U ≤ k-1 s=0 tr( P k-s-1 ) -tr (H) = k-1 s=0 tr( P s ) -tr (H)
Then, we obtain:

S ≤ tr( P s ) -tr (H) = tr P s -H = tr E P 2 (0)P 2 (1) • • • P 2 (s -1) -H = tr E P 2 (0)P 2 (1) • • • P 2 (s -1)Ω .
Since Ω is idempotent and P (s)Ω = ΩP (s), we can express the latter as:

S ≤ tr (E [P (0)ΩP (0) • • • P (s -1)ΩP (s -1)]) = tr (E [P (0)ΩP (0)] • • • E [P (s -1)ΩP (s -1)]) = tr (E [P (k)ΩP (k)] s ) .
Then, the sum in (2.7) is upper bounded by:

U ≤ σ 2 k-1 s=1
tr (E [P (0)ΩP (0)] s ) .

When k → +∞, we get:

lim k→+∞ U ≤ σ 2 tr ∞ s=0 E [P ΩP ] s -H = σ 2 tr (I n -E [P ΩP ]) -1 -H = σ 2 n i=2 1 1 -λ i ( P ) ,
which gives the desired upper bound. ■ Remark 2.2 (Deterministic case) For the deterministic dynamics x(k + 1) = P x(k)+w(t), that is, P (k) = P , the noise index is given by

J noise = σ 2 n n i=2 1 1-λ 2 i ( P ) .
In this case, the upper and lower bounds of Theorem 2.2 coincide and recover the well-known results about deterministic consensus [START_REF] Fagnani | Introduction to Averaging Dynamics over Networks[END_REF][START_REF] Garin | A survey on distributed estimation and control applications using linear consensus algorithms[END_REF]. These questions have also been extensively studied in continuous-time by the closely-related notion of network coherence [START_REF] Patterson | Leader selection for optimal network coherence[END_REF][START_REF] Bamieh | Coherence in largescale networks: Dimension-dependent limitations of local feedback[END_REF].

Noise index in OMAS

When studying OMAS, one is confronted with the issue of agents continuously joining and leaving the system. When the pool of possible participating agents is known, arrivals and departures can be equivalently seen as activations and deactivations [START_REF] Varma | Open multi-agent systems with discrete states and stochastic interactions[END_REF]. In such a case, one can further assume that the interaction network is a subgraph, induced by the active nodes, of a larger network of potential interactions and specific sampling procedures can be used to generate the stochastic matrices.

Let a set V of n nodes be given, together with an undirected graph Ḡ that has V as node set and whose adjacency matrix is denoted by Ā. At each time k, every node has an independent probability p of being active; at each time k, the current graph G(k) is the subgraph of Ḡ that is induced by the set of active nodes. Equivalently, the current adjacency matrix A(k) is generated as

A(k) = Γ(k) Ā Γ(k), (2.10)
where

Γ(k) = diag[γ 1 (k), . . . , γ n (k)
] is a diagonal matrix with n independent identically distributed Bernoulli random variables corresponding to each node. The value of the random variable γ i (k) indicates the state of agent i at time k such that when γ i (k) = 1, the node is active. Given the adjacency matrix A(k), we then define the time-varying doubly stochastic matrix given by:

P (k) = I n -ϵ(D(k) -A(k)), (2.11)
where D(k) is the diagonal matrix of the degrees of the induced graph G(k), 0 < ϵ < 1 dmax( Ḡ) , and d max ( Ḡ) is the maximum degree of the graph Ḡ. After defining the time-varying Laplacian L(k) = D(k) -A(k), it becomes natural to refer to stochastic matrices sampled according to (2.11)-(2.10) as Randomly Induced Discretized Laplacians (RIDL). It is clear that RIDL can be used to study OMAS under the assumption that the activation of each agent is an independent event, as illustrated in Fig. 2.1, where at each iteration we have a different number of active agents n a .

We now study J noise for RIDL matrices. For an exact calculation as per (2.5), it would be necessary to compute the matrix K = E[P (k)Ω ⊗ P (k)]. To this end, we can observe that The second term only depends on P , which is given by

K = E[P (k)Ω ⊗ P (k)] = E P (k) I n - 1 n 11 T ⊗ P (k) = E[P (k) ⊗ P (k)] -E 1 n 11 T ⊗ P (k) = E[P (k) ⊗ P (k)] - 1 n 11 T ⊗ P .
P = I n -ϵp 2 ( D -Ā) = I n -ϵp 2 L, (2.12)
where D and L are the degree and Laplacian matrices of graph Ḡ, but the first term

E[P (k) ⊗ P (k)] = I n 2 -ϵp 2 (I n ⊗ L + L ⊗ I n ) + ϵ 2 E[L(k) ⊗ L(k)]
cannot be written as a function of P or L. Hence, J noise cannot be written as a function of P . Nevertheless, the upper and lower bounds of Theorem 2 can be expressed as functions of the eigenvalues of P .

Proposition 2.1 (Bounds for Discretized Laplacian) For the time-varying stochastic matrix (2.11) generated with (2.10), the lower bound J LB and upper bound J UB in (2.9) become:

J LB = σ 2 ϵp 2 n n i=2 1 2λ i ( L) -ϵp 2 λ 2 i ( L) , (2.13)
J UB = σ 2 ϵp 2 n n i=2 1 2(1 + ϵp -ϵ)λ i ( L) -ϵpλ 2 i ( L) , (2.14
)

where 0 = λ 1 ( L) < λ 2 ( L) ≤ . . . ≤ λ n ( L) are the eigenvalues of L.
Proof: For the lower bound, we consider the expression (2.12). Since L is symmetric, the eigenvalues of P are given by λ i ( P ) = 1 -ϵp 2 λ i ( L) and we obtain (2.13).

For the upper bound, P can be written as:

P = E[(I n -ϵL) 2 ] = I n -2ϵp 2 L + ϵ 2 E[L 2 ] = I n -2ϵp 2 L + ϵ 2 (E[D 2 ] -E[DA] -E[AD] + E[A 2 ]) Since E[γ i (k)] = E[γ 2 i (k)] = p, we have: E[D 2 ] = (p 2 -p 3 ) D + p 3 D2 , E[DA] = p 3 D Ā + (p 2 -p 3 ) Ā, E[AD] = p 3 Ā D + (p 2 -p 3 ) Ā, E[A 2 ] = (p 2 -p 3 ) D + p 3 Ā2
, which yields:

P = I n + 2ϵp 2 (ϵ -ϵp -1) L + ϵ 2 p 3 L2 .
The eigenvalues of P are λ i ( P ) = 1 + 2ϵp 2 (ϵ -ϵp -1)λ i ( L) + ϵ 2 p 3 λ 2 i ( L) and we obtain (2.14).

■ From expressions (2.13) and (2.14) we can observe that for the computation of the bounds for J noise , we only need to know the topology of G P through L.

Furthermore, inspired by the work in [START_REF] Lovisari | Resistance-based performance analysis of the consensus algorithm over geometric graphs[END_REF], we relate J noise with the well-studied Average Effective Resistance (or Kirchhoff Index) of the graph. The latter can indeed be expressed [START_REF] Fagnani | Introduction to Averaging Dynamics over Networks[END_REF]Chapter 5] as a function of Laplacian eigenvalues by

R ave := 1 n n i=2 1 λ i ( L)
.

Proposition 2.2 (Bounds and average effective resistance) For the noise index (2.4) with the time-varying stochastic matrix (2.11) generated by (2.10) we have:

σ 2 2p 2 R ave ϵ ≤ J noise ≤ σ 2 2p 3 (1 -r ϵ ) R ave ϵ ,
where r ϵ = ϵd max ( Ḡ) < 1.

Proof: For the lower bound, we immediately observe from (2.13) that

J LB ≥ σ 2 ϵp 2 n n i=2 1 2λ i ( L) = σ 2 2ϵp 2 R ave .
For the upper bound, we remind λ i ( L) ≤ 2d max ( Ḡ) and obtain:

2(1 + ϵp -ϵ)λ i ( L) -ϵpλ 2 i ( L) = 2ϵ(p -1)λ i ( L) + (2 -ϵpλ i ( L))λ i ( L) ≥ 2ϵ(p -1)λ i ( L) + 2 -ϵp 2d max ( Ḡ) λ i ( L) ≥ 2p(1 -r ϵ )λ i ( L),
which yields:

J UB ≤ σ 2 ϵp 2 n n i=2 1 2p(1 -r ϵ )λ i ( L) = σ 2 2ϵp 3 (1 -r ϵ ) R ave ,
thereby completing the proof. ■

The bounds of Proposition 2.2 can also be equivalently expressed as

σ 2 2p 2 r ϵ d max ( Ḡ)R ave ≤ J noise ≤ σ 2 2p 3 r ϵ (1 -r ϵ ) d max ( Ḡ)R ave .
Since only d max ( Ḡ) and R ave depend on the graph (and therefore on n), this form of the bounds implies that

J noise = Θ(d max ( Ḡ)R ave ) as n → ∞,
where Θ(•) corresponds to the Θ-notation 2 , and proves that the bounds are asymptotically tight for RIDLs. For instance, for sparse graphs with bounded degree, the rate of growth of J noise is equal to the rate of R ave . Instead, for dense graphs where d max ( Ḡ) is linear in n, the rate of growth of J noise is the rate of nR ave .

Particular graph topologies

In this section, we illustrate the behavior of the noise index for RIDL matrices associated to graphs with particular structures, which allow for the explicit computation of the bounds in Theorem 2.2 and Proposition 2.1. We consider ϵ = r ϵ /d max ( Ḡ) with r ϵ ∈ (0, 1) and we perform computations for networks with 3 ≤ n ≤ 100, p = 0.9, and σ 2 = 1.

Sparse graphs

Star graphs: for the star graph S n we have d max (S n ) = n -1 and we choose ϵ = r ϵ /(n -1). The eigenvalues of the Laplacian matrix are given by λ i ( L) = 1 for i = 2, . . . , n -1 and λ n ( L) = n, so that the bounds are:

J LB = σ 2 ϵp 2 n n -2 2 -ϵp 2 + 1 n(2 -ϵp 2 n) ; (2.15)
J UB = σ 2 ϵp 2 n n -2 ϵp + 2 -2ϵ + 1 n(2ϵp + 2 -2ϵ -ϵpn) , (2.16)
implying that

J noise = σ 2 2r ϵ p 2 n + o(n) as n → ∞,
2 For a given function g(n), we denote by Θ(g(n)) the set of functions [START_REF] Cormen | Introduction to Algorithms[END_REF] Θ(g(n)) = {f (n) : there exist positive constants c 1 , c 2 and n 0 such that Line and grids: for the line graph P n , d max (P n ) = 2 and we select ϵ = 0.8/2 for all n. The eigenvalues of the Laplacian matrix are given by λ i ( L) = 2 -

0 ≤ c 1 g(n) ≤ f (n) ≤ c 2 g(n) for all n ≥ n 0 }. (a) J noise (b) Relative error
2 cos π(i-1)
n for i = 2, . . . , n such that the corresponding bounds are:

J LB = σ 2 4ϵp 2 n n-1 i=1 1 1 -ϵp 2 -ϵp 2 cos 2 πi n + (2ϵp 2 -1) cos πi n ;
(2.17)

J UB = σ 2 4ϵp 2 n n-1 i=1 1 1 -ϵ -ϵp cos 2 πi n + (ϵp + ϵ -1) cos πi n , (2.18)
and the rate of growth of the noise index is

J noise = Θ(n).

Dense graphs

Complete graph: since for a complete graph

K n , d max (K n ) = n -1, we select ϵ = r ϵ /(n -1)
. The eigenvalues of the Laplacian matrix of the complete graph K n are λ i ( L) = n for i = 2, . . . , n such that the bounds for J noise are:

J LB = σ 2 (n -1) ϵp 2 n 2 (2 -ϵp 2 n) ; (2.19)
J UB = σ 2 (n -1) ϵp 2 n 2 (2 + 2ϵp -2ϵ -ϵpn) , (2.20)
implying a rate of growth J noise = Θ(1). The limit when n → ∞ is given by:

lim n→∞ J LB = σ 2 p 2 r ϵ (2 -p 2 r ϵ ) ; lim n→∞ J UB = σ 2 p 2 r ϵ (2 -pr ϵ )
.

For the computations we choose r ϵ = 0.8 and we obtain the results in Figures 2.5 Erdős-Rényi graph: each edge between two nodes is generated with a probability p er , independently of any other edge. As an example, figures 2.6(a) and 2.6(b) present the results of the computations for the Erdős-Rényi graph p er = 0.8 and ϵ = 0.8/(n -1).

Conclusion

In this chapter, we analyzed the influence of noise in random consensus with an application to OMAS by considering the random sampling of active nodes. We derived an expression for a noise index that measures the deviation from the consensus point and an upper and lower bound that depend on the spectrum of expected matrices related to the network. For particular cases of OMAS characterized by activations and deactivations of agents, we defined Randomly Induced Discretized Laplacians and expressed the bounds as functions of the eigenvalues of the expected Laplacian matrix and, more specifically, of its average effective resistance. For sparse graphs, such as stars and paths, the noise index grows linearly in n and the lower bound presents a better approximation of the real value. For dense graphs, such as complete, Erdős-Rényi and graphs sampled from graphons, the noise index is bounded in n and the upper bound is a better approximation.

For future work, we identify three main possible directions of the research. Regarding the rich family of graphs sampled from graphons, the eigenvalues of the Laplacian matrix are related to the degree function of the graphon [START_REF] Vizuete | The Laplacian spectrum of large graphs sampled from graphons[END_REF]. It would be interesting to consider approximations of the bounds by using appropriate definitions that include the degree function of the graphon as in the case of the average effective resistance in [START_REF] Vizuete | The Laplacian spectrum of large graphs sampled from graphons[END_REF].

A key assumption for the derivation of the results in this chapter is the independence of active sets across time. However, it would be important to consider the scenario where only few agents change their states from active to inactive or vice versa while the activation state of other agents remain invariant. In this case, independence in the stochastic matrices will not hold and a certain degree of statistical dependence in the generation of the stochastic matrices should be considered.

The metric (2.4) used for the analysis of OMAS was normalized considering the total number of potential agents in the system, i.e., n. We would like to analyze an alternative definition of this index function that could be more informative for OMAS, by considering only the current number of active agents at each iteration, i.e., n a , as the normalization factor.

Chapter 3 Weighted gradient descent with packet losses and agent replacements

The resource allocation (RA) is an important problem in optimization where a budget must be optimally distributed among multiple entities or activities [START_REF] Ibaraki | Resource Allocation Problems: Algorithmic Approaches[END_REF][START_REF] Patriksson | A survey on the continuous nonlinear resource allocation problem[END_REF]. In the case of multi-agent systems, the objective is to find an optimal distribution of the budget among the agents while each agent aims to minimize the whole cost function but it can only access its own local cost and state and those of its neighbors [START_REF] Nedić | Improved convergence rates for distributed resource allocation[END_REF][START_REF] Doan | Distributed resource allocation over dynamic networks with uncertainty[END_REF]. Applications include energy resources [START_REF] Dominguez-Garcia | Decentralized optimal dispatch of distributed energy resources[END_REF], smart grids [START_REF] Dai | Distributed Q-learning algorithm for dynamic resource allocation with unknown objective functions and application to microgrid[END_REF], actuator networks [START_REF] Teixeira | Distributed actuator reconfiguration in networked control systems[END_REF] games [START_REF] Liang | Distributed Nash equilibrium seeking for aggregative games with coupled constraints[END_REF], and distributed computer systems [START_REF] Kurose | A microeconomic approach to optimal resource allocation in distributed computer systems[END_REF]. The majority of the algorithms are designed for ideal scenarios with perfect communication and a symmetric exchange of information through the network that maintain the budget fixed (keeping the budget fixed is indeed the constraint of the optimization problem). However, the RA problem can also be solved in digraphs by introducing additional variables to preserve the constraint, which increases the complexity of the algorithms [START_REF] Xu | A distributed algorithm for resource allocation over dynamic digraphs[END_REF][START_REF] Deng | Distributed continuous-time algorithms for resource allocation problems over weight-balanced digraphs[END_REF]. One of the most important algorithms used to solve the RA problem in a multi-agent system is the weighted gradient descent proposed in [START_REF] Xiao | Optimal scaling of a gradient method for distributed resource allocation[END_REF], where the update rule makes use of a matrix associated to the network that preserves the constraint.

The classical formulation of the RA problem assumes that the set of agents remains invariant during the optimization process; however, this assumption could not be valid for some scenarios. Consider, for instance, the case of distributed energy resources where a fixed amount of energy must be supplied by a network of devices in an optimal way [START_REF] Dominguez-Garcia | Decentralized optimal dispatch of distributed energy resources[END_REF]. Nevertheless, in real life, some of the devices might experience failures with higher probability as the system size increases, or change their operating point due to environmental conditions, generating an open system where the cost functions of agents can change. Optimization in OMAS is a recent topic research, where two different types of formulations have been considered. In [START_REF] Hsieh | Optimization in open networks via dual averaging[END_REF][START_REF] Hsieh | Multi-agent online optimization with delays: Asynchronicity, adaptivity, and optimism[END_REF], the authors study the case of different active agents interacting in the network for the implementation of an optimization algorithm, while the size of the system is fixed and all the possible cost functions are known. In [START_REF] Hendrickx | Stability of decentralized gradient descent in open multi-agent systems[END_REF], the openness of the system is characterized by replacements, where new cost functions can be assigned to the agents, implying a possible change of the global minimizer.

In addition to potential changes in the set of agents, communication networks can suffer from errors in the communications that cause the loss of information. In this case, even if the graph associated with the communication network was originally undirected, interactions between the agents can effectively be asymmetric when the packet loss occurs only in one direction. The symmetry of the system being broken, relevant quantities may fail to be preserved in many scenarios [START_REF] Frasca | On the mean square error of randomized averaging algorithms[END_REF][START_REF] Patterson | Distributed average consensus with stochastic communication failures[END_REF][START_REF] Fagnani | Average consensus with packet drop communication[END_REF]. In the case of the average estimation under packet losses, [START_REF] Frasca | Large network consensus is robust to packet losses and interferences[END_REF] derived upper bounds for the mean squared error to evaluate the deviation from the initial average for general graph topologies. While a possible deviation of the initial average value is not necessarily problematic in the average estimation where the main objective is to reach consensus in most cases, for optimization problems, similar errors can be more consequential, since they may imply the violation of the constraint. For example in the case of the RA problem, agents could perform unnecessary tasks or the sum of the individual tasks may not satisfy the global demand. Indeed, the analysis of the constraint violation in complex scenarios has been the object of many works, specially in the framework of online optimization, where it is assumed that cost functions change in time [START_REF] Yi | Regret and cumulative constraint violation analysis for online convex optimization with long term constraints[END_REF].

In this chapter, we analyze the effects of packet losses in the performance of the weighted gradient descent algorithm to solve the RA problem. In Section 3.1 we introduce the formal formulation of the RA problem, including the assumptions used for the derivation of the results of this chapter. Section 3.2 presents the problem of packet losses in the weighted gradient descent algorithm and the conditions to guarantee convergence of the algorithm to the minimizer at least in expectation. We define appropriate metrics to measure the deviation from the constraint and the error of the expected value of the cost function and we derive upper bounds. In Section 3.3 we analyze the weighted gradient descent in OMAS where replacements of the agents can take place and we show that the constraint error can diverge for this particular scenario. Finally, conclusions and future work are exposed in Section 3.4.

Resource Allocation Problem

In multi-agent systems, the RA problem is formulated as the minimization of an objective function f that is separable in local costs f i : R d → R held by the agents, subject to an equality constraint on the weighted sum of the states x i ∈ R d with respect to a budget composed by the sum of the demands of the agents

d i ∈ R d .
The problem can then be written as

min x=[x T 1 ,...,x T n ] T ∈R nd f (x) = n i=1 f i (x i ) subject to n i=1 a i x i = n i=1 d i , (3.1)
where a i > 0 is the weight of agent i to satisfying the constraint. When the agents are isolated, their objective is to minimize the corresponding local cost functions f i .

In the RA problem (3.1), each agent aims at minimizing the global cost function f while guaranteeing the budget. However, since each agent only has access to local information, they need to exchange the current states of the gradients to achieve the goal.

If the budget is fixed in time n i=1 d i = b, then the constraint in (3.1) can be equivalently expressed as a ⊤ ⊗ I d x = b, where a = [a 1 , . . . , a n ] T and we remind that ⊗ denotes the Kronecker product. The feasible set of (3.1) is thus given by

S a,b := x ∈ R nd | a ⊤ ⊗ I d x = b . (3.2)
For the particular case d = 1, the resource allocation constraint can be expressed as ⟨a, x⟩ = b. We make the following classical assumption on the local cost functions. Since all the local cost functions f i are α-strongly convex and β-smooth, the global cost function f from (3.1) also satisfies f ∈ F α,β . As a consequence, Assumption 3.1 guarantees that the solution of the problem (3.1) denoted by x * is unique. Therefore, for some ζ * ∈ R d , a necessary and sufficient condition for the optimality of x * is ∇f

Assumption 3.1 Each function f i is continuously differentiable, α-strongly con- vex (i.e., f i (x) -α 2 ||x|| 2 is convex) and β-smooth (i.e., ||∇f i (x) -∇f i (y)|| ≤ β ||x -y||, ∀x, y ∈ R d ).
(x * ) = (a ⊤ ⊗ I d ) ⊤ ζ * = a ⊗ ζ * . (3.3)
In this chapter, we restrict to the specific case defined by the following assumptions, while in Chapter 4 we consider the more general case of the RA problem. For this particular setting, the RA problem (3.1) can be formulated as

min x=[x 1 ,...,xn]∈R n f (x) = n i=1 f i (x i ) subject to n i=1 x i = b, (3.4)
and the optimality condition implies that

1 T x * = b and ∇f * = ζ * 1 for some ζ * ∈ R where ∇f * = ∇f (x * ).
The feasible set of (3.4) is given by

S 1,b := x ∈ R n |1 T x = b . (3.5)

Packet losses

We consider that agents interact through a connected undirected graph G = (V, E), where V = {1, . . . , n} is the set of agents and E ⊆ V × V is the set of edges. To the graph G we associate an adjacency matrix A = [a ij ] where a ij = 1 if (i, j) ∈ E and a ij = 0 otherwise. We assume that the graph does not have self-loops (i.e., a ii = 0). The Laplacian matrix of the graph G is defined as

L = D -A where D = diag[d 1 , . . . , d n ]
and d i is the degree of agent i (i.e., d i is the ith row-sum of A). Notice that L is a symmetric matrix that satisfies L = L T .

To solve problem (3.4) at each iteration we apply the weighted gradient descent algorithm with positive step-size h [START_REF] Xiao | Optimal scaling of a gradient method for distributed resource allocation[END_REF]:

x(k + 1) = x(k) -hL∇f (x(k)), (3.6)
where L determines the exchange of information between the agents according to the network topology.

For the RA problem formulated in 3.4 with local cost functions satisfying Assumption 3.1 we provide an upper bound for the step-size h to guarantee linear convergence of the cost function.

f (x(k + 1)) -f (x * ) ≤ (1 -αhλ 2 (L)(2 -hβλ n (L))) (f (x(k)) -f (x * )) , (3.7)
where λ n (L) and λ 2 (L) are the largest and the second smallest eigenvalues of L respectively.

Proof: Since the cost function is β-smooth we have:

f (x(k + 1)) ≤ f (x(k)) + ⟨∇f (x(k)), x(k + 1) -x(k)⟩ + β 2 ||x(k + 1) -x(k)|| 2 = f (x(k)) + ⟨∇f (x(k)), -hL∇f (x(k))⟩ + β 2 ||hL∇f (x(k))|| 2 ≤ f (x(k)) -h L 1/2 ∇f (x(k)) 2 + βh 2 λ n (L) 2 L 1/2 ∇f (x(k)) 2 = f (x(k)) - h 2 (2 -βhλ n (L)) L 1/2 ∇f (x(k)) 2 .
Using a vector projection on 1, the gradient can be expressed as:

∇f (x(k)) = ξ ⊥ 1 + ∇f (x(k)) ⊥ , where ξ ⊥ = 1 n ⟨∇f (x(k)), 1⟩. Since h ≤ 2 βλn(L)
, we can upper bound the last expression by:

f (x(k + 1)) ≤ f (x(k)) - h 2 (2 -βhλ n (L)) ξ ⊥ L 1/2 1 + L 1/2 ∇f (x(k)) ⊥ 2 ≤ f (x(k)) - h 2 (2 -βhλ n (L)) L 1/2 ∇f (x(k)) ⊥ 2 ≤ f (x(k)) - hλ 2 (L) 2 (2 -βhλ n (L)) ||∇f (x(k)) ⊥ || 2 , (3.8)
where we used the fact that L 1/2 1 = 0 since 1 belongs to the kernel of L 1/2 . The expression (3.8) shows the convergence of the algorithm when the step-size satisfies h ≤ 2 βλn(L) . Furthermore, by using the properties of α-strongly convex functions, we have the following inequality [80, Proof of Theorem 4.1]:

f (x(k)) -f (x * ) ≤ 1 2α ||∇f (x(k)) ⊥ || 2 . (3.9)
Finally, by combining (3.8) and (3.9) we obtain:

f (x(k + 1))-f (x * ) ≤ f (x(k))-f (x * ) - hλ 2 (L) 2 (2 -βhλ n (L)) ||∇f (x(k)) ⊥ || 2 ≤ f (x(k))-f (x * ) -αhλ 2 (L)(2 -βhλ n (L)) (f (x(k))-f (x * )) .
■ In this analysis, we consider that the communication network is unreliable and some messages can be lost. Following the approach in [START_REF] Frasca | Large network consensus is robust to packet losses and interferences[END_REF] we assume that at each time instant k ∈ Z ≥0 , we have a matrix

Γ(k) = [Γ(k) ij ] ∈ {0, 1} n×n composed of Bernoulli random variables γ ij
k independent across k, i, j with P(γ ij k = 1) = p and p ∈ (0, 1). Then, for an edge a ij = 1, if γ ij k = 1 then the link a ij is active, otherwise there was a failure of the link a ij and the data transmitted was lost. If a ij = 0, the value of γ ij k does not have any effect.

The topology of the network with possible packet losses at each time instant is thus defined by:

A(k) = A ⊙ Γ(k); (3.10) L(k) = D(k) -A(k), (3.11)
and the weighted gradient descent algorithm becomes

x(k + 1) = x(k) -hL(k)∇f (x(k)). (3.12)
Due to the asymmetric loss of information, the graph topology at each time instant is effectively directed with a non-symmetric Laplacian matrix L(k). Because of this lack of symmetry, the constraint can be violated during an update of the states of the network using (3.12) and hence, the minimizer of the problem (3.1) will not be reached. Nevertheless, it is desirable that algorithm (3.12) solves (3.1) at least in expectation. Since Γ(k) is independent in time and E [L(k)] = pL, the expected dynamics is given by:

E [x(k + 1)] = E [x(k)] -hpLE [∇f (x(k))] . (3.13)
Dynamics (3.13) preserves the constraint in expectation since:

1 T E [x(k + 1)] = 1 T (E [x(k)] -hpLE [∇f (x(k))]) = 1 T E [x(k)] . (3.14)
However, the minimizer x * is not necessarily a stationary point of (3.13), as shown by the following discussion. Let us assume that E [x(k)] = x * , so that we have

E [x(k + 1)] = x * -hpLE [∇f (x(k))] .
Under the assumption that the local functions are β-smooth and α-strongly convex we cannot guarantee that

E [∇f (x(k))] = ∇ (f (E [x(k)])) = ζ1.
For this reason, we restrict the set of possible cost functions to continuous differentiable piecewise quadratic functions, which are β-smooth and α-strongly convex. In this case, the gradient is a vector of piecewise linear functions [∇f (x)] i = m i x i + c i and we have

E [∇f (x(k))] = ∇f (E [x(k)]) . (3.15)
Then, the expected dynamics (3.13) becomes:

E [x(k + 1)] = E [x(k)] -hpL∇f (E [x(k)]) , (3.16)
which corresponds to the original weighted gradient descent, thereby guaranteeing that the constraint is preserved at each time instant and that the minimizer is reached such that

lim k→∞ E [x(k)] = x * and lim k→∞ f (E [x(k)]) = f * ,
where f * = f (x * ). However, a similar equality is not true for the expected value of the cost function. Since the functions are convex, by using Jensen's inequality we obtain:

E [f (x(k))] ≥ f (E [x(k)]
) , and by taking the limit (whose existence will be discussed in the next section) we get:

lim k→∞ E [f (x(k))] ≥ f * , (3.17)
where the inequality is in general strict.

Constraint violation

Although the constraint is preserved in expectation according to (3.14), this fact does not imply that the constraint is also preserved for any particular realization of the process. We thus aim to measure the deviation from the constraint by deriving an upper bound for the constraint violation metric:

J constr := lim k→∞ E 1 T x(k) -b 2 , (3.18)
which corresponds to the asymptotic mean squared error. Before presenting the main result of this subsection, we recall a lemma derived in [START_REF] Frasca | Large network consensus is robust to packet losses and interferences[END_REF].

E L(k) T L(k) = p 2 L 2 + 2p(1 -p)L (3.19) E L(k) T 11 T L(k) = 2p(1 -p)L. (3.20)
Theorem 3.1 Consider the weighted gradient descent given by (3.10), (3.11) and (3.12). Under Assumption 3.5, for any positive scalar h ≤ 2 β(pλn(L)+2-2p) , the constraint violation metric satisfies:

J constr ≤ 4h(1 -p) 2 -hβ(pλ n (L) -2p + 2) (f (x 0 ) -f * ) . (3.21)
Proof: Let us denote by H(k) = 1 T x(k) -b . For (3.12) we denote by {F k } k∈Z ≥0 the filtration of σ-algebras generated by the process x(k). We compute the expectation conditioned upon the filtration generated by x(k):

E H 2 (k + 1)|F k = E 1 T (x(k) -hL(k)∇f (x(k))) -b 2 = H 2 (k) + h 2 ∇f (x(k)) T E L(k) T 11 T L(k) ∇f (x(k)) -hH(k)1 T E [L(k)] ∇f (x(k)) = H 2 (k) + h 2 ∇f (x(k)) T E L(k) T 11 T L(k) ∇f (x(k)),
where we have used the fact that 1 T L = 0. Then, we use (3.20) to obtain:

E H 2 (k + 1)|F k = H 2 (k) + 2h 2 p(1 -p)∇f (x(k)) T L∇f (x(k)) = H 2 (k) + 2h 2 p(1 -p) L 1/2 ∇f (x(k)) 2 .
If we take the total expectation we get:

E E H 2 (k + 1)|F k =E H 2 (k) + 2h 2 p(1 -p)E L 1/2 ∇f (x(k)) 2 .
By using a recursive argument we have that at the time instant k:

E H 2 (k) = H 2 0 + 2h 2 p(1 -p) k-1 i=0 E L 1/2 ∇f (x i ) 2 = 2h 2 p(1 -p) k-1 i=0 E L 1/2 ∇f (x i ) 2 , (3.22)
where we used the fact that at the initial time k = 0, the constraint is not violated H 0 = 0. Now, let us consider the inequality corresponding to β-smooth functions:

f (x(k + 1)) ≤ f (x(k)) + β 2 ||x(k) -x(k + 1)|| 2 + ⟨∇f (x(k)), x(k + 1) -x(k)⟩ = f (x(k)) + h 2 β 2 ||L(k)∇f (x(k))|| 2 -h⟨∇f (x(k)), L(k)∇f (x(k))⟩ = f (x(k)) + h 2 β 2 ⟨∇f (x(k)), L(k) T L(k)∇f (x(k))⟩ -h⟨∇f (x(k)), L(k)∇f (x(k))⟩.
We compute the expectation given the filtration generated by x(k) and we use (3.19) to get:

E [f (x(k + 1))|F k ] ≤ f (x(k)) + h 2 β 2 ⟨∇f (x(k)), E L(k) T L(k) ∇f (x(k))⟩ -h⟨∇f (x(k)), E [L(k)] ∇f (x(k))⟩ = f (x(k)) + h 2 β 2 ⟨∇f (x(k)), (p 2 L 2 + 2p(1 -p)L)∇f (x(k))⟩ -h⟨∇f (x(k)), pL∇f (x(k))⟩ We use the upper bound ||L∇f (x(k))|| 2 ≤ λ n (L) L 1/2 ∇f (x(k))
2 and we obtain:

E [f (x(k + 1))|F k ] ≤ f (x(k))-hp 1 -hβ(1 -p) - hβpλ n (L) 2 L 1/2 ∇f (x(k)) 2 .
We denote ξ = hp 1 -hβ(1 -p) -hβpλn(L)

2
and we compute the total expectation to get:

E [E [f (x(k + 1))|F k ]] ≤ E [f (x(k))] -ξE L 1/2 ∇f (x(k)) 2 . (3.23) Since h ≤ 2 β(pλn(L)+2-2p) we guarantee that E [f (x(k + 1))] ≤ E [f (x(k))] such that: E L 1/2 ∇f (x(k)) 2 ≤ 1 ξ (E [f (x(k))] -E [f (x(k + 1))]) .
If we take the sum over the time instants s we have:

k-1 s=0 E L 1/2 ∇f (x(s)) 2 ≤ 1 ξ (f (x 0 ) -E [f (x(k))]) . Notice that since E [f (x(k))] is monotonically decreasing and E [f (x(k))] ≥ f (x * ),
where x * is the minimizer of f (x) without constraints, the limit lim k→∞ E [f (x(k))] exists and we get:

lim k→∞ k-1 s=0 E L 1/2 ∇f (x(s)) 2 ≤ 1 ξ f (x 0 ) -lim k→∞ E [f (x(k))] ≤ 1 ξ (f (x 0 ) -f * ) , (3.24)
where we used the fact that By inspecting the right-hand side of (3.21), we can make several observations. The upper bound is decreasing in p and

E [f (x(k))] ≥ f (E [x(k)]).
lim p→1 J constr = 0, (3.25)
which is consistent with the fact that the constraint is always preserved with no communication losses. In the case of small values of p, which corresponds to a network exposed to many failures, if we consider h = 2ϵ β(pλn(L)+2-2p) with ϵ ∈ (0, 1) we obtain:

lim p→0 J constr ≤ 2ϵ β(1 -ϵ) (f (x 0 ) -f * ) . ( 3

Error on the expected cost function

Due to the packet losses, f (x(k)) is a random process whose evolution depends on the different realizations of L(k). When the algorithm (3.12) is executed, it is important to know the distribution of f (x(k)), to evaluate the impact of the perturbations on the performance of (3.12) and identify the additional cost corresponding to the deviation from the ideal cost f * . From (3.17) we can observe that this quantity is greater than f * . Since the objective is to remain as close as possible to f * , we define the cost function metric:

J funct := lim k→∞ E [f (x(k))] -f (x * ), (3.27)
which aims to measure this gap due to the packet losses during the application of the weighted gradient descent algorithm. It is important to notice that (3.27) corresponds to the limit of the so called Jensen gap [START_REF] Abramovich | Some new estimates of the 'Jensen gap[END_REF].

Theorem 3.2 Consider the weighted gradient descent given by (3.10), (3.11) and (3.12). Under Assumption 3.5, for any positive scalar h ≤ 2 β(pλn(L)+2-2p) , the cost function metric satisfies:

J funct ≤ 1 - 2 -hβ (pλ n (L) -2p + 2) 2 -hαpλ 2 (L) (f (x 0 ) -f * ) .
Proof: From (3.23), we have that the total expectation of f (x(k + 1)) satisfies:

E [f (x(k + 1))] ≤ E [f (x(k))] -ξE L 1/2 ∇f (x(k)) 2 ,
and by using Jensen's inequality and (3.15) we get:

E [f (x(k + 1))] ≤ E [f (x(k))] -ξ E L 1/2 ∇f (x(k)) 2 = E [f (x(k))] -ξ L 1/2 ∇f (E [x(k)]) 2 .
If we take the sum over the time instants s we obtain:

E [f (x(k))] ≤ f (x 0 ) -ξ k-1 s=0 L 1/2 ∇f (E [x(s)]) 2 ,
and the limit satisfies:

lim k→∞ E [f (x(k))] ≤ f (x 0 ) -ξ lim k→∞ k-1 s=0 L 1/2 ∇f (E [x(s)]) 2 . (3.28)
Now, we consider the dynamics in expectation given by (3.13), which satisfies the inequality of α-strongly convex functions:

f (E [x(k + 1)]) ≥ f (E [x(k)]) + ⟨∇f (E [x(k)]), E [x(k + 1)] -E [x(k)]⟩ + α 2 ||E [x(k + 1)] -E [x(k)]|| 2 = f (E [x(k)]) + ⟨∇f (E [x(k)]), -hpL∇f (E [x(k)])⟩ + α 2 ||hpL∇f (E [x(k)])|| 2 = f (E [x(k)])-hp L 1/2 ∇f (E [x(k)]) 2 + h 2 αp 2 2 ||L∇f (E [x(k)])|| 2
Notice that due to the Min-Max Theorem, for any vector z ∈ R n , it holds

||Lz|| 2 = L 1/2 L 1/2 z 2 ≥ λ 2 (L) L 1/2 z 2 , (3.29)
since L 1/2 z is orthogonal to 1. Hence we obtain: ) we get:

f (E [x(k + 1)]) ≥ f (E [x(k)]) -hp L 1/2 ∇f (E [x(k)]) 2 + h 2 αp 2 λ 2 (L) 2 L 1/2 ∇f (E [x(k)]) 2 = f (E [x(k)]) -hp 1 - hαpλ 2 (L) 2 L 1/2 ∇f (E [x(k)]) 2 We denote δ = hp 1 -hαpλ 2 (L) 2 , which is non-negative if h ≤ 2 αλ 2 (L) . Since 2 β(pλn(L)+2-2p) ≤ 2 αλ 2 (L
L 1/2 ∇f (E [x(k)]) 2 ≥ 1 δ (f (E [x(k)]) -f (E [x(k + 1)])) .
If we take the sum over the time instants s we obtain:

k-1 s=0 L 1/2 ∇f (E [x(s)]) 2 ≥ 1 δ (f (x 0 ) -f (E [x(k)])) ,
and for p > 0 the limit satisfies:

lim k→∞ k-1 s=0 L 1/2 ∇f (E [x(s)]) 2 ≥ 1 δ (f (x 0 ) -f * ) . (3.30)
We use (3.30) in (3.28) to get:

lim k→∞ E [f (x(k))] ≤ f (x 0 ) - ξ δ (f (x 0 ) -f * ) , which yields lim k→∞ E [f (x(k))] -f (x * ) ≤ 1 - ξ δ (f (x 0 ) -f * ) .
■ To analyze the asymptotic behavior of the upper bound for J funct , let us consider a step-size given by h = 2ϵ β(pλn(L)+2-2p) with ϵ ∈ (0, 1). For low values of probabilities corresponding to non-robust networks we have: which implies that small step-sizes, characterized by small values of ϵ, reduce more the gap between the expected value of the cost function and f * . In the case of large values of probabilities, we have the following asymptotic behavior:

lim p→0 J funct ≤ ϵ (f (x 0 ) -f * ) , (3.31)
lim p→1 J funct ≤ ϵ(1 -κ -1 κ -1 λ ) 1 -ϵκ -1 κ -1 λ (f (x 0 ) -f * ) , (3.32)
where κ = β α is the condition number of the cost functions f i and κ λ = λn(L) λ 2 (L) . We can observe that only for the particular case of a complete graph (i.e., λ 2 (L) = λ n (L)) and the particular value κ = 1 (i.e., α = β) we obtain

lim p→1 J funct = 0, (3.33)
which implies a certain level of conservatism in the upper bound since (3.33) should hold for any choice of κ and κ λ in order to recover the behavior of a network without packet losses. Figure 3.2 shows the computation of the upper bound for J funct for a complete, star and line graph with n = 10 agents, ϵ = 0.5, f (x 0 ) -f * = 1 and two values of κ where we can appreciate the behavior of the upper bound and the limit values defined by (3.31) and (3.32).

Packet losses in OMAS

In this subsection, we consider the case when replacements of agents happen in the system, making it open. Each agent i ∈ V gets replaced at random time instants, resulting in the change of its local cost function, and hence of the global minimizer of the system x * . For the sake of simplicity, we assume that when an agent i is replaced, the joining agent that takes its place retrieves its label i and its estimate x i , so that the constraint 1 T x = b is preserved, but receives a new local cost function satisfying Assumptions 3.5. Denoting f k i the local cost function held by agent i at the time instant k, we can then reformulate (3.1) as

min x∈S 1,b f k (x) := n i=1 f k i (x i ). (3.34)
The solution of (3.34) can thus differ from a time instant t k to another, and we denote it by x * (k) := arg min x∈S 1,b f k (x).

Discrete-event modelling

The evolution of the open network presented in the previous section is characterized by the instantaneous occurrence at random time instants of either interactions of agents according to the network topology or replacements. We call such occurrences events, and we define the event set of the system, which contains all events that can happen in the system, as

Ξ = U ∪ R, (3.35)
where U is the update event corresponding to the application of the weighted gradient descent, and R := i∈V R i is the set of all events R i , i.e., the replacement of a single agent i in the system. We assume that two distinct events never happen simultaneously, and following the above characterization, we consider a discrete evolution of the time where each time-step k corresponds to the random time instant at which the k-th event of the system ξ k ∈ Ξ takes place. Moreover, we consider the following assumption that guarantees that replacements and interactions are independent processes, so that at any iteration k, the event ξ k is an interaction, i.e., ξ k ∈ U , with fixed probability p U , and a replacement, i.e., ξ k ∈ R, with fixed probability p R = 1 -p U . Assumption 3.6 An event ξ k is independent of any other event ξ j , with j ̸ = k, and of any information in the system prior to time-step k, such as the estimates or local cost functions.

More precisely, we assume that at each time instant k we can have an update with probability p U or a replacement with probability p R = 1 -p U . If an update occurs, the weighted gradient descent is applied and f k+1 = f k . If a replacement occurs, a single agent i is uniformly randomly selected and receives a new (piecewise quadratic) cost function.

Constraint violation

When the open system does not suffer from packet losses, the constraint is always preserved but the global minimizer x * (k) is time-varying since the cost functions change in time. In this analysis, we assume that the open system also suffers of packet losses such that the weighted gradient descent algorithm becomes:

x(k + 1) = x(k) -hL(k)∇f k (x(k)). ( 3
J constr ≥ 4p U h(1 -p) 2 -αh(pλ 2 (L) + 2 -2p) lim k→∞ k-1 s=0 E [f s (x(s)) -f s (x(s + 1))] .
Proof: Let us recall the notation in the proof of Theorem 3.1 for the error (3.22) we have that when an update U happens

H(k) = 1 ⊤ x(k) -b and the natural filtration {F k } k∈Z ≥0 . From
E H 2 (k + 1)|U = E H 2 (k) + 2h 2 p(1 -p)E L 1/2 ∇f k (x(k)) 2 . (3.37)
When a replacement R occurs, the constraint is not violated at that event

E H 2 (k + 1)|R = E H 2 (k) .
Then, we compute the total expectation and we obtain:

E H 2 (k + 1) = p U E H 2 (k + 1)|U + p R E H 2 (k + 1)|R = E H 2 (k) + 2p U h 2 p(1 -p)E L 1/2 ∇f k (x(k)) 2 .
We use a recursive argument to get:

E H 2 (k) = 2p U h 2 p(1 -p) k-1 s=0 E L 1/2 ∇f s (x(s)) 2 . (3.38)
Since the cost functions f k are α-strongly convex we have:

f k (x(k + 1)) ≥ f k (x(k)) + α 2 ||x(k) -x(k + 1)|| 2 + ⟨∇f k (x(k)), x(k + 1) -x(k)⟩ = f k (x(k)) + h 2 α 2 ⟨∇f k (x(k)), L(k) T L(k)∇f k (x(k))⟩ -h⟨∇f k (x(k)), L(k)∇f k (x(k))⟩.
We compute the expectation given the filtration generated by x(k) and we use (3.19) to get:

E f k (x(k + 1))|F k ≥ f k (x(k))+ h 2 α 2 ⟨∇f k (x(k)), (p 2 L 2 + 2p(1 -p)L)∇f k (x(k))⟩ -h⟨∇f k (x(k)), pL∇f k (x(k))⟩.
By using (3.29), we have the lower bound L∇f k (x(k))

2 ≥λ 2 (L) L 1/2 ∇f k (x(k)) 2
and we obtain:

E f k (x(k + 1))|F k ≥ f k (x(k))-hp 1-hα(1 -p)- hαpλ 2 (L) 2 L 1/2 ∇f k (x(k)) 2 .
We denote

θ = hp 1 -hα(1 -p) -hαpλ 2 (L) 2
and we compute the total expectation to get:

E E f k (x(k + 1))|F k ≥ E f k (x(k)) -θE L 1/2 ∇f k (x(k)) 2 . Since h ≤ 2 β(pλn(L)+2-2p) ≤ 2 α(pλ 2 (L)+2-2p) we guarantee that E f k (x(k + 1)) ≤ E f k (x(k)) such that: E L 1/2 ∇f k (x(k)) 2 ≥ 1 θ E f k (x(k)) -E f k (x(k + 1)) . (3.39)
By using (3.39) in (3.38) we complete the proof. ■ Unlike a closed system, we observe that in general E [f s (x(s + 1))] ̸ = E [f s+1 (x(s + 1))] due to the possible replacements that take place in the system. This mismatch implies that terms are not canceled in the series in the lower bound and the error accumulates in time, so that J constr grows to infinity. This intuition is corroborated by the simulations of E (1 T x(k) -b) 2 for a complete graph with n = 7 agents that are shown in Fig. 3.3. We can observe in the plot that both the constraint violation metric and the lower bound diverge. During the transient state, the accumulation of errors is considerable and the increase of the constraint violation metric is significant, while in the steady state the accumulation of errors is small and the metric increases with a lower rate. This difference in the growth rates can be explained by the distance between the initial state x(0) and the minimizer x * . During the transient state, agents exchange large values to reach the minimizer, such that a loss of information has a high impact in the system. Nevertheless, in the steady state, the states of the agents are close to the minimizer and the values exchanged are small, which implies that the potential loss of information is not significant. When replacements do not take place, the lower bound of Proposition 3.2 is given by:

J constr ≥ 4h(1 -p) 2 -αh(pλ 2 (L) + 2 -2p) f (x(0)) -lim k→∞ E [f (x(k))] .
In Fig. 3.4 we present the behavior of J constr and the lower bound of Proposition 3.36 for a system without replacements considering the same parameters of Fig. 3.3 where we observe a similar behavior for the transient state characterized by a fast increase. However, unlike the case of replacements, during the steady state, J constr converges to a value and remains bounded, which is also guaranteed by Theorem 3.1. When p = 1, corresponding to an ideal system without packet losses, even if replacements take place in the system (i.e., p U > 0), J constr is zero since

E [H 2 (k + 1)|U ] = E [H 2 (k)]
as we can see in (3.37). This can also be derived from (3.14), since the algorithm preserves the constraint independently of the cost functions considered for the computation of x(k + 1).

Conclusion

In this chapter, we analyzed the effect of packet losses in the performance of the weighted gradient descent algorithm. We defined two performance metrics to measure the deviation from the constraint and the error of the expected cost function and we derived an upper bound. We extended the analysis of packet losses to OMAS and we showed that the constraint violation metric may diverge in such a case due to the accumulation of errors during the replacements.

We showed that when the weighted gradient descent algorithm is applied in a system subject only to packet losses, the constraint violation metric remains bounded while in a system subject only to replacements, the constraint violation metric is zero since the constraint is preserved during replacements. However, when updates using the weighted gradient descent algorithm with packet losses and replacements occur in the same system, the constraint violation metric diverges, showing the fragility of the algorithm in OMAS. This result presents a similarity with the case of Lévy flights studied in [START_REF] Wang | Distributed averaging under constraints on information exchange: emergence of Lévy flights[END_REF], where an asymmetric loss of information together with delays and noise can generate complex behaviors.

As for future work, natural continuations of the paper in closed and open systems can be explored. In the closed case, although an upper bound was derived for the cost function metric, the expression presents a certain level of conservatism. It would be important to improve the bound to obtain a tighter result that guarantees a zero value of the metric when p = 1 (no replacements). Regarding the open case, it is necessary to improve the lower bound for the constraint violation metric so that it does not depend on the trajectories but on the properties of the system. Also, the error of the expected cost function should be studied in OMAS to determine if it is bounded or if it could diverge.

Finally, regarding the objective of achieving an optimal resource allocation even under the perturbations generated by packet losses and replacements of agents, the design of more complex algorithms to obtain better performances remains open. One of the most interesting approaches is online optimization, where the problem of time-varying functions is also analyzed from a different perspective and algorithms based on several stages are proposed with the aim of achieving a good performance in average [START_REF] Carnevale | GTAdam: Gradient tracking with adaptive momentum for distributed online optimization[END_REF][START_REF] Carnevale | Distributed online aggregative optimization for dynamic multi-robot coordination[END_REF].

Chapter 4

Random coordinate descent with replacements

In OMAS, a fixed solution for the resource allocation problem (3.1) cannot be obtained as in general, the size of the system is not fixed and the cost functions keep changing, such that the goal of the agents is to track the time-varying solution of (3.1) as well as possible at all times using an optimization algorithm. However, as the size of such systems reaches large values, global optimization methods relying, e.g., on the computation of the whole gradient of f are not suited since the computational complexity would be high and in some cases, it would not be practical to gather the whole gradient as agents may have entered/left in the meantime.

For this reason, it is important to consider optimization algorithms based on local interactions, since they are more flexible. An alternative type of algorithms that allow to considerably reduce the computational complexity is the so-called Coordinate Descent algorithm introduced by Nesterov, where the optimization is performed only along one direction at each iteration [START_REF] Nesterov | Efficiency of coordinate descent methods on huge-scale optimization problems[END_REF]. For multi-agent systems, the selection of one coordinate is equivalent to the choice of a particular edge of the network to perform the optimization. In such algorithms, the sequence of edges is crucial, and hence a randomized choice denoted as Random Coordinate Descent algorithm (RCD) was studied in [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF], where convergence of the cost functions is proved under standard assumptions when only pairwise interactions are considered, so that the algorithm requires only the computation of a pair of local gradients per iteration.

In the context of optimization, problems such as (3.1) typically imply a minimization process on a long period, and thus in many situations, the cost is expected to be paid on a regular basis. In such a setting, a natural way of measuring the performance of an algorithm is to compute its accumulated error with respect to a given strategy over a finite number of iterations. Similar metrics occur in the context of online optimization [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF], where the objective is to minimize the so-called 51 regret, commonly defined as the accumulated error of the estimate x(k) with respect to x * := arg min x T k=1 f k (x), or sometimes with respect to the time-varying solution of (3.1) x * (k) := arg min x f k (x) such as e.g., in [START_REF] Shahrampour | Distributed online optimization in dynamic environments using mirror descent[END_REF].

In this chapter, we study the Random Coordinate Descent algorithm (RCD) to solve the resource allocation (RA) problem in OMAS where agents get replaced during the process, relying on a decoupled analysis of the RCD algorithm and of replacements of single agents. In the first part, we analyze the RCD algorithm using a time-varying optimization framework [START_REF] Simonetto | Time-varying convex optimization: Time-structured algorithms and applications[END_REF], where we focus in the convergence of the algorithm towards the minimizer and its stability. In Section 4.1 we introduce a more general setting for the RA problem where agents can be heterogeneous (i.e., that a i may be different) and hold d-dimensional local cost functions. In addition, we consider an arbitrary distribution for the probabilities associated with the choice of edges during the execution of the RCD algorithm. Section 4.2 presents an analysis of the behavior of the minimizer during a replacement. We derive three different upper bounds for the difference between two possible minimizers using different techniques. In Section 4.3 we prove the linear convergence of the RCD algorithm in a closed system in a norm induced by a matrix associated with the network. In Section 4.4 we analyze the RCD algorithm in an open system subject to replacements and we provide conditions for the convergence of the algorithm inside a ball, whose size depends on parameters of the system. In the second part, we analyze the RCD algorithm using tools similar to those used in Online Optimization. Section 4.5 presents a simplified setting, where homogeneous agents, holding 1-dimensional cost functions, interact in a complete graph. We introduce the performance metrics used in this work, which are inspired by the loss accumulated by the RCD algorithm with respect to the time-varying optimal solution x * (k) over a finite number of iterations, and its gain with respect to the selfish strategy x s (k), which consists in the absence of collaboration between the agents (i.e., x s i (k) = d k i ). We derive upper bounds for the performance metrics and show that the metrics scale linearly with time. In addition, the case of quadratic cost functions is analyzed, for which tighter results are obtained. Finally, the conclusions and future work are exposed in Section 4.6.

Problem statement

Network description and open system

We consider the resource allocation problem (3.1) with a fixed budget b:

min x=[x T 1 ,...,x T n ] T ∈R nd f (x) = n i=1 f i (x i ) subject to n i=1 a i x i = b, (4.1)
where the local cost functions f i : R d → R satisfy Assumption 3.1 (i.e., f i is continuously differentiable, α-strongly convex and β-smooth).

In addition, we assume to have an undirected and connected graph G = (V, E) where the set of nodes is given by V = {1, . . . , n} and the set of edges by E ⊆ V ×V . Each agent i ∈ V has access to a local cost function f i and to a local variable x i . Agents can exchange information at random times through pairwise interactions according to the network G. Whenever an interaction happens in the system, an edge (i, j) ∈ E is selected with some fixed probability p ij > 0 and agents i and j can then exchange information in an undirected manner to update their respective estimates.

Moreover, we consider that replacements of agents happen in the system, making it open. Each agent i ∈ V gets replaced at random time instants, resulting in the change of its local cost function, and hence of the global minimizer x * .

The evolution of the open network presented in this section is characterized by the instantaneous occurrence at random time instants of either pairwise interactions or replacements. We call such occurrences events, and we define the event set of the system, which contains all events that can happen in the system, as

Ξ = U ∪ R, (4.2)
where U = (i,j)∈E U ij is the set of all possible events U ij , i.e., the pairwise interaction between two connected agents i and j, and R := i∈V R i is the set of all events R i , i.e., the replacement of a single agent i in the system. Similarly to Section 3.3.1, we consider a discrete evolution of the time where each time-step k corresponds to the random time instant at which the k-th event of the system ξ k ∈ Ξ takes place. Moreover, we consider that the events ξ k satisfy the Assumption 3.6 that guarantees that replacements and interactions are independent processes, so that at any iteration k, the event ξ k is an interaction, i.e., ξ k ∈ U , with fixed probability p U , and a replacement, i.e., ξ k ∈ R, with fixed probability p R = 1 -p U .

Following the approach in [START_REF] Hendrickx | Stability of decentralized gradient descent in open multi-agent systems[END_REF], we restrict the location of the minimizers of the local functions:

Assumption 4.1 For all i ∈ V , the minimizer of f i denoted as x * i := arg min x f i (x) satisfies x * i ∈ B(0 d , c
) for some c > 0. Moreover, without loss of generality

f i (x * i ) = 0 for all i ∈ V .
Assumption 4.1 guarantees a certain level of uniformity among the local cost functions. In particular, it prevents arbitrary changes of functions, and thus of x * , during replacements.

In the first part of this chapter, we assume that when agent i is replaced, the joining agent that takes its place retrieves its label i and its estimate x i , so that the constraint a ⊤ ⊗ I d x = b is preserved, but receives a new local cost function satisfying Assumptions 3.1 and 4.1 1 . This is consistent with e.g., the case of energy distribution where the sources maintain the same value of generated energy but with a different cost function. Denoting f k i the local cost function held by agent i at the time instant k, we can then reformulate the RA problem (4.1) as

min x∈S a,b f k (x) := n i=1 f k i (x i ), (4.3)
where

S a,b = x ∈ R nd | a ⊤ ⊗ I d x = b
is the feasible set. The solution of (4.3) can thus differ from a time instant t k to another, and we denote it by x * (k) := arg min x∈S a,b f k (x). The objective of the agents is to track x * (k) as well as possible even though replacements happen in the system.

Random Coordinate Descent (RCD) algorithm

To compute the solution of (4.3) we consider the Random Coordinate Descent (RCD) algorithm introduced in [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF]. This algorithm involves the update of the states of only a pair of neighbouring agents at each iteration, so that the computational complexity is cheap. Hence, whenever two agents i and j interact (i.e., in the event U ij ), they perform an RCD update, which is defined as follows for some nonnegative step-size h ≥ 0:

x(k + 1) = x(k) -hQ ij ∇f (x(k)), (4.4)
where Q ij is the nd × nd matrix defined as

Q ij = Q ij ⊗ I d , with Q ij the n × n
matrix filled with zeroes except for the four following entries:

[Q ij ] i,i = a 2 j a 2 i + a 2 j ; [Q ij ] j,j = a 2 i a 2 i + a 2 j ; [Q ij ] i,j = - a i a j a 2 i + a 2 j ; [Q ij ] j,i = - a i a j a 2 i + a 2 j .
With the update rule (4.4) only agents i and j update their estimates while all the other agents keep it the same. For agents i and j, (4.4) essentially amounts to perform a gradient step on the function f i (x i ) + f j (x j ) under the constraint that a i x i + a j x j remains constant. This ensures that the resource allocation constraint is preserved as long as the starting point satisfies it. In particular, in the case of homogeneous agents (i.e., where a = 1 n ), then one shows that

x + i = x i - h 2 (∇f i (x i ) -∇f j (x j ))
, so that the update follows both gradients with equal weight while preserving the constraint.

Remark 4.1 The update rule (4.4) can be formally obtained by solving the following optimization problem, which corresponds to the interpretation given above (we refer to [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF] for details):

arg min

s i ,s j ∈R d :a i s i +a j s j =0 ∇f i (x i ) ∇f j (x j ) , s i s j + β 2
s i s j 2 . (4.5)
Based on the approach of [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF], one can then show that the optimal step-size that solves (4.5) is given by h = 1/β.

We also introduce the following matrix that builds on the definition of the update rule (4.4):

L p = (i,j)∈E p ij Q ij =   (i,j)∈E p ij Q ij   ⊗ I d = L p ⊗ I d . (4.6)
This matrix appears in the dynamics corresponding to the conditional expectation:

E [x(k + 1)|x(k)] = x(k) -hL p ∇f (x(k)),
and will be used for the definition of an appropriate norm for the analysis of the RCD algorithm. Observe that by definition of Q ij and L p , one has

L p a = Q ij a = 0 n , (4.7)
which means that zero is an eigenvalue of both Q ij and L p with corresponding eigenvector a. We denote by λ 2 (L p ) and λ n (L p ) respectively the second smallest and the largest eigenvalues of L p . One can show that λ 2 (L p ) > 0 since the graph G is connected (we refer to Lemma 3.3 of [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF] for a detailed proof).

Remark 4.2 For a graph G = (V, E), when a = 1 n (homogeneous agents) and the probabilities p ij are uniformly distributed, one has L p = 1 2|E| L, where L is the usual Laplacian of the graph. Hence, we refer to L p as a scaled Laplacian, as it enjoys similar properties, especially in terms of eigenvalues.

Effect of replacements

In this section, we bound the distance by which the minimizer of f can change after the replacement of a single agent, i.e., the modification of a single function. Our first two results concern the location of the minimizer: Lemma 4.1 corresponds to an analysis of the location of the global minimizer of the system, and Lemma 4.2 is an analysis of the location of the minimizer held by each individual agent.

f (x) ≥ α 2 ||x -x * || 2 > α 2 κ √ nc + ||b|| ||a|| 2 = β 2 √ nc + ||b|| ||a||
f (x b ) ≤ β 2 ||x b -x * || 2 ≤ β 2 √ nc + ||b|| ||a|| 2 .
Finally, since x b ∈ S a,b there holds f * ≤ f (x b ) by definition. Combining all the inequalities above then yields

f (x) > β 2 √ nc + ||b|| ||a|| 2 ≥ f (x b ) ≥ f * ,
which implies that x cannot be the minimizer of the problem and concludes the proof. Proof: Let us denote x * i the minimizer of f i without constraint which satisfies f i (x * i ) = 0 and ∇f i (x * i ) = 0 d . From β-smoothness of the local functions we have:

■ Lemma 4.2 Let x * := arg min x∈S a,b f (x) = n i=1 f i (x i ). If f i satisfies Assump- tions 3.
||∇f i (x * i )|| 2 ≤ β⟨∇f i (x * i ), x * i -x * i ⟩.
Then, from the optimality condition (3.3), there holds

a 2 i ||ζ * || 2 ≤ β⟨ζ * , a i (x * i -x * i )⟩
. By summing over all the i, we obtain:

||a|| 2 ||ζ * || 2 ≤ β⟨ζ * , n i=1 a i (x * i -x * i )⟩.
We use the Cauchy-Schwarz inequality to get:

||a|| 2 ||ζ * || 2 ≤ β ||ζ * || n i=1 a i (x * i -x * i ) ,
and by using the triangle inequality and the fact that n i=1 a i x * i = b we obtain:

||ζ * || ≤ β ||a|| 2 ||b|| + n i=1 a i x * i ≤ β ||a|| 2 (||b|| + c ||a|| 1 ) ,
which corresponds to (4.9). From α-strongly convexity of the local functions we have:

α ||x i -x * i || ≤ ||∇f (x i )||
. By using the reverse triangle inequality and the optimality condition we get

||x * i || ≤ 1 α ||a i ζ * || + ||x * i || ≤ a i α ||ζ * || + c. (4.11)
Equation (4.10) then follows from combining (4.9) and (4.11). ■ We can now use these two lemmas to characterize the evolution of the distance between the estimate x(k) and the minimizer x * (k) after a replacement event. Without loss of generality, we assume that agent n, and hence f n , is replaced, and for the n

+ 1 functions f 1 , f 2 , . . . , f n-1 , f (1) n , f (2)
n satisfying Assumptions 3.1 and 4.1 we define the minimizer before a replacement x (1) , and after a replacement x (2) as

x (1) (1) n (x n ) ;

:= arg min x∈S a,b n-1 i=1 f i (x i ) + f
x (2) := arg min

x∈S a,b n-1 i=1 f i (x i ) + f (2) n (x n ) . (4.12)
Proposition 4.1 Consider x (1) and x (2) as defined in (4.12), let a + and a -respectively denote the largest and smallest values in a, and let ρ a :=

a 2 + ||a|| 2 -a 2 +
, then there holds

x (1) -x (2) 2 ≤ min{ψ n,κ , χ n,κ , θ n,κ } =: M 2 n,κ , (4.13)
with

ψ n,κ = 4nκ c + c √ κ + ||b|| √ n ||a|| 2 ; (4.14) χ n,κ = 8 a + ||a|| 2 κ(||b|| + c ||a|| 1 ) + c 2 ;
(4.15)

θ n,κ = 4 1 + (κ + 1) 2 4κ ρ a a + ||a|| 2 κ(||b|| + c ||a|| 1 ) + c 2 , (4.16)
Proof: a) Proof of ψ n,κ : Observe that x (1) , x (2) ∈ B(0 nd , R b,κ ) from Lemma 4.1, so that there holds

x (1) -x (2) 2 ≤ 2 x (1) 2 + x (2) 2 ≤ 4R 2 b,κ ,
which yields (4.14).

b) Proof of χ n,κ : We remind that for i = 1, . . . , n -1, there holds ∇f i (x q) , with q = 1, 2. From α-strongly convexity of the local cost functions, there holds that for all i = 1, . . . , n -1:

(q) i ) = a i ζ (
a i ⟨ζ (1) -ζ (2) , x (1) i -x (2) i ⟩ ≥ α x (1) i -x (2) i 2 .
Let us define y (q) ∈ R d(n-1) the vector such that y (q) i = x (q) i for q = 1, 2 and for i = 1, . . . , n -1. Using the fact that n i=1 a i x (q) i = b for q = 1, 2 and summing up the above inequalities over all i = 1, . . . , n -1 yields a n ⟨ζ (1) -ζ (2) , x (2) n -x (1) n ⟩ ≥ α y (1) -

y (2) 2 ,
where we used the fact that n-1 i=1 a i x

(q) i + a n x (q)
n = b. By using Cauchy-Schwarz inequality and triangle inequality we obtain

y (1) -y (2) 2 ≤ a n α ζ (1) + ζ (2)
x (1) n + x (2) n Then, we use (4.9) and (4.10) to get

y (1) -y (2) 2 ≤ 4 a n ||a|| 2 κ (||b|| + c ||a|| 1 ) + c .
(4.17)

By definition we have

x (1) -x (2) 2 = y (1) -y (2) 2 + x (1) n -x (2) n 2 . (4.18)
We apply triangle inequality and (4.10) to obtain

x (1) -x (2) 2 ≤ y (1) -y (2) 2 + x (1) n + x (2) n 2 ≤ y (1) -y (2) 2 + 2 x (1) n 2 + x (2) n 2 ≤ y (1) -y (2) 2 + 4 a n ||a|| 2 κ (||b|| + c ||a|| 1 ) + c (4.19)
Finally, the result (4.15) yields by combining (4.17) and (4.19) and using the fact that a n ≤ a + .

c) Proof of θ n,κ : Since the functions are α-strongly convex and β-smooth, there holds that for all i = 1, . . . , n -1:

a i (1 + κ -1 )⟨ζ (1) -ζ (2) , x (1)
i -x (2) i ⟩ ≥ β -1 a 2 i ζ (1) -ζ (2) 2 + α x (1) i -x (2) i 2
.

By summing up the above inequalities over all i = 1, . . . , n -1 yields

a n (1 + κ -1 )⟨ζ (1) -ζ (2) , x (2) n -x (1) n ⟩ ≥ mβ -1 ζ (1) -ζ (2) 2 + α y (1) -y (2) 2 ,
where m = n-1 i=1 a 2 i . By using Cauchy-Schwarz inequality we obtain:

a n (1 + κ -1 ) ζ (1) -ζ (2) x (2) n -x (1) n ≥ mβ -1 ζ (1) -ζ (2) 2 +α y (1) -y (2) 2 .
This can be written as follows

α y (1) -y (2) 2 ≤ ϕ( ζ (1) -ζ (2) ), (4.20)
where

ϕ(z) = -mβ -1 z 2 + a n (1 + κ -1 ) x (2) n -x (1) n z.
Since ϕ is a concave parabola, there exists ϕ * = max z ϕ(z) < ∞ such that ϕ(z) ≤ ϕ * for all z given by

ϕ * = a 2 n (1 + κ -1 ) 2 x (2) n -x (1) n 2 4mβ -1 ,
and it follows by using (4.20) that

y (1) -y (2) 2 ≤ a 2 n (1 + κ -1 ) 2 αβ -1 x (1) n -x (2) n 2 4m = √ κ + 1 √ κ 2 a 2 + x (1) n -x (2) n 2 4 ||a|| 2 -a 2 + .
Equation (4.16) then follows from (4.18) using the fact that x

n and x

n are upper bounded by (4.10), which thus yields θ n,κ , and the conclusion follows.

■ The bound M 2 n,κ from Proposition 4.1 is obtained by taking the minimum between three quantities: ψ n,κ , χ n,κ and θ n,κ . The first one follows from the largest possible distance existing between two minimizers, defined by the region in which they can be located. The second and third ones rely on the largest possible distance between the local minimizers corresponding to the replaced agents. While χ n,κ and θ n,κ are derived using inequalities associated with α-strongly convex functions, the proof of θ n,κ also involves the use of additional properties corresponding to β-smooth functions and the determination of the maximum value of a concave function. The bound θ n,κ shows a strong dependence on the weights of the agents through the coefficient ρ a , which is not present in the other two bounds. Notice that the bounds χ n,κ and θ n,κ coincide when

(κ + 1) 2 4κ ρ a = 1.
Let ā and a 2 respectively stand for the average value and average of the squared values of a. One can more generally highlight the dependencies of the three quantities with the parameters using standard algebraic manipulations, yielding

ψ n,κ ≤ 4nκ 2c + ||b|| nā 2 = O(nκ); χ n,κ ≤ 8 a 2 + a 2 ||b|| nā + 2c κ 2 = O(κ 2 ); θ n,κ ≤ 4 a 2 + a 2 - κ 2(n-1) + 2 a 2 + a 2 ||b|| nā + 2c κ 2 = O κ 2 + κ 3 n ,
where O(•) corresponds to the O-notation 2 . The linear scaling of ψ n,κ in both n

2 For a given function g(n), we denote by O(g(n)) the set of functions [START_REF] Cormen | Introduction to Algorithms[END_REF] O(g(n)) = {f (n) : there exist positive constants c and n 0 such that

0 ≤ f (n) ≤ cg(n) for all n ≥ n 0 }.
and κ and the higher order scaling of both χ n,κ and θ n,κ in only κ suggest that ψ n,κ is tighter for small values of n and large values of κ, whereas θ n,κ and χ n,κ are tighter otherwise. The main difference between χ n,κ and θ n,κ lies in a multiplicative factor, constant for the former, and depending of the parameters and the values in a for the latter. In general, χ n,κ tends to be tighter than θ n,κ as κ gets large and n small. This difference becomes significant in heterogeneous settings, where it can get tighter than ψ n,κ as well. These behaviors are illustrated in Fig. 4.1. Remark 4.3 The interpretation of the quantities ψ n,κ , χ n,κ and θ n,κ actually strongly depends on the implicit assumption that ||b|| is fixed and ||a|| 1 scales with n (i.e., ā is fixed). This particular modelling choice is arbitrary, and implies that the solution held by an agent x * i (k) becomes smaller for large values of n. Other choices might have different implications on the interpretation, and in particular on the scaling of these quantities. For instance one could choose to either fix ||b|| and ||a|| 1 , or that both ||b|| and ||a|| 1 scale with n, so that the x * i (k) remain mostly the same no matter n (observe that the latter yields the same scalings than those presented above).

The result of Proposition 4.1 can be analyzed with respect to empirical results derived with the PESTO toolbox [START_REF] Taylor | Performance estimation toolbox (PESTO): Automated worst-case analysis of first-order optimization methods[END_REF], which allows computing exact empirical bounds for questions related to convex functions 3 . For the sake of simplicity, the analysis here is only done for the homogeneous case, and consequently does not involve χ n,κ , similar conclusions could however be drawn the same way e.g., using heterogeneous agents.

We can observe in Fig. 4.2, that even though there is some gap between the theoretical result and that obtained using PESTO, the scaling of the bounds with respect to n and κ is well captured. In particular, the top plot shows that M 2 n,κ = θ n,κ when n becomes large, resulting in the convergence of M 2 n,κ towards a constant, consistently with the result obtained with PESTO. In parallel, the bottom plot shows that the bounds from PESTO scale linearly with κ, consistently with the evolution of ψ n,κ , which is the value taken by M 2 n,κ for large values of κ.

Linear convergence of RCD in closed system

We now analyze the effect of the second type of events happening in the system, i.e., pairwise interactions resulting in RCD updates. This corresponds to studying the convergence of the RCD Algorithm in closed system.

Linear convergence and L † p -seminorm

In this section, we derive the convergence rate of the RCD algorithm in terms of the distance to the minimizer with the objective of characterizing the effect of a single RCD step on that expected distance at interaction events. We introduce the following standard definitions [START_REF] Ortega | Iterative Solution of Nonlinear Equations in Several Variables[END_REF].

Definition 4.1 (Linear Convergence) Let {x(k)} be the sequence of points converging to some point x * ∈ R d generated by some algorithm. We say the convergence is linear if there exists r ∈ (0, 1) such that for all k

||x(k + 1) -x * || ≤ r ||x(k) -x * || .
The number r is called the rate of convergence. Definition 4.2 (Exponential Convergence) Let {x(k)} be the sequence of points converging to some point x * ∈ R d generated by some algorithm. We say 3 We consider the setting of Proposition 4.1 with arg min x f i (x) ∈ B(0, 1). We use PESTO to evaluate max x (2) -x (1) 2 , where x (1) and x (2) are defined as in (4.12) and where we impose

i x (1) i = i x (2) i = v b , for some vector v b satisfying ||v b || = b.
The variables of the problem are the functions f i , the vector v b and the decision variables x i , and PESTO derives the performance achieved by the empirical worst-case instance of the problem. the convergence is exponential if there exists r ∈ (0, 1) and some positive constant

C such that for all k ||x(k) -x * || ≤ Cr k .
In optimization, exponential convergence is also known as R-linear convergence.

Clearly, exponential convergence is weaker than linear convergence since it is concerned with the overall rate of decrease in the error, rather that the decrease over each individual iteration of the algorithm [START_REF] Nocedal | Numerical Optimization[END_REF].

In [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF], the author proves linear convergence of the RCD algorithm in expec-tation in terms of the function value, i.e., f (x) -f * . Hence, from the inequalities corresponding to smooth functions and strong convexity [START_REF] Nesterov | Lectures on Convex Optimization[END_REF][START_REF] Bubeck | Convex Optimization: Algorithms and Complexity[END_REF], it is straightforward to prove exponential convergence of the algorithm from [80, Eq. ( 26)]:

E [||x(k) -x * ||] ≤ κ(1 -αλ 2 (L p )) k ||x 0 -x * || . (4.21)
However, due to the alternation of updates and replacements, our analysis in open systems requires the strict contraction of some metric after each iteration. The linear convergence was established in the preliminary work [START_REF] Monnoyer De Galland | Random coordinate descent algorithm for open multi-agent systems with complete topology and homogeneous agents[END_REF] for the Euclidean norm under the assumption of a complete communication graph with homogeneous agents and uniform probabilities p ij . Nevertheless, the following example shows that such contraction no longer holds for the Euclidean norm for general graphs.

Example 4.1 Consider a line graph with 3 agents satisfying the constraint ⟨1, x⟩ = -3 with probabilities p 12 = 0.9, p 23 = 0.1 (and hence p 13 = 0), and whose local cost functions and initial conditions are:

i f i (x i ) x * i x i (0) 1 50(x 1 -2) 2 2 10 2 20(x 2 + 2) 2 -2 7 3 (x 3 + 3) 2 -3 -20
Starting from x(0) the expected result of the RCD operation with step-size h = 1/β = 0.01 is

E ||x(1) -x * || 2 = 437.204 > 434 = ||x(0) -x * || 2 , (4.22)
and hence linear convergence cannot be achieved at least in this case.

For this reason, we propose to study the problem in a different norm associated with the algorithm. Since the RCD is performed along a network of agents, a natural choice is to consider norms induced by associated matrices as in [START_REF] Wang | Laplacian smoothing stochastic gradient Markov Chain Monte Carlo[END_REF]. In this case, we focus on the seminorm induced by the Moore-Penrose inverse of the matrix L p defined in (4.6), denoted by L † p , which we remind is defined as follows for some x ∈ R nd :

∥x∥ L † p := x ⊤ L † p x. (4.23)
We show with the next proposition that this seminorm is a norm on S a,0 , where we recall that S a,0 is the feasible set defined in (3.2) when b = 0 d and corresponds to the kernel of a T ⊗ I d . For the particular case d = 1, S a,0 is the orthogonal complement of a. Proof: By (4.23), ∥x∥ L † p = 0 implies that x must be in the kernel of L † p , defined as ker(L

† p ) = {x ∈ R nd | x = a ⊗ w, w ∈ R d }.
Since x ∈ S a,0 , it must satisfy (a ⊤ ⊗ I d )x = 0 nd and we have:

a ⊤ ⊗ I d x = a ⊤ ⊗ I d (a ⊗ w) = ||a|| 2 w,
which is equal to 0 d only for w = 0 d . ■ If x, y ∈ S a,b , then z = x -y belongs to S a,0 , so that the norm ∥•∥ L † p can be used to measure the distance between two vectors in the context of this work.

Contraction of an iteration in closed system

Let us remind the update rule of the RCD algorithm defined in (4.4) for some positive step-size h as

x(k + 1) = x(k) -hQ ij ∇f (x(k)).
(4.24)

In the following proposition, we analyze the convergence of (4.24) with respect to the norm induced by L † p defined tn the previous section. and for any initial point x ∈ S a,b , then the update rule (4.24) applied on the randomly selected pair of agents (i, j) ∈ E satisfies

E ∥x(k + 1) -x * ∥ 2 L † p ≤ 1 -2hαλ 2 (L p ) + h 2 α 2 λ n (L p ) ∥x(k) -x * ∥ 2 L † p . (4.

26)

Proof: There holds

E ∥x(k + 1) -x * ∥ 2 L † p = (i,j)∈E p ij ∥x(k) -hQ ij ∇f (x(k)) -x * ∥ 2 L † p = ∥x(k) -x * ∥ 2 L † p + h 2 (i,j)∈E p ij ∥Q ij ∇f (x(k))∥ 2 L † p -2h (i,j)∈E p ij ⟨Q ij ∇f (x(k)), L † p (x(k) -x * )⟩. (4.27) Since L p = (i,j)∈E p ij Q ij , it follows that E ∥x(k + 1) -x * ∥ 2 L † p = ∥x(k) -x * ∥ 2 L † p + h 2 (i,j)∈E p ij ∥Q ij ∇f (x(k))∥ 2 L † p -2h⟨L p ∇f (x(k)), L † p (x(k) -x * )⟩. (4.28)
We first treat the second term of the right-hand side of (4.28). Remember from (4.7) that Q ij a = 0 n , and from

(3.3) that ∇f * = a ⊗ ζ * for some ζ * ∈ R d . Hence, since Q ij = Q ij ⊗ I d by definition, there holds Q ij ∇f * = (Q ij ⊗ I d )(a ⊗ ζ * ) = (Q ij a) ⊗ ζ * = 0 nd . (4.29)
It thus follows that

∥Q ij ∇f (x(k))∥ 2 L † p = ∥Q ij (∇f (x(k)) -∇f * )∥ 2 L † p ≤ 1 λ 2 (L p ) Q ij (∇f (x(k)) -∇f * ) 2 , (4.30)
where the inequality follows from the fact that the eigenvalues of L p are exactly those of L p repeated d times (by Theorem 13.12 of [START_REF] Laub | Matrix Analysis For Scientists And Engineers[END_REF]), so that the smallest and largest nonzero eigenvalues of L † p are respectively 1/λ n (L p ) and 1/λ 2 (L p ), yielding for all z ∈ R nd :

∥z∥ 2 L † p ≤ 1 λ 2 (L p ) ||z|| 2 . (4.31)
Therefore, since

Q ij = (Q ij ) ⊤ = (Q ij ) 2
, and using the fact that ||z|| 2 Lp ≤ λ n (L p ) ||z|| 2 for all z ∈ R nd , there holds from (4.30):

(i,j)∈E p ij ∥Q ij ∇f (x(k))∥ 2 L † p ≤ 1 λ 2 (L p ) ||∇f (x(k)) -∇f * || 2 Lp ≤ λ n (L p ) λ 2 (L p ) ||∇f (x(k)) -∇f * || 2 . (4.32)
We now analyze the third term of the right-hand side of (4.28). From (4.29) there holds

L p ∇f * = (i,j)∈E p ij Q ij ∇f * = 0 nd , (4.33) yielding ⟨∇f (x(k)) -∇f * , x(k) -x * ⟩ ≥ β -1 ||∇f (x(k)) -∇f * || 2 1 + κ -1 + α ||x(k) -x * || 2 1 + κ -1 .
Hence, using the result above and (4.31), it follows that

-2h⟨L p ∇f (x(k)), L † p (x(k) -x * )⟩ ≤ -2h β -1 ||∇f (x(k)) -∇f * || 2 1 + κ -1 -2h αλ 2 (L p )∥x(k) -x * ∥ 2 L † p 1 + κ -1 . ( 4
E ∥x(k + 1) -x * ∥ 2 L † p ≤ 1 -2h αλ 2 (L p ) 1 + κ -1 ∥x(k) -x * ∥ 2 L † p (4.35) + h 2 λ n (L p ) λ 2 (L p ) -2h β -1 1 + κ -1 ||∇f (x(k)) -∇f * || 2 . Observe that if h ≤ λ 2 (Lp) λn(Lp) 2 α+β then h 2 λn(Lp) λ 2 (Lp) -2h β -1
1+κ -1 ≤ 0, and hence we can use the definition of α-strongly convex functions to find that

||∇f (x(k)) -∇f * || 2 ≥ α 2 λ 2 (L p )∥x(k) -x * ∥ 2 L † p ,
where we used the fact that ||x(k) -

x * || 2 ≥ λ 2 (L p )∥x(k) -x * ∥ 2 L † p . It follows that E ∥x(k + 1) -x * ∥ 2 L † p ≤ ∥x(k) -x * ∥ 2 L † p + (h 2 α 2 λ n (L p ) -2hαλ 2 (L p ))∥x(k) -x * ∥ 2 L † p ,
which concludes the proof. ■ It is clear that the rate of convergence is less than one if h ≤ 2λ 2 (Lp) αλn(Lp) , which is thus satisfied on all its range of validity since h ≤ 2λ 2 (Lp) (α+β)λn(Lp) ≤ 2λ 2 (Lp) αλn(Lp) . We can then find the step-size which minimizes (4.26) and the corresponding convergence rate.

Corollary 4.1

The optimal rate of convergence in (4.26) under (4.25) is achieved for h * = 2λ 2 (Lp) (α+β)λn(Lp) which yields

E ∥x(k + 1) -x * ∥ 2 L † p ≤ 1 - λ 2 (L p ) 2 λ n (L p ) 1 κ ∥x(k) -x * ∥ 2 L † p . (4.36)
Interestingly, Proposition 4.3 shows that linear convergence can be achieved by the RCD algorithm with respect to the norm induced by L † p with a convergence rate similar to that of classical algorithms based on the gradient descent [START_REF] Nesterov | Lectures on Convex Optimization[END_REF][START_REF] Bubeck | Convex Optimization: Algorithms and Complexity[END_REF].

Remark 4.4 (Complete graph) For the particular case of a complete graph with 1-dimensional homogeneous agents and uniform probabilities, the eigenvalues of L p are λ 2 (L p ) = λ n (L p ) = 1 n-1 and the L † p -norm coincides with the Euclidean norm for all z = x -y, where x, y ∈ S a,b . Then, the result of Proposition 4.3 becomes

E ||x(k + 1) -x * || 2 ≤ 1 - αh n -1 (2 -αh) ||x(k) -x * || 2 .
Since in that case by definition h ≤ 2 α+β ≤ 1 α , it follows that

E ||x(k + 1) -x * || 2 ≤ 1 - αh n -1 ||x(k) -x * || 2 ,
which coincides with [94, Eq. ( 13)].

λn(Lp) 2 α+β ≤ h ≤ κ -1 +κ L κ -1 +1 λ 2 (Lp) βλ 2 n (Lp) , with κ L = λn(Lp) λ 2 (Lp) : E ∥x(k + 1) -x * ∥ 2 L † p ≤ 1 -2βλ 2 (L p )h + h 2 β 2 λ 2 n (L p ) λ 2 (L p ) ∥x(k) -x * ∥ 2 L † p . (4.37)
This result could be used in the rest of this chapter the same way as that of Proposition 4.3 for the corresponding step-size. This development is however omitted as it is only marginal progress.

Homogeneous agents and uniform probabilities

For the particular case of homogeneous agents (i.e., a = 1 n ) and uniform probabilities p ij = p, the matrix L p can be expressed as L p = p 2 L where L is the usual Laplacian matrix. In this case, the matrix L p can be associated to an electrical circuit [START_REF] Dörfler | Electrical networks and algebraic graph theory: Models, properties, and applications[END_REF], and we can use the concept of effective resistance to find an upper bound for the step-size of the algorithm independently of λ 2 (L p ).

Hence, the following proposition provides an alternative bound for the convergence of the RCD algorithm in the specific case described above, and can be used the same way as that of Proposition 4.3 in the remainder of this chapter for that case. However, for the sake of generality, we express the main result in the next section only in terms of Proposition 4.3.

h ≤ 2p λ n (L p ) 2 α + β ,
and for any initial point x ∈ S a,b , then the update rule (4.24) applied on the randomly selected pair of agents (i, j) ∈ E satisfies

E ∥x(k + 1) -x * ∥ 2 L † p ≤ 1 -2hαλ 2 (L p ) + h 2 α 2 λ 2 (L p )λ n (L p ) 2p ∥x(k) -x * ∥ 2 L † p . (4.38)
Proof: Since the matrices Q ij are idempotent, the summation term of the second element in (4.28) can be expressed as:

∥Q ij ∇f (x(k))∥ 2 L † p = ⟨Q ij ∇f (x(k)), Q ij L † p Q ij Q ij ∇f (x(k))⟩
Then we can use an upper bound for the quadratic form and we obtain for each term:

∥Q ij ∇f (x(k))∥ 2 L † p ≤ Q ij ∇f (x(k)) 2 λ max (Q ij L † p Q ij ). Now, the matrix Q ij L † p Q ij is
given by:

Q ij L † p Q ij = Q ij L † p Q ij ⊗ I d , which implies that λ max (Q ij L † p Q ij ) = λ max (Q ij L † p Q ij )
. Then, we have:

Q ij L † Q ij = 1 2 [L † P ] ii + [L † P ] jj -2[L † P ] ij Q ij = 1 2 r ij Q ij ,
where r ij is the effective resistance between the agents i and j. Since there is an edge between i and j, we have r ij ≤ 1/p. Then we have the following upper bound for the largest eigenvalue:

λ max (Q ij L † Q ij ) ≤ 1 2p for all (i, j) ∈ E, (4.39)
and we get: . Hence, we can also find the optimal step-size for the algorithm, and the corresponding convergence rate.

(i,j)∈E p ij ∥Q ij ∇f (x(k))∥ 2 L † p ≤ 1 2p ||∇f (x(k)) -∇f * || 2 Lp ≤ λ n (L p ) 2p ||∇f (x(k)) -∇f * || 2 , ( 4

Corollary 4.2

The optimal rate of convergence for the RCD algorithm is achieved for h * = 4p (α+β)λn(Lp) which yields

E ∥x(k + 1) -x * ∥ 2 L † p ≤ 1 - 2pλ 2 (L p ) λ n (L p ) 1 κ ∥x(k) -x * ∥ 2 L † p . (4.41)
The upper bound for the step-size derived in Proposition 4.4 is better suited for graphs with a small λ 2 (L p ) (also known as algebraic connectivity), that is non-robust networks that can be easily disconnected [START_REF] Jamakovic | On the robustness of complex networks by using the algebraic connectivity[END_REF]. If we denote by µ the eigenvalues of L, which satisfy λ = p 2 µ, we have that for Proposition 4.

Convergence of RCD in open systems

In the following theorem we present the convergence rate of the RCD algorithm in expectation in an open system. The derivation of this result relies on the separate analysis of the effects of replacements and of RCD updates, which we remind is allowed by Assumption 3.6.

E ∥x(k + 1) -x * (k + 1)∥ 2 L † p -Γ η ≤ A η E ∥x(k) -x * (k)∥ 2 L † p -Γ η , (4.42) with A η := 1 -p U αh(2λ 2 (L p ) -αhλ n (L p )) + (1 -p U ) M n,κ η ; (4.43) Γ η := (1 -p U )M n,κ (η + M n,κ )η p U ηαh(2λ 2 (L p ) -αhλ n (L p )) -(1 -p U )M n,κ , (4.44)
where M n,κ = 1 λ 2 (Lp) Mn,κ , with Mn,κ defined in (4.13).

Proof:

Let us denote C(k) := ∥x(k) -x * (k)∥ 2 L † p
. From Assumption 3.6, there holds

E [C(k + 1)] = p U E [C(k + 1)|U ] + (1 -p U )E [C(k + 1)|R] , (4.45)
where U and R respectively stand for the occurrence of an RCD update and a replacement event. Proposition 4.3 then yields for h ≤ λ 2 (Lp) λn(Lp)

2 α+β E [C(k + 1)|U ] ≤ 1 -2hαλ 2 (L p ) + λ n (L p )α 2 h 2 E [C(k)] . (4.46)
Under a replacement event, there holds x(k + 1) = x(k), and hence Proposition 4.1 yields

E [C(k + 1)|R] = E ∥x(k + 1) -x * (k) + x * (k) -x * (k + 1)∥ 2 L † p ≤ E ∥x(k) -x * (k)∥ L † p + ∥x * (k) -x * (k + 1)∥ L † p 2 = E C(k) + M n,κ 2 . ( 4
E [C(k + 1)] ≤ (1 -p U αh(2λ 2 (L p ) -αhλ n (L p ))) E [C(k)] + (1 -p U ) 2M n,κ E C(k) + M 2 n,κ . (4.48)
Since 2x ≤ η + x 2 η holds for all x ≥ 0 and η > 0, there holds C(k) ≤ η 2 + C(k) 2η for any η > 0. Hence, it follows that

E [C(k + 1)] ≤ 1 -p U αh(2λ 2 (L p ) -αhλ n (L p )) + (1 -p U ) M n,κ η E [C(k)] + (1 -p U )M n,κ (η + M n,κ ).
Solving the linear recurrence then yields the conclusion. ■ Let us now define the ratio

ρ R := 1 -p U p U , (4.49)
which characterizes the expected number of replacements happening in the system between two RCD updates. In particular, when ρ R → 0, then the system converges to a closed system, and when ρ R → ∞, then replacements become so frequent that RCD updates are negligible.

Corollary 4.3 Let η := ρ R Mn,κ αh(2λ 2 (Lp)-αhλn(Lp)) , then for any η > η one has A η < 1. Moreover,

lim sup k→∞ E ∥x(k) -x * (k)∥ 2 L † p ≤ Γ η = η(M n,κ + η) 1 -η/η . ( 4

.50)

Proof: The proof follows from studying

A η < 1 ⇐⇒ p U αh(2λ 2 (L p ) -αhλ n (L p )) + (1 -p U )M n,κ < 0,
and from rewriting Γ η as

Γ η = ρ R M n,κ (M n,κ + η)η αh(2λ 2 (L p ) -λ n (L p )αh)η -ρ R M n,κ = ρ R Mn,κ αh(2λ 2 (Lp)-αhλn(Lp)) (M n,κ + η) 1 - ρ R Mn,κ αh(2λ 2 (Lp)-αhλn(Lp))η
, which directly yields the conclusion. ■ Corollary 4.3 guarantees that one can always choose η such that the contraction rate satisfies A η < 1 for any parametrization of the system (as long as ρ R < ∞, i.e., if updates happen). This guarantees convergence to the corresponding Γ η , defined in (4.50). When η → ∞, then the contraction rate becomes minimal, i.e., A η → 1 -p U αh(2λ 2 (L p ) -λ n (L p )αh), and the asymptotic error becomes extremely large, i.e., Γ η → ∞. Observe that for η > η, then Γ η is convex, and one can determine the value of η that minimizes the asymptotic expected error Γ η , as presented in the following corollary.

A η * = 1 -p U αh(2λ 2 (L p ) -αhλ n (L p )) 1 + Mn,κ η 1 + 1 + Mn,κ η ;
(4.51)

Γ η * = (η * ) 2 = η2 1 + 1 + M n,κ η 2 , (4.52)
where we remind η := ρ R Mn,κ αh(2λ 2 (Lp)-αhλn(Lp)) from Corollary 4.3.

Proof: Observe that Γ η is convex for η > η, and there holds d dη Γ| η=η * = 0 with η * = arg min η>η Γ η . Hence we compute

d dη Γ η = η (η -η) 2 (η 2 -2ηη -M n,κ η) = 0, (4.53)
which is satisfied for

η * 1 = η + η2 + M n,κ η ; η * 2 = η -η2 + M n,κ η.
Since η * 2 < η, it must be rejected, and it follows that

η * = η + η2 + M n,κ η = η 1 + 1 + M n,κ η .
We can then compute

A η * = 1 -p U αh(2λ 2 (L p ) -αhλ n (L p )) + (1 -p U )M n,κ η 1 + 1 + M n,κ /η
, and a few algebraic manipulations yield (4.51). Now observe that Γ η = η 2 if and only if

η 2 -2ηη -M n,κ η = 0,
which is equivalent to (4.53) for η > η, so that the solution is η * . Hence, Γ η * = (η * ) 2 , which yields (4.52). ■ The results of Theorem 4.1 and of Corollary 4.4, and in particular the minimal expected asymptotic error and corresponding convergence rate, allow for some interpretation.

When no replacements happen, i.e., ρ R → 0, then the system behaves as a closed system, and we retrieve the corresponding convergence behavior: the expected asymptotic error Γ η → 0 and the contraction rate A η → 1 -αh(2λ 2 (L p )αhλ n (L p )) for all η > 0, consistently with the convergence rate of the RCD in closed system derived in (4.26). By contrast, as ρ R gets larger, i.e., as replacements become more frequent, then the expected asymptotic error increases, and the contraction rate A η * gets closer to 1 (observe that A η * < 1 remains true as long as ρ R < ∞). In the particular limit case where ρ R → ∞, then the minimal expected asymptotic error becomes Γ η * → 4η → ∞ and A η * → 1. The conservatism of the results essentially follow from Proposition 4.1, whose aim is to bound the additive error injected at one single replacement, whereas tighter bounds might be derived on the sum of those additive errors. This highlights one limitation of our model.

Interestingly, within the allowed range of h, the minimal expected asymptotic error Γ η * decays as the step-size h increases, suggesting that choosing h as large as possible leads to the smallest value of Γ η * . This means that the only limitation on the choice of the step-size comes from the analysis of the algorithm in closed system (i.e., in our case h ≤ λ 2 (Lp) λn(Lp) 2 α+β from Proposition 4.3), and that no particular precaution should be taken regarding the open character of the system.

Remark 4.6

The methodology we used in this section can easily be extended to other algorithms than the RCD algorithm. In particular the results of Theorem 4.1, and hence of Corollaries 4.3 and 4.4, can be adapted to any algorithm with linear convergence in closed system, that is, such that

∥x(k + 1) -x * ∥ 2 L † p ≤ K∥x(k) -x * ∥ 2 L † p , (4.54)
with some positive K < 1. In that case, the same convergence rate as that presented in Theorem 4.1 is obtained with

A η := 1 -p U (1 -K) + (1 -p U ) M n,κ η ; (4.55) Γ η := (1 -p U )M n,κ (η + M n,κ )η p U η(1 -K) -(1 -p U )M n,κ . (4.56)
We can show that convergence can be guaranteed in open system following a similar argument as that used to prove Corollary 4.3 if K < 1. Hence, this analysis can be applied e.g., to the results presented in Proposition 4.4 or in Remark 4.5.

To illustrate the results of Theorem 4.1, we consider systems with piecewise quadratic local cost functions f i : for φ i1 , φ i2 ∈ α 2 , β 2 , the cost function f i is given by

f i (x i ) = φ i1 (x i -ν i ) 2 , if x i < ν i φ i2 (x i -ν i ) 2 , if x i ≥ ν i , (4.57)
where ν i is the minimizer of f i satisfying Assumption 4.1. Such function therefore satisfies Assumption 3.1 as well. Observe that no assumption on the way we choose the local cost function f i of a joining agent at replacements is required in the derivation of our results. Hence, we consider two possible cases for that choice: random, where the parameters φ i1 and φ i2 are uniformly randomly chosen in α 2 , β 2 , and adversarial, where these parameters are arbitrarily chosen to maximize the error ∥x(k)-x * (k)∥ 2 This highlights the impact of those parameters in the tightness of the bound M 2 n,κ used to derive our main results.

In Fig. 4.4, we compare the simulated performance in both replacement cases with the upper bound from Corollary 4.4 for a ring graph with homogeneous agents and a complete graph with heterogeneous agents. By contrast with the previous illustrations, the ring graph setting implies a different, sparse, topology which thus reduces the range of validity for the step-size h due to the small value of λ 2 (L) that does not scale with n. Similarly, the heterogeneous setting impacts λ 2 (L) and consequently reduces the range of h, due to the imbalance in L p . Those moreover affect the behavior of the norm ∥•∥ 2 L † p . Furthermore, the heterogeneous setting influences Mn,κ as well, and hence increases the effect of replacements on the bounds. Nevertheless, even though they differ quantitatively, the results of Fig.

Online optimization analysis

In this section we perform an analysis of the RA problem in OMAS by using performance metrics inspired by online optimization. Due to the potential replacements that occur in the system, the local cost functions of the agents can change in time. Time-varying objective functions can remind an alternative line of work called online optimization [START_REF] Hazan | Introduction to Online Convex Optimization[END_REF][START_REF] Shahrampour | Distributed online optimization in dynamic environments using mirror descent[END_REF], where the approach is to minimize a regret function over a finite period of time T , commonly defined as

Reg s T := T k=1 f k (x(k)) -min x T k=1 f k (x), (4.58)
where the superscript s indicates that the minimum of the sum is computed over a static variable x [START_REF] Shahrampour | Distributed online optimization in dynamic environments using mirror descent[END_REF]. The objective in online optimization is thus to determine the sequence of estimates x t that keeps Reg s T as small as possible over the time period under some assumptions about the possible sequences of time-varying cost functions.

In contrast, the objective of the algorithms in our work is to be at all times as close as possible to the instantaneous minimizer while the possible changes of the cost functions do not follow a specific rule [START_REF] Simonetto | Time-varying convex optimization: Time-structured algorithms and applications[END_REF]. Nevertheless, even if the minimization of the regret is not the objective of optimization in OMAS, a related metric called dynamical regret, which is defined in (4.61), can be used along with two other similar metrics to evaluate the performance of the algorithms. These quantities are defined based on the assumption that the cost of performing an operation by an agent (i.e., the cost function) must be paid at each iteration. In this case, it is natural to consider the accumulated cost in time as a measure of performance of the system.

We present a preliminary analysis by considering a special case of the RA problem in a simplified setting.

Simplified setting

Regarding the RA problem defined in (3.1), where the local cost functions satisfy the Assumption 3.1, we consider the following assumptions: Let S 1,n := x ∈ R n ≥0 : 1 ⊤ x = n be the feasible set, we can now express (3.1) in our setting under the assumptions of this section:

min x∈S 1,n f k (x) = n i=1 f k i (x i ). (4.59)

Performance metrics

Natural indexes for measuring the performance of an algorithm in our setting consist in evaluating its accumulated error over a finite number of iterations with respect to a given strategy. We define the two following strategies of interest in the context of the RA problem:

• Perfect collaboration: at each time instant k the agents know the optimal solution of (4.59) denoted x * (k); x(k) defined in (4.64), the optimal solution x * (k), and the selfish strategy x s (k), in a system subject to replacements of agents (i.e., simultaneous departures and arrivals) on average once every 4 RCD steps.

• Selfish players: the agents do not collaborate to minimize f k , and they operate at their individual desired point so that x s (k) = 1 n at all k.

Hence, for any T , the estimate x(k) obtained with a well-designed algorithm is expected to satisfy

T k=1 f k (x * (k)) ≤ T k=1 f k (x(k)) ≤ T k=1 f k (x s (k)). (4.60)
The evolution of these strategies, compared with that of a given algorithm, is illustrated in Fig. 4.5.

We define the following performance metrics to analyze the value provided by the RCD algorithm x(k) with respect to the strategies above: Dynamical Regret:

Reg T := T k=1 f k (x(k)) -f k (x * (k)) ; (4.

61)

Benefit:

Ben T := T k=1 f k (x s (k)) -f k (x(k)) ; (4.62)
Potential Benefit:

Pot T := T k=1 f k (x s (k)) -f k (x * (k)) . (4.63)
The dynamical regret and benefit respectively measure the accumulated error from using a given algorithm with respect to the optimal solution x * (k) and the accumulated gain from using it instead of the selfish strategy x s (k). The potential benefit is independent of the algorithm; it represents the accumulated advantage of the optimal strategy with respect to the selfish one, and satisfies Pot T = Ben T + Reg T .

Remark 4.7 The regret (4.58), also known as static regret, which is commonly used in online optimization, typically compares x(k) with the overall optimal solution taken over all the iterations, i.e., x * = arg min x∈S 1,n T t=1 f k (x). In that sense, it differs from the dynamical regret in (4.61), which compares x(k) with the time-varying instantaneous optimal solution x * (k) = arg min x∈S 1,n f k (x) at each iteration, such as e.g., in [START_REF] Shahrampour | Distributed online optimization in dynamic environments using mirror descent[END_REF].

RCD algorithm and replacements

For the particular case d = 1, the RCD algorithm (4.4) can be expressed as:

x i (k + 1) = x i (k) - 1 β (f ′ i (x i (k)) -f ′ j (x j (k))) x j (k + 1) = x j (k) - 1 β (f ′ j (x j (k)) -f ′ i (x i (k))). (4.64)
Unlike the assumption used in Section 4.1.1 where the new agent retrieves the value of the replaced agent, in this second part of the chapter we assume that whenever an agent in joins the system, it initializes its estimate as

x in = d in = 1, (4.65)
that corresponds to the selfish strategy of the agent, and whenever an agent out leaves the system, it sends a last message to all its neighbours (i.e., all the other agents in our setting) with its current estimate x out and its demand d out so the agents i ̸ = out update their estimates as

x i (k + 1) = x i (k) + x out (k) -x i (k) n = 1 - 1 n x i (k) + 1 n x out (k). (4.66)
We show in the following proposition that RCD iterations, arrivals and departures as they are defined in (4.64) to (4.66) guarantee that as long as the initial estimate x 0 is feasible, then all the estimates remain feasible. Proof: We first consider arrivals: the nonnegativity of x i and preservation of the constraint is a direct consequence of (4.65). In the case of departures, the nonnegativity of x i is a direct consequence of (4.66). Moreover, if the constraint is satisfied at iteration k with n k = n, then under the departure of the agent labelled out we have i̸ =out

x i (k + 1) = n -x out (k) + x out (k) n (n -1) - n -x out (k) n = n -1.
We finally consider iterations of the RCD algorithm. From Assumptions 3.1, 4.2 and 4.3, it follows that for any x i ≥ 0

x i (k)f ′ i (x i (k)) ≥ α |x i (k)| 2 ≥ 0, so that f ′ i (x i (k)) ≥ 0. Moreover, since f i is β-smooth, one has f ′ i (x i (k)) ≤ βx i (k)
, and therefore at each update of the RCD algorithm between agents i and j there holds

x i (k + 1) = x i (k) - 1 2β f ′ i (x i (k)) -f ′ j (x j (k)) ≥ x i (k) - 1 2β (βx i (k)) = x i (k) 2 ,
establishing the nonnegativity of x i (k). A similar analysis can be used for x j (k).

Due to the symmetry of the update rule, the constraint is always preserved and we conclude the proof. ■ Our goal is to analyze the performance of the RCD algorithm (4.64) with the arrival and departure rules (4.65) and (4.66) in the setting described in Section 4.5.1 using the metrics defined in Section 4.5.2 in expectation.

Upper bounds on the performance metrics

We now derive upper bounds on the evolution of the Potential Benefit and the Dynamical Regret respectively defined in (4.63) and (4.61) in expectation. Whereas the former is only related to the problem itself, the latter actually depends on the algorithm we consider.

We first provide the following lemmas, where Lemma 4.3 directly follows from the equivalence of the norms.

Lemma 4.3 Let x ∈ S 1,n , then n ≤ ||x|| 2 ≤ n 2 . Lemma 4.4 Let f (x) = n i=1 f i (x i )
, where all f i satisfy Assumptions 3.1, 4.2 and 4.3, then for any x ∈ S 1,n there holds

α 2 n ≤ f (x) ≤ β 2 n 2 . ( 4

.67)

Proof: From Assumptions 3.1, 4.2 and 4.3, f (0) = f ′ (0) = 0. Hence, since f is β-smooth and using Lemma 4.3, there holds f (x) ≤ β 2 ||x|| 2 ≤ β 2 n 2 which establishes the upper bound. Similarly, since f is α-strongly convex and by using Lemma 4.3, it follows that f (x) ≥ α 2 ||x|| 2 ≥ α 2 n, which establishes the lower bound, and concludes the proof. ■ Lemma 4.4 provides a global upper bound on the difference between any two solutions x, y ∈ S 1,n :

f k (x(k)) -f k (y(k)) ≤ n 2 (nβ -α). (4.68)
This can be used to derive upper bounds on any of the metrics defined in Section 4.5.2, e.g., Ben T ≤ n 2 (nβ -α)T .

Potential Benefit

We first obtain in the following theorem an upper bound on the expected value of the potential benefit, which we remind quantifies the accumulated advantage of using the optimal strategy rather than not collaborating at all. Theorem 4.2 In the setting of Section 4.5.1, there holds

Pot T ≤ n 2 α (κ -1) T, (4.69)
and in particular lim

T →∞ Pot T T ≤ n 2 α(κ -1). (4.70)
Proof: Remember that x s (k) = 1 n by definition, and that f k (0 n ) = 0 and ∇f k (0 n ) = 0 n from Assumption 4.3. Hence, there holds from the β-smoothness

of f k f k (x s (k)) ≤ ∇f k (0 n ) ⊤ (x s (k)) + β 2 ||x s (k)|| 2 = β 2 ||1 n || 2 = β 2 n.
Similarly, since f k is α-strongly convex, we get

f k (x * (k)) ≥ ∇f k (0 n ) ⊤ (x * (k)) + α 2 ||x * (k)|| 2 = α 2 ||x * (k)|| 2 ≥ α 2 n,
where the last inequality follows from Lemma 4.3. Hence

f k (x s (k)) -f k (x * (k)) ≤ β 2 n - α 2 n = n 2 (β -α)
holds, and injecting it into (4.63) yields (4.69). The last result then follows from dividing (4.69) by T . ■ Notice that when the constraint is always preserved, the dynamical regret is nonnegative and the bounds (4.69) and (4.70) also hold for the benefit since Ben T = Pot T -Reg T ≤ Pot T .

Dynamical Regret

We now obtain an upper bound on the expected dynamical regret defined in (4.61), where we remind x(k) is obtained with the RCD algorithm defined in Section 4.5.3. For that purpose, we first introduce the following intermediate quantities:

C(k) := f k (x(k)) -f k (x * (k)); (4.71) ∆f (k) := f k+1 (x(k + 1)) -f k (x(k)); (4.72) ∆f * (k) := f k+1 (x * (k + 1)) -f k (x * (k)). (4.73)
Thus, C(k) corresponds to the instantaneous loss of the RCD algorithm with respect to the optimal solution at iteration k, and ∆f (k) and ∆f * (k) respectively stand for the instantaneous variation at one iteration of the total estimated cost and optimal cost, such that

C(k + 1) = C(k) + ∆f (k) -∆f * (k).
In the following proposition, we study the effect or replacements on ∆f (k) in order to later characterize C(k) in expectation, and consequently the expected dynamical regret. Proposition 4.6 In the setting of Section 4.5.1 the replacement of an agent, denoted R, results in

E [∆f (k) | R] ≤ 5 2 β - 3 2 α. (4.74)
Proof: We analyze the effects of arrivals and departures separately. Let g denote the local cost function of the joining agent at an arrival, then

f k+1 (x(k + 1)) = f k (x(k)) + g(1) and ∆f (k) = f k (x(k)) + g(1) -f k (x(k)) = g(1) ≤ β 2 , (4.75)
where the last inequality follows from Assumptions 3.1 and 4.3, and in particular the β-smoothness of g. Consider now a departure, and let ℓ denote the label of the leaving agent, such that

f k+1 (x(k + 1)) = i̸ =ℓ f k i (x i (k + 1)), with x i (k + 1) = x i (k) + x i (k)+x ℓ (k)
n from (4.66). From the definition of departures, ℓ is uniformly selected among the n agents in the system and by taking the expected value ∆f (k) over the leaving agent, one gets the following, where we omit the reference to time to lighten the notation:

E [∆f (k)] = n ℓ=1 1 n i̸ =ℓ f i x i + x ℓ -x i n -f (x) = 1 n n ℓ=1 f n -1 n x + x ℓ n 1 n -f ℓ (x ℓ ) -f (x) = 1 n n ℓ=1 f n -1 n x + x ℓ n 1 n - n + 1 n f (x).
Since f is β-smooth from Assumption 3.1, one has

f n -1 n x + x ℓ n 1 n ≤ f (x) + 1 n ⟨∇f (x), x ℓ 1 n -x⟩ + β 2n 2 ||x ℓ 1 n -x|| 2 ,
and it follows that

E [∆f (k)] ≤ 1 n 2 n ℓ=1 ⟨∇f (x), x ℓ 1 n -x⟩ + β 2n 3 n ℓ=1 ||x ℓ 1 n -x|| 2 - 1 n f (x). (4.76)
From Assumption 3.1, in particular since f is β-smooth and α-strongly convex, it satisfies α ||x|| 2 ≤ ⟨∇f (x), x⟩ ≤ β ||x|| 2 for any x. Hence, reminding that n ℓ=1 x ℓ = n, the first sum of (4.76) can be upper bounded by

n ℓ=1 ⟨∇f (x), x ℓ 1 n -x⟩ = n⟨∇f (x), 1 n -x⟩ ≤ n(βn -α ||x|| 2 ).
The second sum of (4.76) can be expressed as:

n ℓ=1 ||x ℓ 1 n -x|| 2 = n ℓ=1 n i=1 x 2 ℓ -2x ℓ x i + x 2 i = n ℓ=1 nx 2 ℓ -2nx ℓ + ||x|| 2 = 2n ||x|| 2 -n .
Then (4.76) is upper bounded by

E [∆f (k)] ≤ 1 n βn -α ||x|| 2 + β n 2 ||x|| 2 -n - 1 n f (x).
Lemmas 4.3 and 4.4 yield ||x|| 2 ≤ n 2 and f (x) ≥ α 2 n, so that

E [∆f (k)] ≤ 2β - 3 2 α. (4.77)
The conclusion follows from adding (4.75) and (4.77). ■ We can now use Proposition 4.6 to study the evolution of the expected dynamical regret in the following theorem.

Theorem 4.3 In the setting of Section 4.5.1, there holds

E [Reg T ] ≤ C 0 T k=1 η k + T -1 k=0 η k (M f + (1 -p U )(T -k)θ) , (4.78)
where η = 1 -p U κ(n-1) (with p U the probability that a given event is an update from Assumption 3.6),

M f = n 2 (βn -α) and θ = 5 2 β -3 2 α. Proof: Let γ = 1 -1 κ(n-1)
denote the contraction rate of the RCD algorithm as defined in (4.64) [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF]. Remember that there holds C(k + 1) = C(k) + ∆f (k) + ∆f * (k) for all times k. Hence, at any time-step k one has

E [C(k + 1)] = E [C(k) + ∆f (k) -∆f * (k)] . (4.79)
From Assumption 3.6 the event at iteration k is an update, denoted U k , with probability p U , or a replacement, denoted R k , with probability 1 -p U .

In the case of an update, we have x * (k + 1) = x * (k), so that ∆f * (k) = 0. Hence, we have ∆f (k) = C(k + 1) -C(k), and since E [C(k + 1)|C(k), U k ] ≤ γC(k) with the RCD algorithm from [START_REF] Necoara | Random coordinate descent algorithms for multi-agent convex optimization over networks[END_REF]:

E [∆f (k)|, U k ] = E [C(k + 1) -C(k)|U k ] ≤ (γ -1)E [C(k)] . (4.80)
In the replacement case, there holds

E [∆f (k)|R k ] ≤ θ = 5 2 β - 3 2 α, (4.81)
where θ comes from Proposition 4.6. Injecting (4.80) and (4.81) into (4.79) then yields

E [C(k + 1)] ≤ E [C(k)]+p U (γ -1)E [C(k)]+(1 -p U )θ -E [∆f * (k)] = ηE [C(k)] + (1 -p U )θ -E [∆f * (k)] , (4.82)
where η = 1 + p U (γ -1). Expression (4.82) actually describes the evolution of a discrete-time dynamical system of the type v(k + 1) ≤ Av(k) + Bu(k). Standard results on that framework yield v(k) ≤ A k v(0)+ k-1 j=0 A k-j-1 Bu(j), and we obtain

E [C(k)] ≤ η k C 0 + k-1 j=0 η k-j-1 ((1 -p U )θ -E [∆f * (k)]) . (4.83)
Injecting this last result into (4.61) then yields

E [Reg T ] = T k=1 E [C(k)] ≤ C 0 T k=1 η k + T k=1 k-1 j=0 η k-j-1 ((1 -p U )θ -E [∆f * (k)]) .
After some term re-organization, it becomes

E [Reg T ] ≤ C 0 T k=1 η k + (1 -p U ) T -1 k=0 (T -k)η k θ - T -1 k=0 η k T -k j=0 E ∆f * j .
Finally, using Lemma 4.4, one concludes that

- T -k j=0 E ∆f * j = -E T -k j=0 ∆f * j = -E f T -k+1 (x * (T -k + 1)) -f 1 (x * (1)) ≤ n 2 (βn -α),
and M f = n 2 (βn -α) yields the conclusion. ■ We now analyze the asymptotic behavior of the averaged regret in the following corollary.

Corollary 4.5 Let ρ R := 1-p U p U . In the same setting as that of Theorem 4.3, there holds

lim T →∞ E [Reg T ] T ≤ ρ R (n -1)β 5κ -3 2 . ( 4.84)
Proof: Starting from (4.78), we have

E [Reg T ] T ≤ C 0 T k=1 η k T + T -1 k=0 M f η k T + (1 -p U ) 1 - k T η k θ .
Remember that η = 1-p U κ(n-1) ≤ 1, so that T k=1 η k < T , and lim T →∞

T k=1 η k T = 0. Hence lim T →∞ E [Reg T ] T ≤ lim T →∞ (1 -p U ) T -1 k=0 1 - k T η k θ.
Moreover, for η < 1, one shows ∞ j=0 η j = 1 1-η and ∞ j=0 jη j = η (1-η) 2 . The latter implies that lim T →∞ 1 T T -1 k=0 kη k = 0, and therefore there holds

lim T →∞ E [Reg T ] T ≤ (1 -p U ) 1 1 -η θ = 1 -p U p U κ(n -1)θ,
and the conclusion follows from the definitions of θ in Theorem 4.3, and from ρ R := 1-p U p U . ■ The upper bounds for the potential benefit (4.69) and the dynamical regret (4.78) linearly scale with T . This behavior is rather natural for the former, which does not depend on the algorithm. For the latter, it is most likely unavoidable due to the introduction at each replacement of perturbations of non-decaying magnitude and without any regularity [START_REF] Li | A survey of decentralized online learning[END_REF]. Interestingly, this behavior contrasts with standard results in online optimization, where a sublinear growth in T is desired to cancel the asymptotic averaged regret [85, Ch. 1.1]. However, these results usually apply on the static regret (4.58), where x(k) is compared with an overall timeindependent strategy x * computed over all T iterations, in opposition with x * (k) which is optimal for each iteration. In the case of the dynamical regret, it has been showed that without any regularity in the variations of the cost functions, a sublinear behavior cannot be achieved [START_REF] Li | A survey of decentralized online learning[END_REF].

Moreover, the corresponding asymptotic upper bounds linearly grow with n and α for (4.70), and with n -1 and β for (4.84), consistently with their expected behavior. In particular, the scaling of (4.84) with n-1 follows from the convergence rate of the RCD algorithm γ = 1 -1 κ(n-1) . Interestingly, (4.84) is proportional to ρ R (n -1) = (1 -p U ) n-1 p U , and the bound can thus be seen as the ratio between the probability for a given agent to be involved in a RCD update p U n-1 (involved in γ), and the impact of replacements at the system level 1 -p U , independently of n. This is consistent with the bound on the impact of replacements in (4.74) which is independent of n (by contrast, alternative situations such as e.g., if all agents were to be reset at each replacement are expected to generate an impact growing with n). Hence, for small values of ρ R (i.e., rare replacements), the bound guarantees that the asymptotic dynamical regret remains reasonably bounded, and decays to zero when ρ R → 0, i.e., for closed systems.

Finally, observe that (4.84) is proportional to 5κ-3 2 , consistently with the fact that a larger interval for the possible curvature of the cost functions should generate a larger potential error at replacements. This factor is a potential source of conservatism, e.g., with respect to (4.70) where the scaling is in 1 2 (κ -1). More generally, it is not clear yet whether other algorithms than the RCD might provide tighter bounds.

Remark 4.8 The proofs of Theorem 4.3 and Corollary 4.5 can directly be adapted to any contraction rate γ < 1, and are thus easily generalized to any other algorithm that guarantees linear convergence; in particular lim T →∞

E[Reg T ] T ≤ ρ R θ 1-γ .

The case of quadratic functions

The bound on the expected dynamical regret can be refined for the particular case where all local functions are quadratic, i.e., satisfy the following additional assumption.

Assumption 4.6 (Quadratic functions) The local cost function of any agent i at time t is of the form

f i (x i ) = ϕ i x 2 i , ϕ i ∈ α 2 , β 2 . (4.85)
The parameter ϕ i is randomly chosen according to a distribution with a finite support determined by the interval α 2 , β 2 . Observe that functions satisfying Assumption 4.6 necessarily satisfy Assumption 3.1 as well.

Under Assumption 4.6, we can obtain a tighter bound than that of Proposition 4.6, presented in the following proposition. Proposition 4.7 In the setting of Section 4.5.1, and under Assumption 4.6, the replacement of an agent R results in

E [∆f (k) | R] ≤ 3n 2 -3n + 1 2n 2 (β -α) . (4.86)
Proof: The arrival case is treated the same was as in the proof of Proposition 4.6, resulting in ∆f (k) ≤ β 2 . The departure case follows the same first steps with f i (x) = ϕ i x 2 , and

E [∆f (k)] = n ℓ=1 1 n i̸ =ℓ ϕ i x i + x ℓ -x i n 2 -f (x) = 1 n n ℓ=1 i̸ =ℓ ϕ i x i + x ℓ -x i n 2 -x 2 i -ϕ ℓ x 2 ℓ = 1 n 2 n ℓ=1 i̸ =ℓ ϕ i 1 -2n n x 2 i + 2x i x ℓ n -1 n + x 2 ℓ n - f (x) n .
Using the fact that n ℓ=1 i̸ =ℓ ϕ i x 2 i = (n -1)f (x) and that ϕ i ≤ β 2 for all i, we obtain

E [∆f (k)] ≤ 1 n 2 1 -2n n (n -1) - 1 n f (x) + β 2n 2 n ℓ=1 i̸ =ℓ 2x i x ℓ n -1 n + x 2 ℓ n .
Observe that n ℓ=1 i̸ =l x i x ℓ = n 2 -||x|| 2 and n ℓ=1 i̸ =ℓ x 2 ℓ = (n -1) ||x|| 2 , so that a few algebraic manipulations yield

E [∆f (k)] ≤ - 3n 2 -3n + 1 n 3 f (x) + β n -1 2n 3 2n 2 -||x|| 2 .
Using ||x|| 2 ≥ n (from Lemma 4.3) and f (x) ≥ α 2 n (from Lemma 4.4) then yields

E [∆f (k)] ≤ - 3n 2 -3n + 1 2n 3 α + β 2n 2 -3n + 1 2n 3 ,
and combining with the arrival case concludes the proof. ■ The result above allows us stating the following theorem, which improves Theorem 4.3 and Corollary 4.5 respectively for the case of quadratic functions. Theorem 4.4 In the setting of Section 4.5.1, and under Assumption 4.6, there holds

E [Reg T ] ≤ C 0 T k=1 η k + T -1 k=0 η k (M f + (1 -p U )(T -k)θ) , (4.87)
where η = 1 -p U κ(n-1) (with p U the probability that a given event is an update from Assumption 3.6),

M f = n 2 (βn -α) and θ = (β -α) 3n 2 -3n+1 2n 2
. In particular,

lim T →∞ E [Reg T ] T ≤ ρ R (n -1) 3n 2 -3n + 1 2n 2 β (κ -1) . (4.88)
Proof: The proof follows the exact same steps as those of Theorem 4.3 and Corollary 4.5 where Proposition 4.7 is used instead of Proposition 4.6.

■ The upper bound for the quadratic case is qualitatively better than that of the general case; it was derived based on an additional information of the cost function, thus resulting in tighter bounds. In this case, the dependence of (4.88) is in 3 2 (κ -1), which is consistent with the result derived for the potential benefit (4.70). In particular, for n becoming large, (4.88) and (4.84) become equivalent up to a constant β. Moreover, (4.88) becomes 0 when κ = 1, consistently with the expected behavior of the RCD algorithm for quadratic functions since all the cost functions would then be the same.

Numerical results

To illustrate the results of Theorems 4.2 to 4.3, we consider a system of 5 agents with κ = 10 and ρ R = 0.0125, the latter implies that on average there is one replacement every 80 events. We consider two possibilities: random replacements (RR) where the local function is randomly uniformly chosen among the set of piecewise quadratic functions satisfying Assumptions 3.1, 4.2 and 4.3, and adversarial replacements (AR), where these functions are quadratic functions ϕ i x 2 , with

ϕ i ∈ α 2 , β 2 .
The AR setting is expected to be less favorable than the RR setting, since replacements might result in the largest change of local functions. Notice that the bounds (4.69) and (4.78) are independent of the distribution from which the local cost functions are assigned to the agents when they join the system, so that they hold for any such assignment rule. Fig. 4.6 compares the results of Theorem 4.2 and Corollary 4.5 with simulations for both random and adversarial replacements in the setting described above. Even though the theoretical bounds are conservative, they capture well the qualitative behavior of these metrics. In particular, consistently with Pot T and E [Reg T ] that grow linearly with T , the bounds in the figure do not converge to zero, and a remaining asymptotic error is observed. Our bounds are tighter for the adversarial replacement case than the random replacement case, this suggests that our bounds might be tight for some particular choice of the joining functions at replacements, especially that on the potential benefit. the same setting as described previously, with random replacements (RR). The figure shows that bound (4.88) is tighter for quadratic functions, following the fact that we have access to more information regarding the local cost functions, thus improving the estimation of the effect of replacements on the expected regret.

Conclusion

In this chapter, we studied the Random Coordinate Descent algorithm in a resource allocation problem by analyzing the distance to the minimizer and by using metrics inspired from online optimization. Regarding the behavior of the distance to the minimizer, we proved linear convergence of the RCD algorithm in an appropriate norm for the closed system. We analyzed the algorithm for a general graph topology and possible heterogeneous agents in OMAS when agents can be replaced during the iterations. Under replacement events, we proved that for an appropriate step-size, the algorithm cannot converge to the instantaneous minimizer due to the perturbations generated by the replacements but is stable. We derived an upper bound for the error in expectation which depends on the variation of the minimizer due to replacements and the frequency of these events. Inspired by the similarities with the online optimization, we considered a simple preliminary setting where the budget is homogeneous and the graph is complete, and show that it is not possible to achieve convergence to the optimal solution with the RCD algorithm in expectation in OMAS, but that the error is expected to remain reasonable. We have derived upper bounds on the evolution of the dynamical regret and the potential benefit in expectation and showed that due to the random choice of the new local cost function during replacements, an error is expected to be accumulated with time and cannot be compensated. A natural continuation of this work would be to extend the results to more general graph topologies, improve the current bounds for the three metrics and derive lower bounds that allow us to support the analysis regarding the linear behavior of the metrics.

This chapter opens several questions for future research. Scenarios where the budget in the constraint is time-varying looks more appropriate for OMAS where the occurrence of a replacement may imply a modification of the constraint in the optimization problem. In this case, it would be interesting to analyze possible violations of the constraint as in Chapter 3. In addition, the case of agents replaced with new values of the corresponding weights a i is also linked to this future line of research, where the contributions of new agents to the budget can be different.

A possible varying size of the system is an interesting direction for future work, where agents could join and leave the network independently of the current state of the system. In this scenario, it would be important to specify the potential changes in the constraint due to arrivals and departures since the budget of the system would be defined by the sum of the states of a different number of agents.

In this chapter, the constraint in the RA problem was simply defined as the weighted sum of the states of the agents. However, it would be interesting to analyze scenarios where the constraint has the form Ax = b, where the matrix A is used to model more complex relationships between the state variables of the agents. Finally, block updates at each iteration [START_REF] Richtárik | Iteration complexity of randomized blockcoordinate descent methods for minimizing a composite function[END_REF], looks also as an interesting direction for future work, where at each iteration the optimization is performed along more than one edge, such that the convergence of the algorithm is faster.

Chapter 5 Contagion in OMAS

In this chapter, we consider a continuous time formulation of a SIS epidemic over a network of individuals and perform an analysis of the system when replacements take place. Section 5.1 introduces the general characteristics of epidemic models used to model the spread of diseases. In Section 5.2 we present the nonlinear model of a SIS epidemic over a network and the aggregate function used for the analysis. We also present the linearized version of the model, which will be used for the derivation of an upper bound for the variance of the aggregate function. In Section 5.3 we analyze the SIS epidemic when replacements of agents take place during the evolution of the disease. We consider that the time instants at which replacements occur are determined by a homogeneous Poisson process. The expectation and variance of the aggregate function for the system are studied by deriving an upper bound for their asymptotic values when the network of connections is given. Additionally, we provide an analysis of the expectation of the aggregate function when the graph is sampled from a piecewise Lipschitz graphon. We give some insights into the analysis of more complex scenarios where the connections of the network can change during the replacements according to a piecewise Lipschitz graphon, and we formulate a conjecture on a potential upper bound for the asymptotic value of the aggregate function in expectation. Finally, in Section 5.4 we present the conclusions and future work of this chapter.

Epidemic models

Epidemic models are dynamical systems used to model the spread of a disease in a community of individuals. The importance of epidemic models has grown in the last years, specially due to the COVID-19 outbreak that impacted the entire world. Different types of mathematical models were used to describe the spread of epidemics. For instance, mean field approximations, non-markovian epidemics, per-colation process, partial differential equations, among others, [START_REF] Kiss | Mathematics of Epidemics on Networks[END_REF], which model specific characteristics and properties of the spread of the disease. Depending on the nature of the disease, the models can include different compartments (division of the population); for instance, M for immune class, S for susceptible class, E for exposed class, I for infective class and R for recovered class [START_REF] Hethcote | The mathematics of infectious diseases[END_REF].

One of the most important models is the SIS (Susceptible-Infected-Susceptible) epidemic model, which describes diseases where the agents, after being recovered, can acquire the disease again. The dynamics of this epidemic is determined by the following equation [START_REF] Kiss | Mathematics of Epidemics on Networks[END_REF]:

ẋ(t) = βx(t)(1 -x(t)) -δx(t),
where x(t) is the fraction of a fixed population infected at time t.

In real life, the evolution of the epidemic is not the same in all the infected agents since each individual may be exposed to different external conditions that can have a significant impact in the infection and recovery processes. One of the most important factors that contribute to the propagation of a disease is the network of contacts of an individual, which allow to model the influence of social interactions in the evolution of the disease. SIS models have been extended to the spread of diseases over networks where nodes can have two different interpretations. If the node i is considered as a small population, x i (t) represents the fraction of infected people at time t. If the node i is considered as an individual, x i (t) represents its probability of infection at time t [START_REF] Pastor-Satorras | Epidemic processes in complex networks[END_REF].

Most of the stability results derived using network models depend on the spectral radius of the adjacency matrix associated with the network, assuming that the composition of the system is fixed [START_REF] Nowzari | Analysis and control of epidemics: A survey of spreading processes on complex networks[END_REF]. Nevertheless, in more realistic environments, individuals do not stay in a fixed place for a long time and are constantly moving between different locations. Different works begin to consider the mobility as an important aspect that must be taken into account to obtain more accurate results in the study of the spread of diseases [START_REF] Tizzoni | On the use of human mobility proxies for modeling epidemics[END_REF][START_REF] Hazarie | Interplay between population density and mobility in determining the spread of epidemics in cities[END_REF]. In this case, open multi-agent systems appear as an important tool for the analysis of more realistic epidemics where the set of agents is time-varying. Even, if the total number of individuals in a specific place could remain almost invariant, it is unrealistic to assume that exactly the same agents stay in a fixed location during all the evolution of the epidemic. If we consider a similar rate of departures and arrivals, which preserve the number of individuals, epidemics under these conditions can be modeled as a dynamical system subject to replacements of individuals.

SIS networked model

We consider a system composed by n agents interacting through a connected undirected graph G = (V, E). The SIS epidemic model over a network is given by [START_REF] Nowzari | Analysis and control of epidemics: A survey of spreading processes on complex networks[END_REF]:

ẋi (t) = -δx i (t) + n j=1 a ij βx j (t) (1 -x i (t)) ,
where x i (t) ∈ [0, 1] is the probability of infection of an individual at time t, δ is the recovery rate, β is the infection rate, a ij = 1 if there is an interaction between agents i and j, and a ij = 0 otherwise.

The model of all the network can be expressed as:

ẋ(t) = (βA -δI) x(t) -βX(t)Ax(t), (5.1)
where A = [a ij ] is the adjacency matrix and X(t) = diag[x 1 (t), . . . , x n (t)].

Proposition 5.1 (Stability [START_REF] Nowzari | Analysis and control of epidemics: A survey of spreading processes on complex networks[END_REF]) For the SIS epidemic (5.1), the equilibrium x = 0 is globally asymptotically stable if and only if

λ 1 (A) β δ < 1, (5.2)
where λ 1 (A) is the largest eigenvalue of A.

Since in OMAS the size of the system may change in time due to decoupled arrivals and departures, it is desirable to perform the analysis with the use of an aggregate function V (x(t)) that condensates the associated process in a scalar value, independently of the size of the system. In the case of SIS epidemics, a natural choice is the Lyapunov function

V : [0, 1] n → R ≥0 : V (x(t)) = 1 n ||x(t)|| 2 , (5.3)
which satisfies

lim t→∞ V (x(t)) = 0, (5.4)
when the epidemic is stable. This particular definition of V (x(t)) provides uniform bounds for the aggregate function, independently of the size of the system: 0 ≤ V (x(t)) ≤ 1 for all t.

(5.5) In a closed system, where the evolution of the states of the agents is fully determined by (5.1), the dynamics of the aggregate function V (x(t)) can be upper bounded by [START_REF] Khanafer | Stability properties of infection diffusion dynamics over directed networks[END_REF]:

V (x(t)) = 2 n x T (t) ẋ(t) = 2 n x T (t) (βA -δI) x(t) -βx T (t)X(t)Ax(t) (5.6)
≤ 2 n x T (t) (βA -δI) x(t) ≤ 2(βλ 1 (A) -δ)V (x(t)), (5.7)
where we used the fact that x T (t)X(t)Ax(t) ≥ 0 since x i ∈ [0, 1]. Clearly, the right-hand size of (5.7) is negative if βλ 1 (A) -δ < 0, which corresponds to the stability condition (5.2). Finally, by applying the Comparison Lemma, we have that V (x(t)) satisfies:

V (x(t)) ≤ V 0 + t 0 2(βλ 1 (A) -δ)V (x(τ ))dτ, (5.8)
where V 0 denotes the value of V (x(t)) at the time t = 0.

To illustrate the behavior of the upper bound (5.8), we consider a graph with n = 11 agents represented in Fig. 5.1, whose adjacency matrix is given by

A =                  
0 1 1 1 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 1 0 1 0 0 0 0 0 1 0 1 0 1 1 0 0 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 1 1 0 0 0 1 0 0 1 0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 1 0 0 1 0 0 0 0 0 0 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 1 0 0

                 
(5.9)

Figure 5.2 presents the evolution of V (x(t)) for the graph determined by (5.9) with β = 0.5 and δ = 5.3, where the dashed red line corresponds to the right-hand side of (5.8).

Linearized model

If the probabilities of infection of the agents x i (t) remain close to the equilibrium point x = 0, it is possible to analyze the behavior of the epidemic by using the linearization of (5.1) near the origin [START_REF] Khanafer | An optimal control problem over infected networks[END_REF]:

ẋ(t) = (βA -δI) x(t).
(5.10)

This linearization is exponentially stable [START_REF] Mei | On the dynamics of deterministic epidemic propagation over networks[END_REF] under condition (5.2). Another motivation to consider the linear model in the analysis is that the frequency of occurrence of arrivals, departures and replacements is low compared to the natural evolution of the epidemic. For this reason, the system will remain close to the equilibrium point most of the time and it will only be perturbed at specific time instants. However, the jumps must be small since the analysis is local.

Similar to the nonlinear model, the aggregate function V (x(t)) for (5.10) also satisfies (5.7) and (5.8). However, unlike the nonlinear case, the dynamics of V (x(t)) can also be lower bounded by:

V (x(t)) = 2 n x T (t) ẋ(t) = 2 n x T (t) (βA -δI) x(t) ≥ 2(βλ n (A) -δ)V (x(t)), (5.11)
where λ n (A) is the smallest eigenvalue of A. By applying the Comparison Lemma, the aggregate function V (x(t)) also satisfies:

V (x(t)) ≥ V 0 + t 0 2(βλ n (A) -δ)V (x(τ ))dτ.
The lower bound (5.11) will be used for the derivation of a tighter upper bound for the variance of V (x(t)) when the linearized model (5.10) is considered for the analysis.

Replacements

What is a replacement?

In this subsection, we analyze the case when in addition to the continuous evolution of the states of the individuals according to (5.1), the system also experiments replacements of individuals consisting of a departure of an agent, which stops spreading the infection in the network, followed immediately by the addition of a new infected agent.

In Chapters 3 and 4 we define the replacements of agents as the change of the cost function associated with the agents. However, in the case of a SIS epidemic model, the replacement of an agent is related to the change of the probability of infection associated with the agent 1 . Definition 5.1 We say that an agent j ∈ {1, . . . , n} is replaced at time t if:

x + = [x - 1 , • • • , x - j-1 , x r , x - j+1 , • • • , x - n ] T , x r ̸ = x - j ,
where x -= x(t -) is the state of the system before the replacement, x + = x(t + ) is the state of the system after the replacement and x r is the value of the new agent.

Assumption 5.1 The set of edges E of the network remains invariant for all time t.

In the following proposition we analyze the variation of the aggregate function V (x(t)) under a replacement event where the replaced agent is chosen uniformly and the value of the new agent is determined by a continuous distribution with support in the interval [0, 1], with mean m and variance σ 2 . Proposition 5.2 During the replacement of an agent, the aggregate function V (x(t)) defined in (5.3) satisfies:

E V (x(t + )) -E V (x(t -)) = 1 n σ 2 + m 2 -E V (x(t -)) . (5.12)
Proof: Let us assume that an agent j is replaced. We begin by computing the conditional expectation of V (x(t + )) given x(t -) and the value of the replaced agent Θ. Then, it holds

E V (x(t + ))|x(t -), Θ = 1 n n i=1 1 n j̸ =i x 2 j (t -) + 1 n Θ 2 = 1 n 2 n i=1 x(t -) 2 -x 2 i + Θ 2 = 1 n 2 n i=1 x(t -) 2 - 1 n 2 n i=1 x 2 i (t -) + 1 n 2 n i=1 Θ 2 = 1 n x(t -) 2 - 1 n 2 x(t -) 2 + Θ 2 n = V (x(t -)) - 1 n V (x(t -)) + Θ 2 n .
Then, we compute the total expectation and we obtain

E E V (x(t + ))|x(t -), Θ = E V (x(t -)) - 1 n E V (x(t -)) + σ 2 + m 2 n ,
which completes the proof. ■

E [V (x(t))] = V 0 + t 0 E V (x(τ + )) -V (x(τ -)) λ r dτ,
which can be expressed as the ODE:

dE [V (x(t))] = λ r E V (x(t + )) -V (x(t -)) dt.
(5.13) Proposition 5.3 Under Assumption 5.2 and inter-event dynamics ẋ(t) = 0, the aggregate function V (x(t)) defined in (5.3) satisfies:

E [V (x(t))] = V 0 -σ 2 -m 2 e -λr n t + (σ 2 + m 2 ).
(5.14)

Proof: The result is obtained by applying (5.12) in (5.13) and solving the corresponding ODE:

d dt E [V (x(t))] = - λ r n E [V (x(t))] + λ r n (σ 2 + m 2 ).
■ Proposition 5.3 shows that in a system subject only to replacements, the expected value of the aggregate function is bounded and its asymptotic value is: 3) in a system subject only to replacements with λ r = 30. In the right plot, the solid blue line corresponds to the computation of V (x(t)) considering 10000 realizations of the process while the dashed red line is (5.14).

lim t→∞ E [V (x(t))] = σ 2 + m 2 , ( 5
guarantee a fast convergence. The asymptotic value of V (x(t)) is independent of the rate λ r . Figure 5.4(a) presents the evolution of V (x(t)) for one realization of a replacement process following a Poisson process with intensity λ r = 30 and where the new value of the agent is taken from a uniform distribution U[0, 1].

SIS epidemic with replacements on a given graph

Now, we consider the behavior of the aggregate function in continuous time subject to replacements at time instants determined by a Poisson process. Since the solution to (5.6) is unique and the jumps are determined by a homogeneous Poisson process, the aggregate function V (x(t)) is given by [START_REF] Brockett | Stochastic control[END_REF]:

V (x(t)) = V 0 + t 0 2 n x T (τ ) (βA -δI) x(τ ) -βx T (τ )X(τ )Ax(τ ) dτ + t 0 ∆V k dN τ , (5
V (x(t)) = V 0 + T 1 0 ω(x(τ ))dτ + T 2 T 1 ω(x(τ ))dτ + • • • + Nt k=1
∆V k , which shows that V (x(t)) jumps at the time instants T k determined by the Poisson process, corresponding to the sum, and between these time intervals, V (x(t)) evolves according to (5.6), corresponding to the integrals.

Even if the solution of (5.17) is unique, every realization of the stochastic process can be completely different and the analysis of the asymptotic behavior is not appropriate because the limit lim t→∞ V (x(t)) does not exist due to the continuous replacements. For this reason, we are interested in the statistical properties of (5.17), corresponding to E [V (x(t))].

Theorem 5.1 (Expectation) Consider a SIS epidemic satisfying the stability condition (5.2). Under Assumptions 5.1 and 5.2, the aggregate function V (x(t)) defined in (5.3) satisfies:

lim sup t→∞ E [V (x(t))] ≤ λ r (σ 2 + m 2 ) λ r + 2δn -2βnλ 1 (A) . (5.18)
Proof: We take the expectation of both sides of (5.17) and we obtain:

dE [V (x(t))] =E 2 n x T (t) (βA -δI) x(t) -βx T (t)X(t)Ax(t) dt + λ r n σ 2 + m 2 -E [V (x(t))] dt,
where we applied the Dominated Convergence Theorem to interchange the expectation and the derivative E [dV ] = dE [V ] and we used the fact that

E [∆V ] = 1 n (σ 2 + m 2 -E [V ]
) according to Proposition 5.2. Then we obtain the following ODE:

d dt E [V (x(t))] = E 2 n x T (t) (βA -δI) x(t) -βx T (t)X(t)Ax(t) + λ r n σ 2 + m 2 -E [V (x(t))] ≤ 2(βλ 1 (A) -δ)E [V (x(t))] + λ r n σ 2 + m 2 -E [V (x(t))] (5.19) = 2βnλ 1 (A) -2δn -λ r n E [V (x(t))] + λ r (σ 2 + m 2 ) n .
(5.20)

By the Comparison Lemma, we can guarantee that E [V (x(t))] is upper bounded by the solution of the right-hand side of (5.20). Then, it holds

E [V (x(t))] ≤ E [V 0 ] e 2βnλ 1 (A)-2δn-λr n t + e 2βnλ 1 (A)-2δn-λr n t t 0 λ r (σ 2 + m 2 ) n e - 2βnλ 1 (A)-2δn-λr n τ dτ = E [V 0 ] e 2βnλ 1 (A)-2δn-λr n t + λ r (σ 2 + m 2 ) λ r +2δn-2βnλ 1 (A) e 2βnλ 1 (A)-2δn-λr n t e - 2βnλ 1 (A)-2δn-λr n τ t 0 = E [V 0 ] - λ r (σ 2 + m 2 ) λ r + 2δn -2βnλ 1 (A) e 2βnλ 1 (A)-2δn-λr n t + λ r (σ 2 + m 2 ) λ r + 2δn -2βnλ 1 (A)
.

(5.21)

Finally, the result follows by taking the limit t → ∞ in (5.21). ■ Clearly, the bound (5.18) is a function of the parameters of the SIS epidemic (i.e., β, δ and λ 1 (A)) and the parameters of the replacement process (i.e., λ r , σ 2 and m 2 ). Similarly to (5.15), the bound is proportional to σ 2 + m 2 , but in this case the rate of the Poisson process λ r also appears as a multiplication factor. The upper bound for the rate of convergence of the SIS epidemic in a closed system determined by (5.7) appears in the denominator of (5.18), but this value is increased by the rate of the Poisson process λ r . From (5.21), it can be seen that a large value of λ r will increase the rate of convergence of V (x(t)), such that the convergence will be fast, but the asymptotic value will also increase. When λ r → ∞, the bound (5.18) is given by σ 2 + m 2 , which coincides with the result of Proposition 5.3 since the process will be characterized only by replacements.

Notice that the right-hand side in (5.19) corrresponds to the expectation of a SDE of the form:

dV = f (V )dt + g(V )dN t , (5.22)
which is also known as a Poisson Counter driven Stochastic Differential Equation [START_REF] Brockett | Stochastic control[END_REF][START_REF] Brockett | Stochastic analysis for fluid queueing systems[END_REF]. The solution of (5. [START_REF] Huynh | An integrated trust and reputation model for open multi-agent systems[END_REF]) is a cadlag function corresponding to a Piecewisedeterministic Markov Process (PDMP) and can also be studied using the theory of PDMP [START_REF] Jacobsen | Point Process Theory and Applications: Marked Point and Piecewise Deterministic Processes[END_REF]. Figure 5.5 presents the computation of V (x(t)) for the graph determined by (5.9) with n = 11 agents, β = 0.5, δ = 2 and λ r = 7. The states of the replaced agents are generated with a uniform distribution with m = 1/2 and σ 2 = 1/12. The expectation of V (x(t)) has been computed by considering 10000 realizations of the process.

Second order analysis

Although the asymptotic behavior of the expected aggregate function E [V (x(t))] provides useful information of a system subject to the replacements of agents, it does not give any information about the deviation of single realizations from the expected value.

dϕ(V ) = ⟨ dϕ dV , f (V )⟩dt + [ϕ(V + g(V )) -ϕ(V )] dN t ,
where f (V ) and g(V ) are the mappings in the SDE (5.22) and N t is a Poisson process.

We first analyze the change of E [V 2 (x(t))] during a replacement.

Lemma 5.2 During the replacement of an agent, the aggregate function V (x(t)) defined in (5.3) satisfies:

E V 2 (x(t + )) -E V 2 (x(t -)) ≤ - 2 n E V 2 (x(t -)) + E [Θ 4 ] n 2 + 2E [Θ 2 ] (n -1) + 1 n 2 E V (x(t -)) ,
where Θ is the new value of a replaced agent according to Assumption 5.2.

Proof: During the replacement of an agent x j , the function V 2 satisfies:

V 2 (x(t + )) = V (x(t -)) + Θ 2 n - x 2 j (t -) n 2 = V 2 (x(t -)) + Θ 4 n 2 + x 4 j (t -) n 2 + 2Θ 2 n V (x(t -)) - 2x 2 j (t -) n V (x(t -)) - 2Θ 2 x 2 j (t -) n 2 .
Then, we compute the conditional expectation:

E V 2 (x(t + ))|x(t -), Θ = 1 n n j=1 V 2 (x(t -)) + Θ 4 n 2 + x 4 j (t -) n 2 + 2Θ 2 n V (x(t -)) - 2x 2 j (t -) n V (x(t -)) - 2Θ 2 x 2 j (t -) n 2 = V 2 (x(t -)) + Θ 4 n 2 + 1 n 3 n j=1 x 4 j (t -) + 2Θ 2 n V (x(t -)) - 2V (x(t -)) n 2 n j=1 x 2 j (t -) - 2Θ 2 n 3 n j=1 x 2 j (t -) ≤ V 2 (x(t -)) + Θ 4 n 2 + 1 n 3 n j=1 x 2 j (t -) + 2Θ 2 n V (x(t -)) - 2V 2 (x(t -)) n - 2Θ 2 n 2 V (x(t -)) = V 2 (x(t -)) - 2V 2 (x(t -)) n + Θ 4 n 2 + 1 n 2 + 2Θ 2 n - 2Θ 2 n 2 V (x(t -))
Finally, we compute the total expectation to get the desired result. ■ Proposition 5.4 (Second moment) Consider a SIS epidemic satisfying the stability condition (5.2). Under Assumptions 5.1 and 5.2, the aggregate function V (x(t)) defined in (5.3) satisfies:

lim sup t→∞ E V 2 (x(t)) ≤ λ r E [Θ 4 ](λ r +2δn-2βnλ 1 (A))+λ 2 r E [Θ 2 ](2E [Θ 2 ](n-1)+1) 2n(λ r + 2δn -2βnλ 1 (A)) 2 .
(5.23)

Proof: By applying Itô Lemma with ϕ(V ) = V 2 , we obtain:

dV 2 (x(t)) ≤ ⟨2V (x(t)), 2(βλ 1 (A) -δ)V (x(t))⟩dt + V 2 (x(t + k )) -V 2 (x(t - k )) dN t = 4(βλ 1 (A) -δ)V 2 (x(t))dt + V 2 (x(t + k )) -V 2 (x(t - k )) dN t .
(5.24)

We take the expectation in (5.24) and we get:

dE V 2 (x(t)) ≤ 4(βλ 1 (A) -δ)E V 2 (x(t)) dt + E V 2 (x(t + )) -E V 2 (x(t -)) λ r dt = 4(βλ 1 (A) -δ)E V 2 (x(t)) dt + 2E [Θ 2 ] (n-1)+1 n 2 E [V (x(t))]+ E [Θ 4 ] n 2 - 2 n E V 2 (x(t)) λ r dt = 2 2βnλ 1 (A) -2δn -λ r n E V 2 (x(t)) dt + λ r 2E [Θ 2 ] (n -1) + 1 n 2 E [V (x(t))] dt + λ r E [Θ 4 ] n 2 dt,
where we use the Dominated Convergence Theorem to interchange the expectation and the derivative of V 2 (x(t)) and we apply Lemma 5.2. Let us denote

a V = 2 2βnλ 1 (A)-2δn-λr n , b V = λ r 2E[Θ 2 ](n-1)+1 n 2 and c V = λrE[Θ 4 ] n 2
, so that we have

d dt E V 2 (x(t)) ≤ a V E V 2 (x(t)) + b V E [V (x(t))] + c V (5.25)
The solution of (5.25) is given by:

E V 2 (x(t)) ≤ e a V t E V 2 0 + e a V t t 0 e -a V τ (b V E [V (x(τ ))] + c V ) dτ = e a V t E V 2 0 + b V e a V t t 0 e -a V τ E [V (x(τ ))] dτ + c V e a V t t 0 e -a V τ dτ = e a V t E V 2 0 + b V e a V t t 0 e -a V τ E [V (x(τ ))] dτ - c V a V 1 -e a V t .
(5.26)

We recall that by (5.21), the function E [V ] satisfies:

E [V (x(t))] ≤ d V e h V t + m V , where d V = E [V 0 ] -λr(σ 2 +m 2 ) λr+2δn-2βnλ 1 (A) , h V = 2βnλ 1 (A)-2δn-λr n t, m V = λr(σ 2 +m 2 ) λr+2δn-2βnλ 1 (A) . Then E [V 2 (x(t))] satisfies: E V 2 (x(t)) ≤ e a V t E V 2 0 - c V a V 1 -e a V t + b V e a V t t 0 e -a V τ d V e h V τ + m V dτ = e a V t E V 2 0 - c V a V 1 -e a V t + b V d V e a V t t 0 e (h V -a V )τ dτ + b V m V e a V t t 0 e -a V τ dτ = e a V t E V 2 0 - c V a V 1 -e a V t + b V d V h V -a V e h V t -e a V t - b V m V a V 1 -e a V t = E V 2 0 + c V a V - b V d V h V -a V + b V m V a V e a V t + b V d V h V -a V e h V t - c V a V + b V m V a V . (5.27)
Finally, the result follows by taking the limit t → ∞ in (5.27), where we recall that a V and h V are negative. ■ The bound (5.23) is derived using (5.18), which corresponds to the asymptotic behavior of E [V (x(t))]. Unlike the previous result (5.18), the value of (5.23) is also affected by the fourth moment of Θ, while (5.18) only depends on the second moment (E [Θ 2 ] = σ 2 + m 2 ). Proposition 5.4 can be used to derive an upper bound for the variance:

lim sup t→∞ Var(V (x(t))) = lim sup t→∞ E V 2 (x(t)) -(E [V (x(t))]) 2 ≤ lim sup t→∞ E V 2 (x(t)) .
Although the result is a bit conservative since we are neglecting the effect of -(E [V (x(t))]) 2 , we have an estimate of the deviation of the realizations from the expected value. However, a stronger bound can be obtained for the linearized model (5.10) by using the lower bound (5.11).

Corollary 5.1 (Variance) Consider the linearized model of a SIS epidemic satisfying the stability condition (5.2). Under Assumptions 5.1 and 5.2, the aggregate function V (x(t)) defined in (5.3) satisfies:

lim sup t→∞ Var(V (x(t))) ≤ λ r E [Θ 4 ](λ r +2δn-2βnλ 1 (A))+λ 2 r E [Θ 2 ](2E [Θ 2 ] (n-1)+1) 2n(λ r + 2δn -2βnλ 1 (A)) 2 - λ 2 r (E [Θ 2 ]) 2 (λ r + 2δn -2βnλ n (A)) 2 .
(5.28)

Proof: We compute the expectation of the aggregate funciton and we obtain

d dt E [V (x(t))] = E 2 n x T (t) (βA -δI) x(t) -βx T (t)X(t)Ax(t) + λ r n σ 2 + m 2 -E [V (x(t))] ≥ 2(βλ n (A) -δ)E [V (x(t))] + λ r n σ 2 + m 2 -E [V (x(t))] = 2βnλ n (A) -2δn -λ r n E [V (x(t))] + λ r (σ 2 + m 2 ) n , (5.29)
where we used (5.11) for the lower bound. By the Comparison Lemma, we can guarantee that E [V ] is lower bounded by the solution of the right-hand side of (5.29). We follow similar steps as in the proof of Theorem 5.1 and we obtain:

lim inf t→∞ E [V (x(t))] ≥ λ r (σ 2 + m 2 ) λ r + 2δn -2βnλ n (A) . (5.30)
The result follows by applying (5.23) and (5.30) in the definition of the variance:

Var(V (x(t))) = E V 2 (x(t)) -(E [V (x(t))]) 2 ,
and remembering that E [V 2 (x(t))] ≥ (E [V (x(t))]) 2 . ■ Notice that the right-hand side of (5.30) is 0 only when λ r = 0, so that this lower bound can be used to prove that the linearized model of a SIS epidemic under replacements is not stochastically stable in the sense that it does not satisfy [START_REF] Costa | Continuous-Time Markov Jump Linear Systems[END_REF]:

∞ 0 E [V (x(t))] dt < ∞.
In Figure 5.6, we present the computations of the second moments of V (x(t)) for the graph determined by (5.9) with n = 11 agents, β = 0.5, δ = 2 and λ r = 7 considering 10000 realizations of the process and a uniform distribution with m = 1/2 and σ 2 = 1/12 for the states of the replaced agents. where the dashed red line is the upper bound (5.23) and the dash-dotted yellow line corresponds to the upper bound (5.28) considering the linearized model of the SIS epidemic. Notice that Var(V (x(t))) does not go to zero with t and that the upper bound obtained using the linearized model is also valid for the nonlinear system for this particular network topology and parameters of the infection.

SIS epidemic with replacements on a graph sampled from graphon

From Proposition 5.1 and Theorem 5.1, it is clear that the analysis of the stability and behavior of a SIS epidemic in OMAS is linked to the largest eigenvalue λ 1 (A) of the adjacency matrix of the graph. However, in the case of large graphs, it is difficult to obtain a complete representation of the network because of the presence of noise and errors in data and the constant evolution of links and nodes. Also, even if it is possible to obtain a good knowledge of network topology, their sheer size prevents the full simulation or analysis of the dynamics, or the computation of relevant network properties, because of limitations in computational resources.

In this subsection, we consider that the graphs over which the epidemic takes place are sampled from graphons using the Definition A.1 and we derive results equivalent to Proposition 5.1 and Theorem 5.1 based on the spectrum of the graphon operator. (For a description of the most important properties of graphons see Section A.2). In the perspective approximating graphs with graphons, it is convenient to consider sequences of graphs parametrized by their size n. Therefore, it is reasonable to assume that parameters δ or β be also dependent on n: upon need, we shall emphasize this dependence by writing δ n and β n . Proposition 5.5 ([115]) Consider a SIS epidemic over a graph G with n nodes sampled from a piecewise Lipschitz graphon W . Then for n large enough, the SIS epidemic will reach the disease-free state with probability at least 1 -ν if:

δ n > nβ n (λ 1 (T W ) + ϕ n ), (5.31)
where λ 1 (T W ) is the largest eigenvalue of the graphon operator T W .

Proof: Condition (5.2) is equivalent to δ > βλ 1 (A). By Lemma A.1, with n large enough and probability at least 1 -ν:

|nλ 1 (A) -λ 1 (T W )| ≤ ϕ n , λ 1 (A) ≤ n(λ 1 (T W ) + ϕ n ),
which yields the desired result. ■ Since almost surely λ 1 (A) grows linearly with n (except for the trivial graphon W = 0), in order to ensure a uniform bound on λ 1 (A)β n /δ n , we will assume that the infection strength β n /δ n satisfies

β n δ n = 1 n β δ (5.32)
for some positive constants β and δ.

Assuming that δ n = δ is constant and β n = n -1 β means that, as the graph grows in size, the healing rate (which depends on each individual) remains constant, whereas the infection rate decreases. This natural scaling law is also chosen in [START_REF] Gao | Spectral representations of graphons in very large network systems control[END_REF]. Indeed, on dense graphs this assumption means that in larger graphs, even though there are more potential interactions, the average strength of the connections is suitably adjusted: this fact accounts for natural limitations in the rates of contact between individuals.

In the case of replacements, it is natural to assume that the frequency of these events increases with the number of agents in the system. For this reason, we assume:

λ r = n λr , (5.33)
for some positive constant λr .

Theorem 5.2 Consider a SIS epidemic over a graph G with n nodes sampled from a piecewise Lipschitz graphon W with δ n = δ and β n = n -1 β satisfying condition (5.31). Under Assumptions 5.1 and 5.2, for n large enough, with probability at least 1 -ν, the aggregate function V (x(t)) defined in (5.3) satisfies:

lim sup t→∞ E [V (x(t))] ≤ λr (σ 2 + m 2 ) λr + 2 δ -2 β(λ 1 (T W ) + ϕ n ) . ( 5

.34)

Proof: The aggregate function V (x(t)) satisfies:

dV (x(t)) = 2 n x T (t) (β n A -δ n I) x(t) -β n x T (t)X(t)Ax(t) dt + ∆V k dN t (5.35)
We compute the expectation and we obtain

d dt E [V (x(t))] = E 2 n x T (t) (β n A -δ n I) x(t) -β n x T (t)X(t)Ax(t) + λ r n σ 2 + m 2 -E [V (x(t))] ≤ E 2 n x T (t) (β n A -δ n I) x(t) + λ r n σ 2 + m 2 -E [V (x(t))] ≤ 2(β n λ 1 (A) -δ n )E [V (x(t))] + λ r n σ 2 + m 2 -E [V (x(t))] (5.36)
Then, with probability at least 1 -ν:

d dt E [V (x(t))] ≤ 2 (β n n(λ 1 (T W )+ϕ n )-δ n ) E [V (x(t))]+ λ r n σ 2 +m 2 -E [V (x(t))] = 2β n n 2 (λ 1 (T W ) + ϕ n ) -2δ n n -λ r n E [V (x(t))] + λ r (σ 2 + m 2 ) n , = 2 β(λ 1 (T W ) + ϕ n ) -2 δ -λr E [V (x(t))] + λr (σ 2 + m 2 ), (5.37)
which replaces 5.20. The rest of the proof follows the same step as in the proof of Theorem 5.1 ■ Remark 5.2 (Given graph) In the case of a given graph with a Poisson rate λ r given by (5.33), the upper bound (5.18) is equivalent to:

lim sup t→∞ E [V (x(t))] ≤ λr (σ 2 + m 2 ) λr + 2δ -2βλ 1 (A) ,
which is independent of the number of agents n. λr+2 δ-2 βλ 1 (T W ) . Hence, by choosing α > 1 and applying Borel-Cantelli Lemma, we obtain that for n large enough, with probability 1:

E [V (x(∞))] n ≤ λr (σ 2 + m 2 ) λr + 2 δ -2 β(λ 1 (T W ) + ϕ n ) → λr (σ 2 + m 2 ) λr + 2 δ -2 βλ 1 (T W ) as n → ∞, where E [V (x(∞))] n := lim sup t→∞ E [V (x(t))].
Remark 5.4 (Stochastic latent variables) An alternative procedure for the generation of the complete weighted graph Ḡ from graphons is the use of stochastic latent variables proposed in [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF]. All the results of this section easily extend to this stochastic sampling by considering the appropriate values of the probability ν and the bound ϕ n , which are established in [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF].

To illustrate the results on graphons, we consider the stochastic block model graphon W SB with pixel diagram in Fig. 5.7, where the values of the blocks are:

W SB =      
0.9 0.7 0.6 0.5 0.2 0.7 0.4 0.1 0.3 0.1 0.6 0.1 0.5 0.9 0.8 0.5 0.3 0.9 0.5 0.5 0.2 0.1 0.8 0.5 0.7

      .
The largest eigenvalue of the graphon operator T W SB is ρ W SB = 0.5275. We assume β n = β/n, with β = 1.5, and ν = 0.02, which satisfies condition (A.2c) for n ≥ 40. The value of the rate of replacements is λr = 2 while δ n = δ = β(λ 1 (T W ) + ϕ 40 ) is selected, satisfying condition (5.31). We compute the value of the upper bound E [V (x(∞))] n for 40 ≤ n ≤ 1000 and we present the results in Figure 5.8 where the dashed red line corresponds to the limit when n → ∞. ] n while the dashed red line corresponds to the limit when n → ∞.

SIS epidemic with replacements with new neighborhoods determined by graphons

Now, we consider the case of a SIS epidemic where during a replacement, the topology of the network may also change. The dynamics of the system is given by

ẋ(t) = (β n A(t) -δ n I) x(t) -β n X(t)A(t)x(t), (5.38)
where the adjacency matrix A(t) is time-varying.

We assume that at time t = 0, the graph of the network is sampled from a graphon W according to Definition A.1 and the network interconnections are preserved between the switching times determined by a Poisson process with intensity λ r , which implies that the dynamics (5.38) also remains invariant. However, we consider that during the replacement of an agent the topology of the network also changes.

Assumption 5.3 During the replacement of an agent, the connections of all the agents are sampled again according to the piecewise Lipschitz graphon W .

The motivation for this particular choice of the change of the network topology is that stochastic block model graphons are used to model community structures in social networks. In this case, the replacement of an agent may impact all the network such that connections between individuals may also change but preserve the same pattern determined by the graphon.

By definition, the aggregate function V (x(t)) satisfies:

dV (x(t)) = 2 n x T (t) (β n A(t) -δ n I) x(t) -β n x T (t)X(t)A(t)x(t) dt + ∆V k dN t ,
where A(t) is sampled from a graphon. Unlike Theorem 5.2 which is based on the computation of the expectation given a simple graph G sampled from a graphon, in this section, we compute the expectation of all the graphs that can form the sequence G t on which the process takes place. Then, the expectation is given by

d dt E [V (x(t))] = E 2 n x T (t) (β n A(t) -δ n I) x(t) -β n x T (t)X(t)A(t)x(t) + λ r n σ 2 + m 2 -E [V (x(t))] . (5.39)
We consider an approximation of (5.39) by assuming independence between the evolution of x(t) and the adjacency matrix A(t) such that we have:

d dt E [V (x(t))] = E 2 n x T (t) β n Ār -δ n I x(t) -β n x T (t)X(t) Ār x(t) + λ r n σ 2 + m 2 -E [V (x(t))] (5.40) ≤ E 2 n x T (t) β n Ār -δ n I x(t) + λ r n σ 2 + m 2 -E [V (x(t))] ≤ 2(β n λ 1 ( Ār ) -δ n )E [V (x(t))] + λ r n σ 2 + m 2 -E [V (x(t))] , (5.41)
where Ār = E [A(t)]. Notice that Ār = Ā -D, where Ā is the adjacency matrix of the complete weighted graph Ḡ according to Definition A.1 and D = diag[W (u i , u i )]. Then by applying Lemma A.1, we have that for n large enough:

λ 1 ( Ār ) = Ār ≤ Ā + D ≤ nλ 1 (T W ) + nϑ n + 1, where ϑ n = 2 √ L 2 -K 2 +Kn n
. Under this approximation, the stability condition for the SIS epidemic in expectation is given by:

δ n > nβ n (λ 1 (T W ) + ϑ n + 1/n).
(5.42)

By following similar steps in the proof of Theorem 5.2 with (5.41) instead of (5.36), we formulate the following conjecture.

Conjecture 5.1 Consider a SIS epidemic over a sequence of graphs G t with n nodes sampled from a piecewise Lipschitz graphon W with δ n = δ and β n = n -1 β satisfying condition (5.42). Under Assumptions 5.2 and 5.3, for n large enough, the aggregate function V (x(t)) defined in (5.3) satisfies:

lim sup t→∞ E [V (x(t))] ≤ λr (σ 2 + m 2 ) λr + 2 δ -2 β(λ 1 (T W ) + ϑ n + 1/n) . (5.43)
Theorem 5.2 holds with probability 1 -ν because of the upper bound for λ 1 (A) in terms of the spectral radius of the graphon operator T W , as per Lemma A.1. However, in the case of Conjecture 5.1, it can be seen from (5.40), that the potential result is based on an epidemic process with replacements over a time-invariant graph with adjacency matrix Ār , which is sampled from the graphon using deterministic latent variables and does not involve the use of random variables.

Figure 5.9 presents the computation of V (x(t)) for graphs sampled from the stochastic block model W SB with n = 50 agents, β = 1.5, δ = β(λ 1 (T W ) + ϑ 50 + 1/50) and λr = 2 where the states of new agents are drawn from a uniform distribution with m = 1/2 and σ 2 = 1/12, and the expected values were calculated considering 1000 realizations of the process. We can observe that the behavior of the approximated model (5.40) in dashed red line is similar to the behavior of the real model (5.39) in solid blue line, and both seem to converge to the same limit when t → ∞. Furthermore, the upper bound (5.43) in dash-dotted yellow line is valid for both processes, which supports the conjecture.

Conclusion

In this chapter, we analyzed a SIS epidemic in OMAS when replacements occur at time instants determined by a Poisson process. We defined an aggregate function to analyze the epidemic and derived upper bounds for the expectation and variance when t → ∞. We derived a similar result in expectation when graphs are sampled from graphons. When the connections of the network can change according to a graphon during the occurrence of replacements, we performed a preliminary analysis based on the assumption of independence of some random variables and formulated a conjecture for the asymptotic value in expectation of the aggregate function, which is corroborated by simulations.

Due to the versatility of graphons to generate different graph topologies with a given number of nodes, they are an important tool to be used in the case of OMAS over dynamic networks, which are more realistic, since it is unlikely that in a large network, during the replacement of an agent, the connections remain invariant. Graphons also appear to be a promising tool for analyzing OMAS under arrivals of agents, where the connections of the new agents can be well defined by generating edges based on the graphon and the values of the latent variables associated with the agents of the network. For future work, it would be important to extend the results by considering potential departures and arrivals of agents, which would imply a time-varying size of the network. In this case, the arrivals and departures of agents can be studied as birth and death events respectively and the theory of birth-death process [START_REF] Renshaw | Modelling Biological Populations in Space and Time[END_REF] can provide useful tools for the analysis of more general situations. This approach has been used in [START_REF] Monnoyer De Galland | Modelling gossip interactions in open multi-agent systems[END_REF], where a linear death process has been considered for the departure of agents. In this line of research, it is necessary to determine if the aggregate function used for replacement remains valid for possible changes on the number of agents of the system or if a new definition is needed to condensate the behavior of the system in a more appropriate way.

In the case of graphs sampled from graphons, we only considered the case of deterministic latent variables. However, the use of stochastic latent variables could provide a more appropriate sampling procedure during arrivals where the connections of the current agents in the network are preserved and new connections are established only with new agents.

Chapter 6 Conclusion and Future Work

Open multi-agent systems (OMAS) allow to model more realistic dynamic networks where agents can enter, leave or being replaced in the system. Unlike single changes in the network topology, these variations in the set of agents make the analysis more complex since the dimension of the system may change and traditional concepts associated to dynamical systems might not be applied in many cases. In this thesis, we analyzed OMAS considering different problems associated with classical multi-agent systems. In all of these problems, we selected a set of assumptions to make the problem realistic and tractable at the same time, and performed an analysis based on the study of an appropriate scalar quantity independent of the size of the system.

The problem of consensus was analyzed in Chapter 2. We considered a randomized algorithm where the interaction matrix is drawn from a common distribution and we studied the effect of additive noise in the performance of the algorithm. To perform the analysis, we considered the mean squared error and derived a closed form and bounds based on the eigenvalues of the expected matrices. This approach allowed us to study the problem in an OMAS when we consider a finite superset of agents, where at each time instant, only a subset of agents is active, determined by independent Bernoulli random variables associated with each agent in the system. We expressed the bounds derived for the general case as a function of the parameters associated with the OMAS and provided some expressions for particular graph topologies.

In Chapter 3, we focused on the weighted gradient descent algorithm to solve the well-known problem of resource allocation (RA) in multi-agent systems where the transmission of information is subject to packet losses. First, we proposed a set of assumptions to guarantee that the problem can be solved at least in expectation, which justifies the application of the algorithm. Then, we proposed two metrics to measure the deviation from the constraint and the ideal expected cost, and we derived bounds for the metrics which are proportional to the difference between 119 the initial cost function and the ideal cost function. Next, we extended the results to the case of OMAS, where the cost functions of the agents can be replaced with a certain probability. Unlike the closed system, we showed that the constraint violation error diverges and does not remain bounded due to the accumulation of errors. This result showed that a well-defined metric, which is bounded in a closed system, could diverge in an OMAS.

Chapter 4 continued the analysis of the RA problem in OMAS but considering the application of the RCD algorithm to perform the optimization. In this chapter we analyzed the problem from a perspective of dynamical systems, focusing on the stability of the time-varying minimizer. We proved that for a general graph topology, it is not possible to guarantee linear convergence in a closed system, and we performed an analysis based on a norm induced by a matrix associated with the network. Some bounds for the difference between the minimizers where proved using different approaches related to the properties of β-smooth and α-strongly convex functions. We proved that the RCD algorithm applied to an OMAS subject to replacements will converge inside a region determined by the parameters of the system. Additionally, by considering a simple setting, we proposed an alternative approach based on online optimization, where the states of the agents may also change during the replacements. We defined three metrics based on a cooperative and selfish strategy and derived some bounds for them which grow linearly with time. We provided a detailed analysis to justify this linear behavior, which is mainly due to the accumulation of errors during the replacements, similar to the constraint violation metric studied in Chapter 3.

In Chapter 5, we analyzed a SIS epidemic in OMAS where agents can be replaced during the evolution of the dynamics. Unlike the previous chapters, we used a continuous time formulation of the system, which made the analysis more complex since the states of the agents can change between replacements. Inspired by Lyapunov functions, we defined an appropriate aggregate function for the analysis of the system and used a stochastic differential equation (SDE) to model its evolution. We performed an analysis in expectation of the aggregate function and derived an upper bound for its variance. An extension of the results was provided when the graphs are sampled from graphons. In the case of possible changes of the connections during replacements following a sampling procedure from graphons, we presented a preliminary analysis of this scenario and formulated a conjecture about the asymptotic behavior of the aggregate function that was corroborated by simulations.

Future research

Although different types of problems and applications of OMAS were analyzed in this thesis, there are still many problems, questions, extensions and improvements that can be considered to develop a further understanding of the main characteristics of OMAS and the best tools to address more complex scenarios.

Regarding Chapter 2, the results were derived by considering independence in the stochastic matrices. However, scenarios where only few agents change their states from active to inactive or vice versa still need to be explored. In addition, it would be important to analyze an alternative performance metric where the normalization term includes only the number of active agents, which could provide a more informative result for OMAS.

In Chapter 3, considering that the upper bound derived for the cost function metric presents a certain level of conservatism, it would be important to improve this bound to obtain a tighter result. Regarding open systems, it would be important to improve the lower bound for the constraint violation metric. Finally, an analysis of the cost function metric for OMAS needs to be developed.

Chapter 4 analyzed the RCD algorithm in OMAS only in the case of possible replacements of agents while the budget remains the same. The case of timevarying budget and weights is a natural continuation of this work. A possible time-varying size of the system is also an interesting direction for future work, where the size of the system can change in time due to arrivals and departures of agents. In the case of the metrics inspired from online optimization, it would be important to extend the results for more general graph topologies and develop the corresponding lower bounds.

In Chapter 5, we analyzed a SIS epidemic in OMAS subject to replacements. It would be natural to extend the work to the case of scenarios where arrivals and departures can take place such that the size of the system can change in time [START_REF] Monnoyer De Galland | Modelling gossip interactions in open multi-agent systems[END_REF]. In the case of graphons, there is a conjecture that still needs to be verified to validate the results obtained in simulations. Also, for graphs sampled from graphons, it would be important to consider more complex sampling procedures, to study scenarios where the connections of the new agent do not follow a specific pattern and the connections of the other agents are preserved.

Finally, regarding OMAS in general, many questions and problems remain open. For instance, in this thesis we have considered only two approaches for the analysis of OMAS: time-invariant finite superset and fixed state-space. In the case of multi-mode multi-dimensional systems, preliminary works only considered the case of a fixed number of spaces with different dimensions [START_REF] Verriest | Multi-mode multi-dimensional systems[END_REF][START_REF]Pseudo-continuous multi-dimensional multi-mode systems[END_REF], which does not seem suitable for a general OMAS where the dimension of the system could take an arbitrary value. For this reason, further work is needed in this field to extend this approach for systems with arbitrary dimensions. Regarding infinite dimensional dynamical systems, during this thesis two works focused on the use of graphons for the study of epidemics [START_REF] Vizuete | Graphon-based sensitivity analysis of SIS epidemics[END_REF] and the analysis of the spectrum of the Laplacian matrix corresponding to sampled graphs [START_REF] Vizuete | The Laplacian spectrum of large graphs sampled from graphons[END_REF] were developed. Since most of the dynamics in multi-agent systems are formulated in terms of the Laplacian matrix, it is important to extend the theory of graphons for this matrix through the definition of an appropriate operator. In this line of research, the authors in [START_REF] Bonnet | Consensus formation in firstorder graphon models with time-varying topologies[END_REF] studied the case of consensus in first-order graphon models, while in [START_REF] Bramburger | Pattern formation in random networks using graphons[END_REF], the authors analyzed the spectrum of the Laplacian matrix for particular types of graphons. However, the relation between the spectrum of the Laplacian operator associated with an arbitrary graphon and the eigenvalues of the Laplacian matrix of a sampled graph is still under analysis.

In most of the works, OMAS have been analyzed only in the case of undirected graphs. In [START_REF] Xue | Stability of multi-dimensional switched systems with an application to open multi-agent systems[END_REF], the authors analyzed the case of switching between multi-agent systems of different dimensions where agents interact through digraphs. However, the approach used in this work disregards the impacts of the processes associated with arrivals, departures and replacements, limiting the potential applications of the work. Since most of the results derived form OMAS rely on the symmetry of the Laplacian matrix, it seems possible that new techniques of analysis are required for the extension of the results to the case of digraphs. In addition, the definition of appropriate graph topologies is critical in this case, specially in the case of arrivals where the connections of the new agents could play an important role in the preservation of the connectivity of the network.

Agents with more general dynamics like double integrators or nonlinear dynamics (Kuramoto oscillators, Lotka-Volterra, etc.) should be analyzed in the context of arrivals/departures of agents, which would extend the analysis of OMAS to more complex phenomena. In this line of research, the approach used in Chapter 5 seems promising for these cases, where Lyapunov functions could be used to study the systems in terms of a SDE. Alternatively, the theory of stochastic hybrid systems appear as an important framework that can be used for the analysis of OMAS when the objective is the determination of conditions to guarantee the stability of OMAS [START_REF] Teel | Stability analysis for stochastic hybrid systems: A survey[END_REF][START_REF] Teel | Stochastic hybrid systems: A modeling and stability theory tutorial[END_REF].

Heterogeneity is an important characteristic that should be explored in future works. In simple settings, it is assumed that all the agents that interact in the system have the same dynamics. However, in many scenarios, this assumption could not be appropriate as in the case of robotics, where different types of robots need to interact to achieve a desired goal [START_REF] Torres | On resilience and heterogeneity in robotic networks[END_REF]. Several approaches could be used for this analysis, including the study of emergent dynamics [START_REF] Panteley | Synchronization and dynamic consensus of heterogeneous networked systems[END_REF][START_REF] Maghenem | Singular-perturbations-based analysis of dynamic consensus in directed networks of heterogeneous nonlinear systems[END_REF], which has provided promising results in the case of switching systems with heterogeneous agents [START_REF] Adhikari | An emerging dynamics approach for synchronization of linear heterogeneous agents interconnected over switching topologies[END_REF].

Finally, as it was remarked in [START_REF] Proskurnikov | A tutorial on modeling and analysis of dynamic social networks. Part II[END_REF][START_REF] Ravazzi | Learning Hidden Influences in Large-Scale Dynamical Social Networks: A Data-Driven Sparsity-Based Approach, in Memory of Roberto Tempo[END_REF], OMAS look as a promising tool to handle opinion dynamics where agents are constantly entering and leaving the system as in the case of social networks, forums, discussions, etc. Models based on opinion-dependent connectivity like the Hegselmann-Krause model [START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF][START_REF] Blondel | Continuous-time average-preserving opinion dynamics with opinion-dependent communications[END_REF][START_REF] Ceragioli | Continuous and discontinuous opinion dynamics with bounded confidence[END_REF] and k-nearest neighbors [START_REF] Rossi | Asynchronous opinion dynamics on the k-nearestneighbors graph[END_REF] seem to be more suitable for OMAS where the interactions are determined only by the values of the opinions. This type of interconnections simplifies the definition of the network topology in OMAS in the case of arrivals and departures, since it is not necessary to define the connections of new agents or reconfigure the connections of the network when agents leave the system.

A.2 Graphons

We denote by W the space of all bounded symmetric measurable functions W : [0, 1] 2 → R. The elements of this space are called kernels given their connection to integral operators. The set of all kernels W ∈ W such that 0 ≤ W ≤ 1 is denoted by W 0 and its elements are called graphons, whose name is a contraction of graph-function. By analogy with degrees in finite graphs, the degree function of a graphon is defined as d W (x) := 1 0 W (x, y) dy. In order to consider differences between graphons, we shall sometimes work in the set W 1 of kernels W such that -1 ≤ W ≤ 1.

Every function W ∈ W defines an integral operator T W : L 2 [0, 1] → L 2 [0, 1] by:

(T W f ) (x) := 1 0 W (x, y)f (y) dy.
This operator is compact and has a discrete spectrum with 0 as the only accumulation point. Every nonzero eigenvalue has finite multiplicity [START_REF] Lovász | Large Networks and Graph Limits[END_REF]. A graphon W is said to have finite rank if the spectrum of the associated operator contains a finite number of nonzero eigenvalues [START_REF] Lovász | Large Networks and Graph Limits[END_REF].

A step graphon is a graphon defined as a step function. A function is called a step function if there is a partition S 1 ∪ • • • ∪ S k of [0, 1] into measurable sets such that W is constant on every product set S i × S j where the sets S i are the steps of W . This type of graphon is also called stochastic block model graphon because of its relation to stochastic block models [START_REF] Airoldi | Stochastic blockmodel approximation of a graphon: Theory and consistent estimation[END_REF].

Step graphons are finite rank graphons with a rank at most equal to the number of steps. Also, graphons expressed as a finite sum of products of integrable functions have finite rank [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF].

Each graph G has an associated step graphon W G obtained by considering a uniform partition of A graphon is usually visualized with a pixel picture, where each point (x, y) ∈ [0, 1] 2 is colored with a grey level representing W (x, y). For a step graphon asso-

A.2.1 Norms

In the study of kernels, various norms are relevant to consider [START_REF] Lovász | Large Networks and Graph Limits[END_REF][START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF][START_REF] Janson | Graphons, cut norm and distance, couplings and rearrangements[END_REF]. For 1 ≤ p < ∞, we define the L p norm of a kernel as ∥W ∥ L p := For W ∈ W 1 , we have the following inequalities between L p norms and the cut norm:

∥W ∥ □ ≤ ∥W ∥ L 1 ≤ ∥W ∥ L 2 ≤ ∥W ∥ 1/2 L 1 ≤ 1.
By considering the operator T W associated to a kernel W ∈ W, we can define the operator norm:

|||T W ||| := sup f ∈L 2 [0,1] ∥f ∥ L 2 =1 ∥T W f ∥ L 2 .
For graphons, the operator norm is equal to the largest eigenvalue of the operator:

|||T W ||| = λ 1 (T W ).
For the elements of W 1 , the cut and operator norms are related by:

∥W ∥ □ ≤ |||T W ||| ≤ √ 8∥W ∥ 1/2 □ .

A.2.2 Sampling and Approximation

A graphon W can be used to generate random graphs using a sampling method [START_REF] Lovász | Large Networks and Graph Limits[END_REF].

Definition A.1 (Sampled Graph [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF]) Given a graphon W and a size n ∈ Z >0 , we say that the graph G is sampled from W if it is obtained through:

1. Complete weighted graph Ḡ: fixing deterministic latent variables {u i = i n } n i=1 , we generate the complete weighted graph Ḡ with n vertices, whose adjacency matrix is defined as: Ā(i, j) = W (u i , u j ) for all i, j in {1, . . . , n}.

2.

Simple graph G: taking n vertices {1, . . . , n} and randomly adding undirected edges between vertices i and j independently with probability W (u i , u j ) for all i > j.

Definition A.2 (Piecewise Lipschitz graphon [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF]) Graphon W is said to be piecewise Lipschitz if there exists a constant L and a sequence of non-overlapping intervals I k = [α k-1 , α k ) defined by 0 = α 0 < • • • < α K+1 = 1, for a finite

A.3 Poisson processes

A stochastic process (N t ) t≥0 is a counting process if N t represents the total number of events that have occurred up to time t. More specifically, a stochastic process is said to be a counting process if [START_REF] Ross | Stochastic Processes[END_REF]:

1. N t ≥ 0 2. N t ∈ Z ≥0
3. If s < t, then N s ≤ N t 4. For s < t, N t -N s equals the number of events that have occurred in the interval (s, t]

A counting process has independent increments if the numbers of events that occur in disjoint time intervals are independent. A counting process has stationary increments if the distribution of the number of events that occur in any interval of time depends only on the length of the time interval [START_REF] Ross | Stochastic Processes[END_REF].

Definition A.4 (Poisson process [START_REF] Ross | Stochastic Processes[END_REF]) The counting process (N t ) t≥0 is said to be a homogeneous Poisson process with rate λ > 0 if 1. N 0 = 0 2. The process has independent increments 3. The number of events in any interval of length t is Poisson distributed with mean λt. That is, for all s, t ≥ 0, P {N t+s -N s = n} = e -λt (λt) n n! , n = 0, 1,

The last property implies that a Poisson process has stationary increments and that E [N t ] = λt. By construction, a homogeneous Poisson process only has a finite number of jumps on any interval and is a Markov counting process (i.e., has the Markov property), as the exponential distribution has the memoryless property. Abstract: In this thesis we address several problems formulated in an open multi-agent system (OMAS) scenario where the set of agents can change in time, independently of the evolution of the dynamics associated with the system. We analyze the case of OMAS formulated in a fixed finite size network and we use two approaches for the analysis of the systems. In the first approach we consider scenarios characterized by activation/deactivation of agents such that at each time instant a different subset of active agents can interact in the system. In the second approach we study scenarios characterized by replacements of agents where at a specific time instant, an agent can be replaced while the rest of the agents remain the same. In this case, all the agents are able to interact at all time.

Three different problems are considered in this thesis: randomized consensus, resource allocation problem and epidemics. First, we analyze the problem of randomized consensus subject to additive noise where different subset of agents exchange information at each iteration. We define a noise index based on the expected mean squared error and we derive upper bounds. Then, we consider the resource allocation problem where agents can be replaced during the implementation of an optimization algorithm. For this problem, we analyze two different algorithms: weighted gradient descent and random coordinate descent. For the weighted gradient descent, we evaluate the performance of the algorithm in an OMAS subject to packet losses by defining appropriate performance metrics. For the random coordinate descent algorithm, we study the convergence to the minimizer in an OMAS and we provide an alternative analysis using tools inspired from online optimization. Finally, we study a SIS epidemic in continuous time subject to replacements of agents during the evolution of the disease. We perform the analysis using an aggregate function and deriving upper bounds for its asymptotic behavior.

3 . 3 47 3. 4 48 4. 1 63 4. 3 75 4. 4 76 LIST OF FIGURES xv 4 . 5 117 A. 1

 3347448163375476451171 and the relative error for a sequence of graphs sampled from the graphon W (x, y) = 1 -xy with growing n. . . . . 2.8 Computations of the relative error for different values of p for graphs with n = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Computation of the upper bound for J constr for a complete, star and line graph with n = 10, β = 1, ϵ = 0.5 and f (x 0 ) -f * = 1. . . . . . 3.2 Computations of the upper bound for J funct for a complete, star and line graph with n = 10, ϵ = 0.5, f (x 0 ) -f * = 1 and two values of κ. Simulation of E (1 T x(k) -b) 2 and of the lower bound in Proposition 3.2 for a complete graph under replacements with n = 7, α = 1, β = 10, p = 0.8, p U = 0.5 by considering 10000 realizations of the process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Simulation of E (1 T x(k) -b) 2 and of the lower bound in Proposition 3.2 for a complete graph with no replacements with n = 7, α = 1, β = 10, p = 0.8, p U = 1 by considering 10000 realizations of the process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Bound of Proposition 4.1 with respect to the system size n for b = 1, c = 1, respectively for κ = 50 with homogeneous agents (a i = 1 for all i) on the left, and κ = 2 with heterogeneous agents (a 1 = 10, a i = 1 for i > 1) on the right. The figures show all three quantities ψ n,κ , θ n,κ and χ n,κ as well as the final bound M 2 n,κ for both cases. . . 61 4.2 Evolution of the upper bound on x (1) -x (2) 2 respectively with respect to n for several values of κ (top) and with respect to κ for several values of n (bottom). For each plot the bound obtained in Proposition 4.1 (right) is compared with the empirical upper bound derived using PESTO in the same settings (left). The top-right plot also shows the asymptotic value expected to be reached by θ n,κ as n → ∞ based on (4.16). . . . . . . . . . . . . . . . . . . . . . . . . Performance of the RCD algorithm in a complete graph constituted of respectively n = 5 agents with κ = 5 (left) and n = 30 agents with κ = 1.2 (right), with p U = 0.95 and b = 1, and where each local objective function is defined by (4.57). The solid blue and red dashed lines represent the actual performance of the algorithm averaged over 500 realizations of the process, respectively for the random and adversarial replacements cases. The yellow dotted line is the upper bound (4.42) obtained from Corollary 4.4. . . . . . . . Performance of the RCD algorithm with n = 5 agents, κ = 1.2, p U = 0.95 and b = 1, respectively in (left) a ring graph with homogeneous agents (i.e., a i = 1 for all i) and (right) a complete graph with heterogeneous agents (i.e., a 1 = 10, a i = 1 for i > 1), and where each local objective function is defined by (4.57). The solid blue and red dashed lines represent the actual performance of the algorithm averaged over 500 realizations of the process, respectively for the random and adversarial replacements cases. The yellow dotted line is the upper bound (4.42) obtained from Corollary 4.4. . . . . . . . Evolution of the function value f k evaluated with the RCD algorithm x(k) defined in (4.64), the optimal solution x * (k), and the selfish strategy x s (k), in a system subject to replacements of agents (i.e., simultaneous departures and arrivals) on average once every 4 RCD steps. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Evolution of the averaged asymptotic expected Potential Benefit (on the left) and dynamical regret (on the right) in a system of 5 agents with ρ R = 0.0125 and κ = 10. Each plot compares the upper bounds, respectively from (4.70) and (4.84), with simulated results, either with random replacements (RR) or adversarial replacements (AR). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 Evolution of the expected averaged regret for a system of 5 agents holding quadratic functions with ρ R = 0.0125 and κ = 10. The solid blue line and the dash-dotted red line respectively correspond to the upper bounds of Theorems 4.3 and 4.4 respectively. The dotted yellow line corresponds to simulation where we consider random replacements (RR). . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Graph G. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Evolution of the aggregate function V (x(t)) defined in (5.3) in continuous time for the graph determined by (5.9) with n = 11 agents, β = 0.5 and δ = 5.3. The solid blue line corresponds to the simulation of V (x(t)) while the dashed red line is the upper bound (5.8). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Sample path of V (x(t)). . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Evolution of the aggregate function V (x(t)) defined in (5.3) in a system subject only to replacements with λ r = 30. In the right plot, the solid blue line corresponds to the computation of V (x(t)) considering 10000 realizations of the process while the dashed red line is (5.14). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.5 Evolution of the aggregate function V (x(t)) defined in (5.3) for the graph determined by (5.9) with n = 11 agents, β = 0.5, δ = 2 and λ r = 7. In the right plot, the solid blue line corresponds to the computation of E [V (x(t))] considering 10000 realizations of the process while the dashed red line is the upper bound (5.18). . . . . 5.6 Evolution of the second moments of the aggregate function V (x(t)) defined in (5.3) for the graph determined by (5.9) with n = 11 agents, β = 0.5, δ = 2 and λ r = 7. In the left plot, the solid blue line corresponds to the estimation of E [V 2 (x(t))] while the dashed red line is the upper bound (5.23). In the right plot, the solid blue line corresponds to the estimation of Var(V (x(t))) while the dashed red line is the upper bound (5.23) and the dash-dotted yellow line is the upper bound (5.28) considering the linearized model. The simulated values were obtained considering 10000 realizations of the process. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 110 5.7 Pixel diagram of the stochastic block model graphon. . . . . . . . . 113 5.8 Upper bound E [V (x(∞))] n as a function of n for the stochastic block model graphon W SB with β = 1.5, δ = β(λ 1 (T W ) + ϕ 40 ) and λr = 2. The solid blue line corresponds to E [V (x(∞))] n while the dashed red line corresponds to the limit when n → ∞. . . . . . . . 114 5.9 Evolution of the aggregate function V (x(t)) defined in (5.3) for graphs sampled from the stochastic block model graphon W SB with n = 50 agents, β = 1.5, δ = β(λ 1 (T W )+ϑ 50 +1/50) and λr = 2. The solid blue line corresponds to the estimation of E [V (x(t))] while the dashed red line corresponds to the approximation (5.40) of the real process and the dash-dotted yellow line is the conjectured upper bound (5.43). The expected values were computed by considering 1000 realizations of the process. . . . . . . . . . . . . . . . . . . . . Graph G, adjacency matrix A and step graphon W G . . . . . . . . . 127 A.2 Sampling graphon procedure. . . . . . . . . . . . . . . . . . . . . . 128

Figure 1 . 1 :

 11 Figure 1.1: Replacement of an agent in OMAS, while preserving the topology.

Figure 1 . 2 :

 12 Figure 1.2: Time-invariant finite superset of agents characterized by different active agents.

Figure 1 . 3 :

 13 Figure 1.3: Multi-mode multi-dimensional systems as a tool to study OMAS.

Figure 2 . 1 :

 21 Figure 2.1: Random interactions with different numbers of active agents n a in a network with n = 18 possible agents.

Figure 2 . 2 :

 22 Figure 2.2: Computation of J noise and the relative error for a star graph S n with growing n. In the left plot, the solid blue line corresponds to the lower bound (2.15) and the dashed red line is the upper bound (2.16).

Fig. 2 .

 2 3(a) and 2.3(b) present the results of the computations for the line graph. The maximum degree for a 2-D grid is d max (2-D) = 4 and for a 3-D grid is d max (3-D) = 6. The spectrum of the Laplacian matrix for a m-D grid are given by: (a) J noise (b) Relative error

Figure 2 . 3 :

 23 Figure 2.3: Computation of J noise and the relative error for a line graph P n with growing n. In the left plot, the solid blue line corresponds to the lower bound (2.17) and the dashed red line is the upper bound (2.18).

Figure 2 . 4 :

 24 Figure 2.4: Computation of J noise for 2-D and 3-D grids with growing n.

  (a) and 2.5(b): the latter implies that the upper bound is equal to the noise index.

  (a) J noise (b) Relative error

Figure 2 . 5 :

 25 Figure 2.5: Computation of J noise and the relative error for a complete graph K n with growing n. In the left plot, the solid blue line corresponds to the lower bound (2.19) and the dashed red line is the upper bound (2.20).

Figure 2 . 6 :

 26 Figure 2.6: Computation of J noise and the relative error for a sequence of Erdős-Rényi graphs with growing n.

Figure 2 . 7 :

 27 Figure 2.7: Computation of J noise and the relative error for a sequence of graphs sampled from the graphon W (x, y) = 1 -xy with growing n.

Figure 2 . 8 :

 28 Figure 2.8: Computations of the relative error for different values of p for graphs with n = 100.

Assumption 3 .

 3 1 provides an upper and a lower bound to the curvature of the functions. The value κ := β α ≥ 1 is called the condition number of the functions. The set of the functions satisfying Assumption 3.1 is denoted by F α,β .

Assumption 3 .

 3 2 (1-D functions) The local cost function of any agent at any time is one-dimensional: f t i : R → R. Assumption 3.3 (Fixed demand) The total demand associated with the system at any time is n i=1 d i = b. Assumption 3.4 (Homogeneous agents) The weight associated with any agent i at any time is a i = 1.

Proposition 3 . 1

 31 Under Assumptions 3.1 to 3.4, for any initial point x 0 ∈ S 1,b , the sequence {x(k), k > 0}, produced by (3.6) with h ≤ 2 βλn(L) , satisfies:

Definition 3 . 1 (

 31 Piecewise quadratic function) A continuous differentiable function f : R → R is piecewise quadratic if there exists finitely many quadratic functions {q i } I i=1 such that f (x) ∈ {q i (x)} I i=1 for all x ∈ R. Assumption 3.5 (Piecewise quadratic functions) The local cost function of any agent at any time is piecewise quadratic.

Lemma 3 . 1

 31 For L(k) given by (3.10) and (3.11):

Figure 3 . 1 :

 31 Figure 3.1: Computation of the upper bound for J constr for a complete, star and line graph with n = 10, β = 1, ϵ = 0.5 and f (x 0 ) -f * = 1.

. 26 )

 26 From(3.26), we can observe that large values of ϵ, corresponding to large values of the step-size, increase the potential violation of the constraint. Figure3.1 presents the computation of the upper bound for the constraint violation metric in a complete, star and line graph with n = 10 agents, β = 1 and ϵ = 0.5, where we can observe the behavior of the upper bound and the limit values determined by(3.25) and(3.26).

Figure 3 . 2 :

 32 Figure 3.2: Computations of the upper bound for J funct for a complete, star and line graph with n = 10, ϵ = 0.5, f (x 0 ) -f * = 1 and two values of κ.

Figure 3 . 3 :

 33 Figure 3.3: Simulation of E (1 T x(k) -b) 2 and of the lower bound in Proposition 3.2 for a complete graph under replacements with n = 7, α = 1, β = 10, p = 0.8, p U = 0.5 by considering 10000 realizations of the process.

Figure 3 . 4 :

 34 Figure 3.4: Simulation of E (1 T x(k) -b) 2 and of the lower bound in Proposition 3.2 for a complete graph with no replacements with n = 7, α = 1, β = 10, p = 0.8, p U = 1 by considering 10000 realizations of the process.

Lemma 4 . 1 8 )

 418 Let x * := arg min x∈S a,b n i=1 f i (x i ). If all f i satisfy Assumptions 3.1 and 4.1, then x * ∈ B(0 nd , R b,κ ) with R b,κ = √ where c was defined in Assumption 4.1. Proof: We prove by contradiction. Let us consider some x ∈ S a,b such that x / ∈ B(0, R b,κ ), and let x * = arg min x∈R nd f (x) denote the global minimizer without constraint. We have ||x|| > R b,κ by definition and ||x * || ≤ √ nc since x * ∈ B(0 d , c) n , and it thus follows that ||x -x * || > R b,κ -√ nc. Hence, since f is α-strongly convex, and reminding that f (x * ) = 0 from Assumption 4.1, there holds

2 .

 2 Now let us define x b := 1 ||a|| 2 (a ⊗ I d )b, which is a feasible point with norm ||x b || = ||b|| ||a|| . Since f is β-smooth, and since f (x * ) = 0 from Assumption 4.1, there holds

  1 and 4.1 for all i = 1, . . . , n, then for ζ * defined in (3.3) there holds ||ζ * || ≤ β ||a|| 2 (||b|| + c ||a|| 1 ) ; (4.9) and ||x * i || ≤ a i ||a|| 2 κ (||b|| + c ||a|| 1 ) + c. (4.10)

Figure 4 . 1 :

 41 Figure 4.1: Bound of Proposition 4.1 with respect to the system size n for b = 1, c = 1, respectively for κ = 50 with homogeneous agents (a i = 1 for all i) on the left, and κ = 2 with heterogeneous agents (a 1 = 10, a i = 1 for i > 1) on the right. The figures show all three quantities ψ n,κ , θ n,κ and χ n,κ as well as the final bound M 2 n,κ for both cases.

Figure 4 . 2 :

 42 Figure 4.2: Evolution of the upper bound on x (1) -x (2) 2 respectively with respect to n for several values of κ (top) and with respect to κ for several values of n (bottom). For each plot the bound obtained in Proposition 4.1 (right) is compared with the empirical upper bound derived using PESTO in the same settings (left). The top-right plot also shows the asymptotic value expected to be reached by θ n,κ as n → ∞ based on (4.16).

Proposition 4 . 2

 42 The seminorm ∥•∥ L † p is a norm on S a,0 .

Proposition 4 . 3

 43 Let f (x) := n i=1 f i (x i ) and x * := arg min x∈S a,b f (x). Under Assumption 3.1, for any positive scalar h ≤ λ 2 (L p ) λ n (L p )

Remark 4 . 5 (

 45 Alternative rate) Starting from(4.35) in the proof of Proposition 4.3, one can use a similar argument to derive the following alternative convergence rate, valid for λ 2(Lp)

Proposition 4 . 4

 44 Let a function f (x) := n i=1 f i (x i ) and x * := arg min x∈S b f (x). Under Assumption 3.1, for any positive scalar

2 α+β

 2 3 the step-size must satisfy h ≤ µ 2 µn while for Proposition 4.4 the step-size is upper bounded by h ≤ 4 µn 2 α+β .

Theorem 4 . 1

 41 In the setting described in Section 4.1, the iteration rule (4.24) with h ≤ λ 2 (Lp) λn(Lp) 2 α+β generates a sequence of estimates x(k) satisfying the following for all k and for any η > 0:

Corollary 4 . 4

 44 When η = η * = η 1 + 1 + Mn,κ η , then the convergence of the RCD algorithm in open system is guaranteed with minimal asymptotic error Γ η * , and there holds:

  among 100 realizations of such uniform random choice.In Fig.4.3, we show the evolution of the expected errorE ∥x(k) -x * (k)∥ 2 L †p simulated for a network with interconnections defined by a complete graph, homogeneous agents and uniform probabilities. We consider two parametrizations of κ and n in both random and adversarial replacement cases, and we compare the simulations with (4.42) using the values given by Corollary 4.4. The figure shows that convergence is indeed guaranteed for the RCD in the presented settings and that the result of Corollary 4.4 shows some conservatism, which is inherited from Proposition 4.1. It is interesting to point out that these settings respectively make use of M 2 n,κ = θ n,κ for n = 30, κ = 1.2, and of M 2 n,κ = ψ n,κ for n = 5, κ = 5, consistently with the description of M 2 n,κ in the homogeneous case of Section 4.2.

Figure 4 . 3 :

 43 Figure 4.3: Performance of the RCD algorithm in a complete graph constituted of respectively n = 5 agents with κ = 5 (left) and n = 30 agents with κ = 1.2 (right), with p U = 0.95 and b = 1, and where each local objective function is defined by (4.57). The solid blue and red dashed lines represent the actual performance of the algorithm averaged over 500 realizations of the process, respectively for the random and adversarial replacements cases. The yellow dotted line is the upper bound (4.42) obtained from Corollary 4.4.

4

 4

  .4 are qualitatively similar to the case of the complete graph with homogeneous agents presented in Fig. 4.3.

Figure 4 . 4 :

 44 Figure 4.4: Performance of the RCD algorithm with n = 5 agents, κ = 1.2, p U = 0.95 and b = 1, respectively in (left) a ring graph with homogeneous agents (i.e., a i = 1 for all i) and (right) a complete graph with heterogeneous agents (i.e., a 1 = 10, a i = 1 for i > 1), and where each local objective function is defined by (4.57). The solid blue and red dashed lines represent the actual performance of the algorithm averaged over 500 realizations of the process, respectively for the random and adversarial replacements cases. The yellow dotted line is the upper bound (4.42) obtained from Corollary 4.4.

Assumption 4 . 2 (Assumption 4 . 4 (Assumption 4 . 5

 424445 Nonnegative 1-D functions) The local cost function of any agent at any time is nonnegative one-dimensional: f k i : R ≥0 → R ≥0 . Assumption 4.3 (Minimizer) The local cost function of any agent satisfies arg min x∈R p ≥0 f i (x) = 0 p and f i (0 p ) = 0. Homogeneous demand) The demand associated with any agent i at any time k is d k i = 1 such that b = n. The graph G = (V, E) is complete.

Figure 4 . 5 :

 45 Figure 4.5: Evolution of the function value f k evaluated with the RCD algorithmx(k) defined in (4.64), the optimal solution x * (k), and the selfish strategy x s (k), in a system subject to replacements of agents (i.e., simultaneous departures and arrivals) on average once every 4 RCD steps.

Proposition 4 . 5 (

 45 Well-posedness) The event set (4.2) guarantees that if x 0 ∈ S 1,n , then x(k) ∈ S 1,n for all k.

Figure 4 . 6 :

 46 Figure 4.6: Evolution of the averaged asymptotic expected Potential Benefit (on the left) and dynamical regret (on the right) in a system of 5 agents with ρ R = 0.0125 and κ = 10. Each plot compares the upper bounds, respectively from (4.70) and (4.84), with simulated results, either with random replacements (RR) or adversarial replacements (AR).

Fig. 4 .Figure 4 . 7 :

 447 Figure 4.7: Evolution of the expected averaged regret for a system of 5 agents holding quadratic functions with ρ R = 0.0125 and κ = 10. The solid blue line and the dash-dotted red line respectively correspond to the upper bounds of Theorems 4.3 and 4.4 respectively. The dotted yellow line corresponds to simulation where we consider random replacements (RR).

Figure 5 . 1 :

 51 Figure 5.1: Graph G.

Figure 5 . 2 :

 52 Figure 5.2: Evolution of the aggregate function V (x(t)) defined in (5.3) in continuous time for the graph determined by (5.9) with n = 11 agents, β = 0.5 and δ = 5.3. The solid blue line corresponds to the simulation of V (x(t)) while the dashed red line is the upper bound (5.8).

Figure 5 . 3 :

 53 Figure 5.3: Sample path of V (x(t)).

. 15 )

 15 which corresponds to the second moment of the distribution used to generate the new values of the replaced agents. The rate of the Poisson process λ r only has an influence in the rate of convergence of V (x(t)) such that a large value of λ r will

Figure 5 . 4 :

 54 Figure 5.4: Evolution of the aggregate function V (x(t)) defined in (5.3) in a system subject only to replacements with λ r = 30. In the right plot, the solid blue line corresponds to the computation of V (x(t)) considering 10000 realizations of the process while the dashed red line is(5.14).

Figure 5 .

 5 4(b) presents the evolution of E [V (x(t))] computed considering 10000 realizations of the process.

Figure 5 . 5 :

 55 Figure 5.5: Evolution of the aggregate function V (x(t)) defined in (5.3) for the graph determined by (5.9) with n = 11 agents, β = 0.5, δ = 2 and λ r = 7. In the right plot, the solid blue line corresponds to the computation of E [V (x(t))] considering 10000 realizations of the process while the dashed red line is the upper bound (5.18).

Figure 5 .

 5 6(a) shows the evolution of E [V 2 (x(t))] where the dashed red line corresponds to the upper bound(5.23). Finally, Fig.5.6(b) presents the evolution of Var(V (x(t)))

Figure 5 . 6 :

 56 Figure 5.6: Evolution of the second moments of the aggregate function V (x(t)) defined in(5.3) for the graph determined by (5.9) with n = 11 agents, β = 0.5, δ = 2 and λ r = 7. In the left plot, the solid blue line corresponds to the estimation of E [V 2 (x(t))] while the dashed red line is the upper bound (5.23). In the right plot, the solid blue line corresponds to the estimation of Var(V (x(t))) while the dashed red line is the upper bound (5.23) and the dash-dotted yellow line is the upper bound (5.28) considering the linearized model. The simulated values were obtained considering 10000 realizations of the process.

Figure 5 . 7 :

 57 Figure 5.7: Pixel diagram of the stochastic block model graphon.

Figure 5 . 8 :

 58 Figure 5.8: Upper bound E [V (x(∞))] n as a function of n for the stochastic block model graphon W SB with β = 1.5, δ = β(λ 1 (T W ) + ϕ 40 ) and λr = 2. The solid blue line corresponds to E [V (x(∞))] n while the dashed red line corresponds to the limit when n → ∞.

Figure 5 . 9 :

 59 Figure 5.9: Evolution of the aggregate function V (x(t)) defined in (5.3) for graphs sampled from the stochastic block model graphon W SB with n = 50 agents, β = 1.5, δ = β(λ 1 (T W ) + ϑ 50 + 1/50) and λr = 2. The solid blue line corresponds to the estimation of E [V (x(t))] while the dashed red line corresponds to the approximation (5.40) of the real process and the dash-dotted yellow line is the conjectured upper bound(5.43). The expected values were computed by considering 1000 realizations of the process.

af

  [0, 1] into the intervals B n i , whereB n i = [(i -1)/n, i/n) for i = 1, . . . , n -1 and B n n = [(n -1)/n, 1] such that: W G (x, y) := ij 1 B n i (x)1 B n j (y),where 1 A (x) is the indicator function. The operator associated to the step graphon is(T W G f )(x) := (y) dy for any x ∈ B n iand the spectrum of T W G consists of the normalized spectrum of the graph (i.e., λ i (T W G ) = λ i (A)/n), together with infinitely many zeros.

Figure A. 1 :

 1 Figure A.1: Graph G, adjacency matrix A and step graphon W G .

[0, 1 ] 2 |W

 12 (x, y)| p dx dy 1/p and its cut norm by ∥W ∥ □ := sup S,T ⊆[0,1] S×T W (x, y) dx dy .

Figure A. 2 :

 2 Figure A.2: Sampling graphon procedure.

Titre:

  Contributions aux systèmes multi-agents ouverts : consensus, optimisation et épidémies Mots clés: Systèmes multi-agents ouverts, optimisation distribuée, consensus, agents et systèmes autonomes Résumé: Dans cette thèse, nous abordons différents problèmes formulés dans un scénario de système multi-agents ouvert (OMAS) où l'ensemble des agents peut changer dans le temps, indépendamment de l'évolution de la dynamique associée au système. Nous considérons des OMAS formulés à l'aide d'un réseau fixe de taille finie, et utilisons deux approches distinctes pour les analyser. Dans la première approche, nous considérons des scénarios caractérisés par l'activation et la désactivation d'agents, de sorte qu'à chaque instant un sous-ensemble différent d'agents actifs peut interagir dans le système. Dans la deuxième approche, nous étudions des scénarios caractérisés par des remplacements d'agents où, à un instant donné, un agent peut être remplacé alors que les autres ne changent pas. Dans ce cas, tous les agents peuvent interagir à tout moment. Trois problèmes différents sont considérés dans cette thèse : le consensus randomisé, le problème d'allocation des ressources et les épidémies. Premièrement, nous analysons le problème du consensus randomisé soumis à un bruit additif où différents sous-ensembles d'agents échangent des informations à chaque itération. Nous définis-sons un indice de bruit basé sur l'erreur quadratique moyenne attendue et nous dérivons des bornes supérieures. Ensuite, nous considérons le problème d'allocation de ressources où les agents peuvent être remplacés lors de l'implémentation d'un algorithme d'optimisation. Pour ce problème d'optimisation, nous analysons deux algorithmes différents : la descente de gradient pondérée (weighted gradient descent) et la descente de coordonnées aléatoires (random coordinate descent). Pour la descente de gradient pondérée, nous évaluons les performances de l'algorithme dans un OMAS soumis à des pertes de paquets en définissant des métriques de performances appropriées. Pour l'algorithme de descente de coordonnées aléatoires, nous étudions la convergence vers le minimiseur dans un OMAS et nous proposons une analyse alternative à l'aide d'outils inspirés de l'optimisation en ligne. Enfin, nous étudions une épidémie SIS en temps continu sujette à des remplacements d'agents lors de sa propagation. Nous effectuons l'analyse en utilisant une fonction d'agrégation et en dérivant des bornes supérieures pour son comportement asymptotique. Title: Contributions to open multi-agent systems: consensus, optimization and epidemics Keywords: Open multi-agent systems, distributed optimization, consensus, agents and autonomous systems

  . Main contributions of this thesis . . . . . . . . . . . . . . . . . . . . 1.4 List of publications . . . . . . . . . . . . . . . . . . . . . . . . . . . Randomized consensus with additive noise . . . . . . . . . . . . . . 2.1.1 Upper and lower bounds on the noise index . . . . . . . . . 2.2 Noise index in OMAS . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Particular graph topologies . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Sparse graphs . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.2 Dense graphs . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Weighted gradient descent with packet losses and agent replacements 3.1 Resource Allocation Problem . . . . . . . . . . . . . . . . . . . . . . 3.2 Packet losses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Constraint violation . . . . . . . . . . . . . . . . . . . . . . 3.2.2 Error on the expected cost function . . . . . . . . . . . . . . Packet losses in OMAS . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Discrete-event modelling . . . . . . . . . . . . . . . . . . . . 3.3.2 Constraint violation . . . . . . . . . . . . . . . . . . . . . . 3.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Replacement of an agent in OMAS, while preserving the topology. . 1.2 Time-invariant finite superset of agents characterized by different active agents. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Multi-mode multi-dimensional systems as a tool to study OMAS. . 2.1 Random interactions with different numbers of active agents n a in a network with n = 18 possible agents. . . . . . . . . . . . . . . . . 2.2 Computation of J noise and the relative error for a star graph S Computation of J noise and the relative error for a line graph P n with growing n. In the left plot, the solid blue line corresponds to the lower bound (2.17) and the dashed red line is the upper bound (2.18). 2.4 Computation of J noise for 2-D and 3-D grids with growing n. . . . 2.5 Computation of J noise and the relative error for a complete graph K Computation of J noise and the relative error for a sequence of Erdős-Rényi graphs with growing n. . . . . . . . . . . . . . . . . . . . . . 2.7 Computation of J noise

	x	CONTENTS
	Contents 3.3 List of Figures
	Abstract 1.1	i
	Aperçu de la thèse	iii
	1 Introduction
	1.2 Randomized consensus with activation/deactivation of agents sub-
	ject to noise disturbance
	vii 2.1 ix

1 Classical multi-agent systems . . . . . . . . . . . . . . . . . . . . . 1.2 Open multi-agent systems (OMAS) . . . . . . . . . . . . . . . . . . 1.2.1 The case of replacements . . . . . . . . . . . . . . . . . . . . 1.2.2 Main challenges in OMAS . . . . . . . . . . . . . . . . . . . 1.2.3 Potential approaches to OMAS: pros and cons . . . . . . . . 1.3 n with growing n. In the left plot, the solid blue line corresponds to the lower bound (2.15) and the dashed red line is the upper bound (2.16). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 n with growing n. In the left plot, the solid blue line corresponds to the lower bound (2.19) and the dashed red line is the upper bound (2.20). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6

  .36) Similarly to the closed system we analyze the constraint violation metric J constr defined in(3.18).

	Proposition 3.2 Consider the weighted gradient descent given by (3.10), (3.11)
	and (3.36). Under Assumption 3.5, for any positive scalar h ≤ constraint violation metric satisfies:	2 β(pλn(L)+2-2p) , the

A neighbor of an agent i in a multi-agent system is an agent that is directly connected to i.

R. Vizuete and C. Monnoyer de Galland equally contributed to this paper and can be considered as the first author.

Graph G P is formed by adding edges between two nodes i, j in V when Pij > 0, for all i, j, including possible self-loops.

We define o(g(n)) as the set[START_REF] Cormen | Introduction to Algorithms[END_REF] o(g(n)) = {f (n) : for any positive constant c > 0, there exists a constantn 0 > 0 such that 0 ≤ f (n) < cg(n) for all n ≥ n 0 }.

The sequence starts with a single node. At the n-th iteration, a new node is added and it is connected to node i with probability 1 -i/n. Then, every pair of nonadjacent nodes is connected with probability 2/n.

This assumption about the preservation of the value of x i by the new agent will be relaxed in Section 4.5 corresponding to the second part of this chapter where tools inspired by online optimization are used for the analysis.

The state x i (t) in a SIS epidemic can also be interpreted as the fraction of a population infected at time t. However, under this interpretation, the replacement of x i may not represent the real behavior of a epidemic since a change of the state could imply that a considerable fraction of the population became healthy abruptly.

Acknowledgments

Acknowledgments vii

Remark 5.1 (Other systems) The proof of Proposition 5.2 did not use the equation (5.1), only the definition of V (x(t)). Therefore, this proof remains valid for other systems.

Along this work, we make the following assumption about the occurrence of replacements in the evolution of the epidemic. Assumption 5.2 (Replacement process) The replacements instants are determined by a homogeneous Poisson process with rate λ r > 0. During a replacement, an agent j ∈ {1, . . . , n} is chosen uniformly and it is assigned a new state x r determined by a random variable Θ, which takes values according to a continuous distribution with support in the interval [0, 1], with mean m and variance σ 2 .

Pure replacements ẋ(t) = 0

To analyze the behavior of the aggregate function during the jumps, first, we consider the setting where V (x(t)) does not change between the jumps. Notice that the jumps of the state x(t) corresponds to the jump of the aggregate function V (x(t)). In this case, the value of V (x(t)) is given by:

where ∆V k is the size of the jump of the function V (x(t)) at the T k jump time of the Poisson process N t and

The aggregate function satisfies

for almost all t, except in a countable number of time jumps. Figure 5.3 presents a sample path of V (x(t))

where the jumps are determined by the Poisson process N t . Then, we have that the aggregate function can be expressed as [START_REF] Privault | Stochastic Finance: An Introduction with Market Examples[END_REF]:

Appendix A

Mathematical background

A.1 Graphs

A graph can be defined as a pair G = (V, E) where V = {1, . . . , n} is a set of vertices or nodes and E ⊆ {(i, j) ∈ V × V : i ̸ = j} is the set of edges. An undirected graph has edges with no direction and a simple graph does not contain self-loops or multiedges. A graph is connected if each node is reachable from any other node via a path (a sequence of adjacent edges).

The adjacency matrix of a graph A = [a ij ] ∈ R n×n is defined by a ij > 0 if (i, j) ∈ E and a ij = 0 otherwise. This matrix is real symmetric non-negative with real eigenvalues [START_REF] Van Mieghem | Graph Spectra for Complex Networks[END_REF]. If the graph has weights w ij in the edges, we design the entries of the adjacency matrix as a ij = w ij . If we do not consider weights in the edges, we define the entries of the adjacency matrix as a ij = 1 for all the edges. We consider the order of the eigenvalues of the adjacency matrix as:

Another matrix associated with a graph is the Laplacian matrix defined as

, and all the eigenvalues λ i (L) have nonnegative real part. Furthermore, 0 is always an eigenvalue of the Laplacian matrix. We consider the order of the eigenvalues of the Laplacian matrix as:

Specifically, this matrix is symmetric and positive semidefinite if the graph is undirected. If the graph is connected, we have λ 2 (L) > 0. The Laplacian matrix can be expressed as:

where D = diag[d 1 , . . . , d n ] and d i is the degree of the node i (i.e., d i is the ith row-sum of A).

125 non-negative integer K, such that for any k, l, any set I kl = I k × I l and pairs (x 1 , y 1 ) and (x 2 , y 2 ) ∈ I kl we have that:

Definition A.3 (Large enough n [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF]) Given a piecewise Lipschitz graphon W (as per Definition A.2) and ν < e -1 , n is large enough if n satisfies the following conditions:

(α k -α k-1 ), (A.2a)

Lemma A.1 (Theorem 1 [START_REF] Avella-Medina | Centrality measures for graphons: Accounting for uncertainty in networks[END_REF]) Let W be a piecewise Lipschitz graphon (as per Definition A.2) and G a graph with n nodes sampled from W . Then for n large enough with probability 1:

and with probability at least 1 -ν:

By considering a constant value of ν, the difference between the graphons in the operator norm is bounded by: |||T

Lemma A.2 Let W be a piecewise Lipschitz graphon and G a graph with n nodes sampled from W . Then for n large enough, with probability at least 1 -ν:

Proof: First, we have

It is easy to see that W G is a step graphon that takes values in {0, 1} and we can apply the inequality derived in [START_REF] Janson | Graphons, cut norm and distance, couplings and rearrangements[END_REF]Remark 10.8], so that:

Applying the relation ∥W -

L 1 yields the desired result. ■