A priori data-driven robustness guarantees on strategic deviations from generalised Nash equilibria

Automatica 167 (2024) 111746 Contents lists available at ScienceDirect Automatica journal homepage: www.elsevier.com/locate/automatica A priori data-driven robustness guarantees on strategic deviations from generalised Nash equilibria✩ Georgios Pantazis a,Filiberto Fele b,∗,Kostas Margellos c aDelft Center for Systems and Control, TU Delft, 2628 CD Delft, Netherlands bDepartment of Systems Engineering and Automation, University of Seville, Spain cDepartment of Engineering Science, University of Oxford, OX1 3PJ, UK article info Article history: Received 11 April 2023 Received in revised form 24 November 2023 Accepted 15 April 2024 Available online 1 June 2024 Keywords: Generalised equilibrium problem Randomised methods Stochastic game theory Multi-agent systems abstract In this paper we focus on noncooperative games with uncertain constraints coupling the agents’ decisions. We consider a setting where bounded deviations of agents’ decisions from the equilibrium are possible, and uncertain constraints are inferred from data. Building upon recent advances in the so called scenario approach, we propose a randomised algorithm that returns a nominal equilibrium such that a pre-specified bound on the probability of violation for yet unseen constraints is satisfied for an entire region of admissible deviations surrounding it—thus supporting neighbourhoods of equilibria with probabilistic feasibility certificates. For the case in which the game admits a potential function, whose minimum coincides with the social welfare optimum of the population, the proposed algorithmic scheme opens the road to achieve a trade-off between the guaranteed feasibility levels of the region surrounding the nominal equilibrium, and its system-level efficiency. Detailed numerical simulations corroborate our theoretical results. ©2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). 1. Introduction The study of noncooperative games plays a significant role in a panoply of applications ranging from smart-grids (Saad, Han, Poor, & Basar,2012) to communication (Scutari, Facchinei, Pang, & Palomar,2014) and social networks (Acemoglu & Ozdaglar, 2011). In these setups, agents can be modelled as self-interested entities that interact with each other and make decisions based on possibly misaligned individual criteria, while being subject to constraints (local or global) that restrict the domain of their choices. Even though a variety of systems can be analysed by means of deterministic game-theoretic tools (Grammatico, Parise, Colombino and Lygeros,2016;Paccagnan, Gentile, Parise, Kamgarpour, & Lygeros,2019;Scutari et al.,2014), in many applications the decision making procedure is affected by uncertainty. A number of results in the literature have explicitly addressed uncertainty in a noncooperative setting. Specifically, Bopardikar, Borri, Hespanha, Prandini, and Di Benedetto (2013) follow a randomised approach for the special case of stochastic zero-sum ✩The material in this paper was not presented at any conference. This paper was recommended for publication in revised form by Associate Editor Gurdal Arslan under the direction of Editor Florian Dorfler. ∗Corresponding author. E-mail addresses: [email protected] (G. Pantazis), [email protected] (F. Fele), [email protected] (K. Margellos). games. Most results rely on specific assumptions on the probability distribution (Couchman, Kouvaritakis, Cannon, & Prashad, 2005;Singh, Jouini, & Lisser,2016) and/or the geometry of the uncertainty set (Aghassi & Bertsimas,2006;Hayashi, Yamashita, & Fukushima,2005;Nishimura, Hayashi, & Fukushima,2009). To circumvent these limitations, recent developments adopt a data-driven perspective, focusing on the connection of game theory with the so called scenario approach (Calafiore & Campi, 2006). This is based on the idea that an optimisation problem with constraints parametrised by an uncertain parameter — with fixed but possibly unknown support set and probability distribution — can be approximated by drawing samples of that parameter and solving the problem subject to the constraints generated by those samples only; this approximation is known as the scenario program. Standard results in the scenario approach (Calafiore,2010;Campi & Garatti,2008;Campi, Garatti, & Prandini,2009) provide certificates on the probability that a new yet unseen constraint will violate the randomised solution obtained by the scenario program. While the aforementioned results apply to uncertain convex optimisation problems, the works Campi and Garatti (2018) and Campi, Garatti, and Ramponi (2018) paved the way towards the provision of data-driven robustness guarantees to solutions of more general nonconvex problems. In Paccagnan and Campi (2019) and Fele and Margellos (2019,2021), these theoretical advancements were leveraged for the first time in a game-theoretic https://doi.org/10.1016/j.automatica.2024.111746 0005-1098/©2024 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 context, for the formulation of distribution-free probabilistic feasibility guarantees for randomised Nash equilibria. These works provide guarantees for one specific equilibrium point (often assumed to be unique); this was extended in Pantazis, Fele, and Margellos (2020,2022), by providing a posteriori feasibility guarantees for the entire domain. Besides the game-theoretic context, alternative methodologies for set-oriented probabilistic feasibility guarantees have been proposed in the seminal works (Alamo, Tempo, & Camacho,2009;de Farias & Roy,2003), which a priori characterise probabilistic feasibility regions constructed out of sampled constraints using statistical learning theoretic results. More recently, the so called probabilistic scaling (Alamo, Mirasierra, Dabbene, & Lorenzen,2019;Mammarella, Mirasierra, Lorenzen, Alamo, & Dabbene,2022) has been proposed to obtain a posteriori guarantees on the probability that a polytope generated out of samples is a subset of some chance-constrained feasibility region. Following an approach similar to Pantazis et al. (2020), the works (Fabiani, Margellos, & Goulart,2020,2022) deliver tighter guarantees by focusing on variational-inequality (VI) solution sets. The results above follow a standard approach in the gametheoretic literature, where a strict behavioural assumption — the so called rationality — is imposed on the players’ decision making. Namely, the players are viewed as rational agents wishing to maximise their profit (expressed by some given cost function). However, studies have shown that this is unrealistic in practice (Camerer,2003;Kahneman & Tversky,1979;Plott & Smith,2008;Wang, Saad, Mandayam, & Poor,2016) and that agents usually exhibit a boundedly rational behaviour (Rubinstein, 1998), i.e., their decisions can deviate from rationality due to individual biases, behavioural inertia, restricted computational power/time, etc. The consequences of this become relevant in engineering applications, as the human role in technical systems evolves beyond mere users and consumers to active agents, operators, decision-makers and enablers of efficient, resilient and sustainable infrastructures (Lamnabhi-Lagarrigue et al.,2017). To bridge this gap between real-world applications and the cognate literature, here we study games with uncertain constraints, where deviations from a nominal equilibrium are explicitly considered. We follow a randomised approach to approximate the coupling constraints by means of data. In this more general setting, where deviations are considered, providing guarantees for a single solution is devoid of any meaning: indeed, repetition of the game might lead to a different solution in a neighbourhood around the nominal equilibrium, irrespective of the employed dataset. Technically speaking, this renders the identification of the data samples that support the solution (cf. sample compression (Margellos, Prandini, & Lygeros,2015)) a challenging task. Focusing on the class of generalised Nash/Wardrop equilibrium seeking problems (GNEs/GWEs) (Facchinei & Kanzow, 2010), we contribute to the provision of data-driven robustness guarantees for the collection of possible deviations from the equilibrium as follows: (1) Adopting a scenario-theoretic paradigm, we establish a methodology for the provision of a posteriori probabilistic feasibility guarantees for a region around the randomised equilibrium of the game under study. This result (Theorem 1), complements (Fabiani et al.,2020,2022;Pantazis et al.,2020,2022), that instead focus on the entire feasibility region. Focusing on a circumscribed region around a GNE/GWE allows offering tighter probabilistic bounds, while the results of Fabiani et al. (2020,2022) and Pantazis et al. (2020,2022), can be obtained as a limiting case of Theorem 1. (2) We design a data-driven algorithm that returns a GNE/GWE and show that all points in a predefined admissible region surrounding it enjoy a priori probabilistic feasibility guarantees. This result (Theorems 2 and 3), unlike Theorem 1, offers an a priori statement valid for a region that is tunable by the user, modelling possible deviations from a nominal equilibrium that a designer wishes to take into account when incentivising a certain operation profile. A distinctive feature of this result is that it provides a priori guarantees for a set rather than single points (Fele & Margellos,2019,2021;Paccagnan & Campi,2019). These guarantees depend on a prespecified quantity, which in turn can affect the location of the nominal equilibrium and the size of the region for which these probabilistic guarantees hold. As such this region is tunable, unlike Grammatico, Zhang, Margellos, Goulart and Lygeros (2016) where a priori guarantees for a set of solutions are provided, but this set is not controlled by the user and could be arbitrarily narrow. Moreover, the results of Grammatico, Zhang et al. (2016) do not focus on games and follow a fundamentally different approach. Furthermore, when the game under study admits a potential function — whose minimum coincides with some social welfare optimum — our methodology provides a new perspective for trading off the probabilistic feasibility of the region surrounding the nominal equilibrium and its system-level efficiency. (3) We propose an equilibrium seeking algorithm as the machinery to obtain a region surrounding a GNE/GWE over which the aforementioned feasibility guarantees hold. The algorithm relies on a primal–dual scheme and is inspired by seminal developments in Facchinei and Pang (2003). However, the mapping that characterises the algorithm updates differs from those typically encountered in the literature (e.g., see Facchinei and Pang (2003, Ch. 12)). This requires showing that the ad-hoc mapping enjoys certain continuity and co-coercivity properties, thus extending the proof-line of Facchinei and Pang (2003) (see Lemmas 2 & 3, and proof of Theorem 2), a task which is interesting per se. Our contributions compared to the cognate literature are summarised in Table 1. The rest of the paper is organised as follows. In Section 2we provide fundamentals of game theory and the scenario approach. In Section 3.1 we show how the feasibility guarantees for a region around the game solution can be a posteriori quantified. In Section 3.2 we propose a data driven algorithm and prove its convergence to an equilibrium such that the considered neighbourhood of strategic deviations can satisfy prespecified probabilistic feasibility requirements. An illustrative example in Section 4corroborates our theoretical analysis. Section 5concludes the paper and presents future research directions. To streamline the presentation of our results, some proofs are deferred to the Appendix. 2. Preliminaries Notation: All vectors are column unless otherwise indicated. Rn +is the nonnegative orthant in Rn. For an n×nmatrix A, we write A≻0 (A⪰0) when x⊺Ax >0 (x⊺Ax ≥0), for any x∈Rn. We denote by 0q×rthe q×rnull matrix, by Irthe r×r identity matrix, and by 1rthe vector of rones; dimensions can be omitted when clear from the context. eqis the unit vector whose qth element is 1 and all other elements are 0, ∥·∥pthe p-norm operator, and (·)rdenotes the rth component of its vector argument. Bp(x, ρ)= {y∈Rd: ∥y−x∥p< ρ}is the open p-normed ball centred at xwith radius ρ; when pis omitted, any 2 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 Table 1 Classification of results and comparison with cognate literature. Problem class Solution sets Type of feasibility guarantees Tuning References Affine feasibility problems Entire feasible set A posteriori ✗Pantazis et al. (2020, Thm. 6) Convex feasibility problems Entire feasible set A posteriori ✗Pantazis et al. (2022, Thm. 2) Uncertain GNEs GNE solution set A posteriori ✗ Uncertain GNEs/GWEs Subset of feasible deviations around GNE/GWE A posteriori ✗Theorem 1 Uncertain VIs Unique solution A priori/A posteriori ✗Fele and Margellos (2019, Thm. 5), Fele and Margellos (2021, Thm. 8) and Paccagnan and Campi (2019, Cor. 1) Convex feasibility problems (Arbitrary) inner approx. of feasible set A priori ✗Grammatico, Zhang et al. (2016, Thm.2) Uncertain GNEs/GWEs Tunable subset of feasible deviations around GNE/GWE A priori ✓Theorem 3 choice of norm is valid. For a set S,|S|denotes its cardinality, while 2Sdenotes its power set, i.e., the collection of all subsets of S. Finally, given D≻0, projK,D[x] := arg miny∈K(y−x)⊺D(y−x) is the skewed projection of xonto the set K. 2.1. Games with uncertain constraints We consider a population of agents with index set N= {1,...,N}. The decision vector xiof each agent i∈Ntakes value in the set Xi⊆Rn, while x=(xi)N i=1∈X=∏N i=1Xi⊆RnN is the global decision vector that is formed by concatenating the decisions of the entire population. The vector x−i∈Rn(N−1) comprises all agents’ decisions except for those of agent i. In our setup, the cost incurred by agent i∈Nis expressed by a real-valued function Ji(xi,x−i) that depends on local decisions as well as on the decisions from other agents j∈N\ {i}. In the following, with a slight abuse of notation, we can exchange xfor (xi,x−i) to single out agent i’s decision from the global decision vector. Furthermore, we consider uncertain constraints coupling the agents’ decisions. These can be expressed in the form1 Cδ= {x∈X:g(x, δ)≤0}, δ ∈∆,(1) where g:RnN ×∆→Rdepends on some uncertain parameter δtaking values in a support set ∆according to a probability measure P. Feasible collective decisions under this setup can be found by letting every agent i∈Nsolve the following optimisation program, where the decisions x−iof all other agents are given, G:min xi∈Xi Ji(xi,x−i) subject to xi∈⋂ δ∈∆ Ci δ(x−i)⎫ ⎪ ⎬ ⎪ ⎭∀i∈N;(2) here, Ci δ(x−i)= {xi∈Xi:g(xi,x−i, δ)≤0}is the projection of the coupling constraint on Xifor fixed x−iand uncertain realisation δ∈∆. The collection of coupled optimisation programs in (2) for all i∈Nconstitutes an uncertain noncooperative game; we denote it as G. Note that (2) follows a worst-case paradigm, taking into account all possible coupling constraints that can be realised by variations of the uncertain parameter δ∈∆. This can render the solutions of Grather conservative. Furthermore, it is in general not possible to compute a solution for Gwithout an accurate knowledge of, and/or additional assumptions on, the support set ∆and the probability distribution P. To circumvent these limitations, we follow a data-driven paradigm and approximate Gby means of a finite number of samples drawn from ∆, namely 1This formulation can describe deterministic and/or local constraints as special cases. the K-multisample δK=(δ1, . . . , δK)∈∆K. In the sequel, we hold on to the standing assumption that these samples are independent and identically distributed (i.i.d.). Apart from this, no other knowledge on the support set ∆and the probability distribution Pof the uncertain parameter is required. Then, for a given multi-sample δK,(2) can be rewritten as GK:min xi∈Xi Ji(xi,x−i) subject to xi∈ K ⋂ k=1 Ci δk(x−i)⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭∀i∈N.(3) Instead of considering all possible uncertainty realisations δ∈∆ as in (2), we let the data encoded in δKlead agents to their decision by solving (3). We refer to the collection of coupled optimisation programs in (3) as the scenario game GK(the subscript Kimplies dependence on the drawn multi-sample δK). Under standard assumptions, a solution to the scenario game GKexists and the problem is, in contrast to G, tractable using state-of-the-art equilibrium seeking algorithms. 2.2. Variational inequalities and game equilibria Notably — under certain assumptions detailed next — solutions to the game GKcan be retrieved as solutions to a variational inequality (VI), for specific choices of the mapping F:X→ RnN (Facchinei & Kanzow,2010, Thm 3.9): VIK:Find x∗∈ΠKsuch that (x−x∗)⊺F(x∗)≥0 for any x∈ΠK,(4) where ΠK:= X∩⋂K k=1Cδkdenotes the problem domain. A classic game solution concept, which encounters wide application in the literature, is the Nash equilibrium (NE) (Nash,1950). At a NE, no agent can decrease their cost by unilaterally changing their decision. Formally, this can be stated as follows. Definition 1. A point x∗=(x∗ i,x∗ −i)∈ΠKis called a generalised Nash equilibrium (GNE) of GKif, for all i∈N, Ji(x∗ i,x∗ −i)≤Ji(yi,x∗ −i),∀yi∈Xi∩ K ⋂ k=1 Ci δk(x∗ −i). For our analysis, we rely on the following conditions: Assumption 1. For all i∈N,Ji(xi,x−i) is convex and continuously differentiable in xifor any fixed x−i. Assumption 2. (1) For any multi-sample δK∈∆K, the domain ΠKis nonempty. 3 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 (2) The set X=∏N i=1Xiis compact, polytopic and convex. (3) For any δ∈∆,gis an affine function of the form g(x, δ)= a(δ)⊺x−b(δ), where a:∆→RnN and b:∆→R. Note that convexity of the cost function with respect to the agent’s local decision is crucial for the design of tailored algorithms with theoretical convergence guarantees for Nash equilibrium seeking. Under these assumptions, we can determine a GNE as in Definition 1 by solving (4) with F(x)=FNE(x):= ⎡ ⎢ ⎣∇x1J1(x1,x−1) . . . ∇xNJN(xN,x−N) ⎤ ⎥ ⎦.(5) A class of problems of common interest can be modelled by the so called aggregative games (Acemoglu & Jensen,2013;Jensen,2010; Kukushkin,2004), where the cost incurred by agents depends on some aggregate measure — typically the average — of the decision of the entire population. Such a cost can be expressed in (3) by the function Ji(xi, σ(x)), where the aggregate σ:RnN →Rnis defined as the mapping x↦→ 1 N∑N i=1xi. A solution frequently linked to this class of games is the Wardrop equilibrium (WE), a concept akin to the NE but specifically defined in the context of transportation networks (Beckmann, McGuire, & Winsten,1956). The variational WEs of GKcan be expressed by using F(x)= FWE(x):= [∇xiJi(xi,z)|z=σ(x)]i∈N; notice that in this case the second argument of Jiis fixed and set to σ(x), consistently with the notion of WE where agents neglect the impact of their decision on others. We restrict the considered class of variational mappings as follows: Assumption 3. The mapping Fis (1) α-strongly monotone, i.e., (x−y)⊺(F(x)−F(y)) ≥α∥x−y∥2 for any x,y∈X, (2) LF-Lipschitz continuous., i.e., ∥F(x)−F(y)∥ ≤ LF∥x−y∥for any x,y∈X. Assumptions 1 and 3are standard in the game-theoretic literature (Facchinei & Pang,2003;Scutari et al.,2014). Assumption 2 is relatively mild; the affine form of the constraints is exploited in the proposed algorithm (see Section 3) for the convergence to an equilibrium bearing the desired robustness properties. We point out that in general only a subset of solutions to GK can be retrieved through (4): these are referred to as variational equilibria, and enjoy favourable properties over nonvariational ones, as with the former the coupling constraints’ burden is equally split among agents (Harker,1991;Kulkarni & Shanbhag, 2012). The following lemma, adapted from Facchinei and Pang (2003, Thm. 2.3.3), formalises the connection between the solutions to VIKand the GNEs (or GWEs) of GK. Lemma 1. Under Assumptions 1–3, VIKhas a unique solution that is also an equilibrium of GK. For the considered class of VIs, several algorithms from the literature can be employed to obtain a variational equilibrium of GK; see, e.g., Facchinei and Kanzow (2010) and Paccagnan et al. (2019). We remark that, even if not explicitly shown for ease of notation, any solution x∗to GKis itself a function of the drawn multisample δK∈∆K. Probabilistic feasibility guarantees for the unique solution of VIKcan then be provided both in an a priori and a posteriori fashion by resorting to the results in Fele and Margellos (2019,2021) and Paccagnan and Campi (2019). However, these results are tailored to the provision of probabilistic feasibility guarantees for a single point (namely the solution of a VI): any strategic deviation from the equilibrium is not covered by such guarantees. We cover this issue in Section 3. First, we provide some background on the scenario approach. 2.3. Basic concepts in the scenario approach A fundamental notion in the scenario approach is the probability of violation of an uncertain constraint. Definition 2. (1) The probability of violation V:RnN → [0,1]of a point x∈ΠKis defined as the probability that a new yet unseen sample δ∈∆will give rise to a constraint Cδ(as defined in (1)) such that x/∈Cδ, i.e., V(x):= P{δ∈∆:x/∈Cδ}. (2) The probability of violation V:2RnN → [0,1]of a set S⊆ΠKis defined as the worst-case Vamong all the points in S, i.e., V(S)=supx∈SP{δ∈∆:x/∈Cδ}. A data-driven decision-making process can be formally characterised by a mapping — the algorithm — that takes as input the data encoded by the samples and returns a set of decisions. Definition 3. An algorithm is a function A:∆l→2RnN × RnN that takes as input an l-multisample and returns the pair (x∗,S∗ l(x∗)), namely, a point x∗and a solution set S∗ lparametrised by x∗. In the setting we consider, we have x∗∈S∗ l(x∗); however, this ought not to be the case in general. In the following, we interpret the above definition as context-dependent, in that the size lof the input multisample is admitted to vary—all else remaining fixed for a given algorithm A. A key notion, strongly linked to that of algorithm, is the minimal compression set (Margellos et al.,2015). This concept springs from the observation that typically only a subset of the sampled data is relevant to a decision or set of decisions, and all other samples are redundant. Definition 4 (Compression Set).Consider an algorithm Aas in Definition 3. A subset of samples I⊆δKis called a compression for A(δK) if A(I)=A(δK).2As multiple subset of samples can exist that fulfil this property, the ones with the minimal cardinality are called minimal compression sets. If we feed the algorithm with the set of samples corresponding to a compression, then the same decision will be returned as when we feed the algorithm with the entire multi-sample. As established in Margellos et al. (2015), the compression set is related to the notion of support samples (Campi & Garatti,2008) and that of essential constraints (Calafiore,2010). Under certain non-degeneracy assumptions these concepts coincide. 3. Probabilistic feasibility of sets around equilibria 3.1. A first a posteriori result Returning to the scenario game GKin (3), we now consider a more general setup where agents are allowed to deviate from x∗ following, e.g., unmodelled changes in their cost functions; while we suppose that these deviations are feasible with respect to the local constraints, we want to study the feasibility as regards the coupling constraints obtained through sampling. Specifically, the region in which agents’ strategies can deviate from the nominal equilibrium is assumed to lie within a predefined open ball B(x∗, ρ), where ρ > 0 is a fixed radius that denotes the maximum possible distance of agents’ deviations from x∗; the latter 2With some abuse of notation, in the remainder the symbol δKis interpreted as either the i.i.d. sample vector δK∈∆K, or the set comprising its components, i.e., δK= {δ1, . . . , δK} ⊆ ∆, depending on the context. 4 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 Fig. 1. Region S∗ K(in green) obtained as the intersection of the set of deviations B∞(x∗, ρ) around the equilibrium x∗(red dot) with the domain ΠK. The samples producing the constraints in blue are in the compression set of x∗, while those associated with the red constraint are in the compression set of S∗ K; discarding these does not change x∗. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) is assumed to be unique as per Lemma 1. As such, the region of interest is S∗ K=ΠK∩B(x∗, ρ). This is depicted in Fig. 1 using the ∞-norm (any other norm could have been used instead): an algorithm A(see Section 2.3) takes as input a multi-sample δKand returns the region S∗ Karound the solution x∗∈R2of a game with two players whose decisions are defined as scalar quantities. For this pictorial example, ΠK is shaped exclusively by sampled coupling constraints. Any compression set as per Definition 4 for Amust be associated with the solid blue constraints (these form a compression for x∗), and with the dashed red constraint that intersects B(x∗, ρ)—as its removal would change S∗ K. We can quantify the number of samples that form a compression set for the algorithm that returns S∗ Kin an a posteriori fashion as established in Theorem 1. To this end, for a fixed confidence β∈(0,1), let the violation level be defined as a function ϵ: {0,...,K}→[0,1]satisfying (Campi et al.,2018, Eq. (7)) ϵ(K)=1 and K−1 ∑ i=0(K i)(1 −ϵ)K−i=β. (6) Theorem 1. Under Assumptions 1–3, let algorithm Areturn a pair (x∗,S∗ K(x∗)). Fix a confidence parameter β∈(0,1) and a violation level ϵ: {0,...,K} → [0,1]that satisfies (6). We have that PK{δK∈∆K:V(S∗ K)> ϵ(s∗+M)}≤β, where s∗is the number of samples that form a compression set for the equilibrium mapping x∗and M the number of uncertain constraints of ΠKthat intersect S∗ K. Proof. Let (x∗,S∗ K) be the solution returned by Afor some given δK, according to Definition 3. We aim at determining a compression set for A(δK), and use its cardinality to reach the theorem’s conclusion by means of Theorem 2 in Pantazis et al. (2022). This would be the union of: (i) the samples that form a compression set for x∗— i.e., solving the problem using only these would result in the same equilibrium obtained by using all samples —, and (ii) any other sample (not in the compression set of x∗) whose removal can still lead to a change of the region S∗ K. Case (i): Determining a (possibly non-minimal) compression set for x∗can be achieved, as suggested in Campi et al. (2018), by progressively removing samples till a subset that leaves the solution unchanged is determined. We denote its cardinality by s∗. With reference to Fig. 1, this set would be associated with the blue constraints active at x∗. Case (ii): We need to count the samples whose removal does not change x∗but yields a larger region S∗ K(red constraint in Fig. 1). Their number can be upper bounded by the Mfacets of ΠKthat intersect S∗ K. Hence, the number of samples that form a compression set for A(δK) is bounded by s∗+M. Existence of a compression set I with a bound on its cardinality is sufficient for the application of Theorem 2 in Pantazis et al. (2022). The fact that for the minimal compression set |I∗|≤|I| ≤ s∗+Malways holds leads then to the statement of this theorem. ■ It is important to stress that the application of Theorem 1 is agnostic on the choice of the equilibrium seeking algorithm. To use the result of Theorem 1, one needs to quantify (an upper bound of) the number of samples s∗that form a compression set for x∗and (an upper bound of) the number Mof coupling constraints that correspond to facets of S∗ K. While s∗≤nN under Assumptions 1 and 3,3an upper bound for Mcan in general only be achieved a posteriori, i.e., once δKis sampled. In the next section we show how we could obtain a priori bounds for the same quantity. 3.2. A priori probabilistic certificates Consider the scenario game GKand suppose that bounded deviations from the solution are allowed. We model such deviations as a ball of radius ρaround the equilibrium, as in Section 3.1. In contrast to the a posteriori nature of the result therein, our goal here is to achieve an a priori bound. Namely, we aim at establishing the main statement of Theorem 1 with a prespecified violation level, which does not depend on the given multi-sample δK. In other words, we seek a statement — holding with known confidence — of the form V(ΠK∩B(x∗, ρ)) <¯ϵ, with ¯ϵ∈(0,1) a priori fixed. To achieve this we build upon the previous conclusions, which expose a link between the probability of constraint violation and the number Mof facets of ΠK(each originated from some uncertainty sample) that B(x∗, ρ) intersects. In particular, a monotonic relationship follows from (6): the smaller Mthe better, i.e., less conservative, the theoretical feasibility guarantees on constraint violation for the strategies belonging to the feasible region S∗ K surrounding the equilibrium. Also, a smaller value of Mcan result in a larger region for which the guarantees of Theorem 1 hold— due to a smaller portion of B(x∗, ρ) being cut off by intersection with ΠK. This motivates us to study the role of Mas a modulating parameter for the robustness of the feasibility certificates offered for the region S∗ K, as well as the extent of deviation from the nominal equilibrium covered by such certificates. 3.2.1. GNE-seeking algorithm with a priori robustness guarantees We consider an iterative scheme to determine a solution of VIKin (4). In particular, since the problem involves coupling constraints, we build our Algorithm 1upon a primal–dual scheme, where constraint satisfaction is achieved by the use of Lagrange multipliers; similar developments hold for both GNE and GWE problems. To this end, we define the augmented vector y:= (x, µ)∈RnN+mby stacking the global decision vector xand the Lagrange multipliers µ=(µℓ)m ℓ=1∈M⊆Rm +. The set Mdenotes the domain of µ; in the sequel we impose some structure on M once some necessary theoretical ingredients are introduced. As deterministic constraints do not play a role in the evaluation of the robustness guarantees, suppose for ease of exposition that ΠK only comprises uncertain coupling constraints. Let A∈Rm×nN and b∈Rmsuch that ΠK= {x∈X:Ax ≤b},(7a) ∥aℓ∥2=1, ℓ =1,...,m,(7b) 3By arguments similar to those in Pantazis et al. (2022) and Schildbach, Fagiano, and Morari (2012) it can be shown that a tighter bound s∗≤nholds for the game GKin case coupling constraints only concern the aggregate variable. 5 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 Algorithm 1 A priori robust GNE seeking algorithm Require: y(0), ρ ∈R+,δK∈∆K,M≤m, ξ ≥0 1: κ←0 2: repeat 3: y(κ+1) ←projX×M,D[y(κ)−D−1T(y(κ), ρ, M)] 4: κ←κ+1 5: until ∥y(κ+1) −y(κ)∥ ≤ ξ 6: y∗←y(κ+1) 7: return y∗and ΠK∩B(x∗, ρ) where aℓ ⊺denotes the ℓth row of A. Eq. (7) is the irredundant H-representation of the polytopic feasibility region ΠKdefined in (4), where the rows of matrix Aare unit vectors. Property (7b) is key to the second statement in Lemma 2. It entails no loss of generality, since for any A,bforming an equivalent H- representation of ΠK,(7) can be obtained by normalising each row of Aand the corresponding component of bby the row-vector norm. Thus, the pair (A,b) encodes the set of randomised coupling constraints that constitute facets of ΠK.4 The main step of Algorithm 1(step 3) is a projected gradient descent (ascent) update for x(µ) through the mapping T: RnN+m+1×N→RnN+mgiven by T(y, ρ, M):= [F(x)+A⊺µ −(Ax −b+Q(µ, M)ρ)].(8) Tfollows from the primal–dual conditions of the game solution; see Facchinei and Kanzow (2010, Sec. 4.2), Facchinei and Pang (2003, Sec. 1.4.1). Fis the pseudo-gradient mapping defined as in Section 2.2,A,bare as in (7), and ρ:= cρ1m, where cis a constant scaling factor (see Section 3.2.2) and Ma nonnegative integer. In the second block-row of (8), the m−Mleast relevant (based on the multipliers value) coupling constraints are tightened by an amount cρthrough the mapping Q:Rm +×N→ {0,1}m×m. Finally, the asymmetric projection matrix D ≻0 includes the step-size parameter τ > 0 and is defined as D:= [1 τInN 0 −2A1 τIm].(9) Note that the constraint tightening performed in the second block-row of Tis equivalent to preventing B(x∗, ρ) from intersecting these constraints. In other words, Qensures that the number of facets of ΠKintersecting B(x∗, ρ) is at most M, which in turn enables to obtain an a priori estimate of the number of samples that form a compression for S∗ Kand hence on V(ΠK∩ B(x∗, ρ)); this is formalised by Theorems 2 and 3. Since m−M coupling constraints are tightened, smaller values for Mcan result in a more robust and possibly larger region S∗ K; however, they can also move the location of the nominal equilibrium x∗to a somewhat less efficient point towards the interior of ΠK. As we will demonstrate numerically in the sequel, this is the case with potential games (Facchinei, Piccialli, & Sciandrone,2011). 3.2.2. Constraint tightening via mapping Q We define the mapping Qas Q(µ, M):= P⊺(µ)R(M),(10) where 4Formally, A:∆K→Rm×nN and b:∆K→Rmare mappings from the K-multisample to the space of real m×nN matrices and m-dimensional vectors, respectively. •P:Rm→ {0,1}m×mreturns a permutation matrix such that P(µ)µis the vector composed by the elements of µarranged in decreasing order. •R:N→ {0,1}m×mtakes as input the number of coupling constraints M≤mwe allow B(x∗, ρ) to intersect with and returns as output the matrix R(M)=[0m×M 0M×m−M Im−M].(11) Compatibly with the definition of P(·), R(M)P(µ)ρ= (0M ⊺cρ1m−M ⊺)⊺=R(M)ρ, where the last equality holds since all components of ρare equal. As discussed in Section 3.2.1,Q(·,M) allows to tighten the constraints corresponding to the smallest m−Mmultipliers. For this, we use the radius of the sphere that circumscribes B(x∗, ρ). This is ρℓ=cρ∥aℓ∥2=cρ, where the last equality is due to (7b); depending on the choice of norm, c=1 if B(·, ρ) is expressed by a p-norm with p≤2, c=√notherwise. Conversely, at most Mconstraints can intersect B(x∗, ρ) upon convergence of the algorithm. Let L(M)⊆ {1,...,m}contain the indices of the Mlargest multipliers. Then, ℓ∈L(M)⇔(Q(µ, M)ρ)ℓ=0, and the second block row of Tin (8) expresses {aℓ ⊺x≤bℓif (Q(µ, M)ρ)ℓ=0, aℓ ⊺x≤bℓ−cρif (Q(µ, M)ρ)ℓ=cρ. (12) Illustrative example: Let ΠKresult from the intersection of 3 hyperplanes and allow B(·, ρ) to intersect at most M=1 of them. From (11),R(M)=[0 0 0 0 1 0 0 0 1 ]. At iteration κof Algorithm 1, let the multiplier vector µ(κ)=(µ(κ) ℓ)3 ℓ=1be such that µ(κ) 2> µ(κ) 1> µ(κ) 3. Then, P(µ(κ))=[0 1 0 1 0 0 0 0 1 ]is the permutation matrix such that P(µ(κ))µ(κ)=(µ(κ) 2µ(κ) 1µ(κ) 3)⊺. So Q(µ(κ),M)ρ= P⊺(µ(κ))R(M)ρ=(cρ0cρ)⊺, where P⊺(·) applies the correct ordering to the vector R(M)ρ. Suppose µ(j) 2> µ(j) 1> µ(j) 3holds for all j≥κ. Then, at convergence, it follows from (12) that B(x∗, ρ) will not be intersecting the constraints associated to µ1and µ3, whereas it could be intersecting the hyperplane associated to µ2. 3.3 Convergence analysis and main result Due to Q, mapping Tis discontinuous on X×Rm. To circumvent this, we restrict the multipliers to the set Mon which we impose some structure granting continuity of Ton X×M. To this end, let Z:= [ζ, +∞)∪ {0}, for some small ζ > 0, i.e., Z⊂R contains all nonnegative scalars which take value greater than ζ when nonzero. Assumption 4. Let Λbe an arbitrarily large compact set. M admits the form M:= {µ∈Λ:(P(µ)µ)ℓ+1≤(P(µ)µ)ℓ−ζ, ∀ℓ=1,...,m−1}∩Zm.(13) Recalling that P(µ)µrearranges the multipliers in descending order, the set Mcontains all vectors where the difference between every pair of strictly positive components—and the distance of the smallest of these from zero — is no less than ζ. We note that (13) is the union of q=m! + m+1 disjoint convex subsets of Rm +, each of which we denote as Mj, i.e., M= ⋃q j=1Mj;Fig. 2 illustrates this set for m=3. It is therefore possible to compute the projection in line 3 of Algorithm 1by, e.g., projecting on Mj, for j=1,...,q, and then setting y(κ+1) to be the solution among these that results in the minimum distance 6 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 Fig. 2. Domain Mof the Lagrange multipliers associated to the coupling constraints, for ζ=0.2 and m=3. This results in q=10 convex subsets, including the origin and a portion of the axes. from y(κ)−D−1T(y(κ), ρ, M). Still, the projection on Mcan be computationally intensive if qis large. Imposing on Mthe structure of (13) endows Twith the desired nonexpansiveness properties that are exploited in the proof of Lemma 3. In the numerical implementation of the algorithm, ensuring µ∈Mcan possibly introduce small perturbations in the multipliers — compared to standard formulations where µ∈Rm + — which in turn could produce a slight violation of the constraints (this can be controlled through the magnitude of ζ). We note that Mis compact by construction due to the intersection with the compact set Λin (13) which can, however, be arbitrarily large thus not impacting the result numerically. Compactness is used in the proof of Theorem 2;Remark 1 discusses cases where this requirement can be lifted. Lemma 2. Define T as in (8)–(11), where A,b satisfy (7). Then, for any µ, µ′∈M,µ= µ′, there exists an integer 0≤h≤M such that (µ−µ′)⊺(Q(µ, M)−Q(µ′,M))ρ≤ −hζcρ. (14) Lemma 3. Consider T as in (8)–(11), where A,b satisfy (7) and M=⋃q j=1Mjas in (13). For each j =1,...,q, let VI(X× Mj,T)denote the VI problem defined by the map T restricted to the subdomain X ×Mj. Under Assumptions 1–3 the following holds: (1) T is continuous on X ×M. (2) Let D as in (9) and set τ > 0such that τ < min{−L2 F+√L4 F+4α2∥A∥2 2α∥A∥2, −ρ(1+∥A∥2)+√ρ2(1+∥A∥2)2+16ζ2∥A∥2 4ζ∥A∥2}.(15) Then, for any j =1,...,q, Algorithm 1converges to a solution of VI(X×Mj,T), when the gradient step in line 3 is projected on the corresponding subdomain, for any y(0) ∈ X×Mj. Continuity of the mapping is essential for the theoretical convergence of Algorithm 1. The second part of Lemma 3 provides an admissible range of values for τsuch that Algorithm 1converges to a solution of VI(X×Mj,T) if at each iteration the projection in line 3is performed on the (convex) subdomain Mj⊂M, j∈ {1,...,q}. The stepsize τis chosen such that conditions standard in NE seeking are satisfied and oscillations among multiple equilibria are avoided. Still, we are interested in establishing convergence on the entire domain M, so at each iteration the projected solution might belong to a different subdomain. This does not trivially follow from the second part of Lemma 3; therefore, by Lemmas 2 and 3we establish an additional condition on τsuch that Algorithm 1retrieves a solution of VI(X×M,T). Theorem 2. Consider Assumptions 1–4. Fixed 0≤M≤m, assume the domain ΠKis nonempty for any of the (m M)combinations of constraints tightened as in (12). Let D ≻0be defined as in (9), where τsatisfies (15) and τ < −(¯ C+¯ R)+√(¯ C+¯ R)2+2ζ¯ R 2¯ R,(16) where ¯ R=max{supx∈Xsupµ∈M∥2A(F(x)+A⊺µ)∥,supx∈X∥Ax − b∥}and ¯ C=cρ√m−M. Then Algorithm 1converges to a solution of VI(X×M)for any initial condition y(0) ∈X×M. Note that as µ(κ)→µ∗, we have Q(µ(κ))→Q(µ∗)=: Q∗. Then, the solution returned by Algorithm 1is the equilibrium of a variant of GKwith m−Mtightened constraints (follows from (12) with Q(µ) replaced by Q∗). Remark 1 (Relaxing Compactness).Theorem 2 still holds when Λ=Rmin the definition of Min (13) if for all multi-samples, (i)Ais full row-rank, or (ii) all elements of Aare positive. (i) To show this, consider mapping Tand matrix Din (9). The multipliers’ update involves projecting (weighted according to D) on M, the term µ(k)−2τ2A(F(x(k))+A⊺µ(k)) +τ(Ax(k)−b+Q(µ(k),M)ρ).(17) Since Xis compact, there exists a subsequence {κi}i∈Nsuch that limi→∞ κi= ∞, limi→∞ x(κi)=x, for some x∈X. It suffices to show that the sequence of multipliers {µ(κi)}∞ i=1remains bounded (all arguments in the proof of Theorem 2 from (C.1) onwards remain unaltered). For the sake of contradiction, assume that there exists at least one element of µ(κi)that tends to infinity across the considered subsequence. Let then µ(κi)=(µ(κi) ∞µ(κi) F), where based on our contradiction hypothesis limi→∞ ∥µ(κi) ∞∥ = ∞ while limi→∞ ∥µ(κi) F∥<∞. (Taking the first elements of µ(κi)to be the ones that tend to infinity is only to simplify notation and is without loss of generality.) Let then A= [A∞⊺AF ⊺]⊺,b=(b∞bF) be the corresponding partition of Aand b, respectively, where A∞,b∞are non-empty by hypothesis. To have ∥µ(κi) ∞∥ → ∞, we need the terms that are integrated in the multipliers’ update, i.e., last two terms in (17), to be positive for all i(in fact across a subsequence), which since τ > 0 is equivalent to A∞A⊺µ(κi)<−A∞F(x(κi)) +1 2τ(A∞x(κi)−b∞+(Q(µ(κi),M)ρ)∞),(18) where (·)∞denotes the elements of its argument corresponding to µ(κi) ∞. Notice that A∞A⊺µ(κi)=A∞A∞⊺µ(κi) ∞+A∞AF ⊺µ(κi) F. As such, we have A∞A∞⊺µ(κi) ∞<−A∞AF ⊺µ(κi) F−A∞F(x(κi)) +1 2τ(A∞x(κi)−b∞+(Q(µ(κi),M)ρ)∞).(19) However, limi→∞ x(κi)=x∈X, and (Q(µ(κi),M)ρ)∞≤cρfor all i, while by Lemma 3,Fis continuous over the domain of multipliers 7 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 satisfying (13). Moreover, µ(κi) Fcontains the components of µ(κi) that remain finite. Therefore, the limit as i→ ∞ of the righthand side of (19) is finite. Due to the assumed full row-rank structure of A, matrix A∞A∞⊺is invertible, hence (19) implies limi→∞ ∥µ(κi) ∞∥<∞, establishing a contradiction showing that the subsequence {µ(κi)}remains bounded. (ii) If all elements of Aare positive, and since aℓ ⊺aℓ=1, for all ℓ=1,...,m, all arguments of case (i) remain the same with the only difference that we directly have that ∥A∞A∞⊺µ(κi) ∞∥ ≥ ∥µ(κi) ∞∥. The next result accompanies the region S∗ K=ΠK∩B(x∗, ρ) of strategic deviations from the equilibrium x∗with a priori probabilistic feasibility guarantees that can be tuned by means of M. It should be noted that Theorem 2 establishes that there exists a choice of τto guarantee convergence of Algorithm 1. The admissible range of values for τis explicit via (15),(16), but difficult to quantify due to ¯ R. Numerical evidence suggests that selecting a small enough value is sufficient for convergence. Theorem 3. Consider Assumptions 1–4. Let x∗and S∗ K=ΠK∩ B(x∗, ρ)be returned by Algorithm 1; fix ϵ∈(0,1) and M. We then have that PK{δK∈∆K:V(S∗ K)≤ϵ} ≥1− nN+M−1 ∑ i=0(K i)ϵi(1 −ϵ)K−i.(20) By Definition 2, Theorem 3 guarantees that for any point in S∗ K, the probability of constraint violation is bounded by ¯ϵ, with confidence at least 1 −∑nN+M−1 i=0(K i)ϵi(1 −ϵ)K−i. The dependence of this term on Mgives us an additional degree of freedom in trading the robustness of the solution for its associated probabilistic confidence. The choice of Mcan also have an effect on the size of S∗ K, as well as on the location of x∗, thus resulting in a trade-off between performance and robustness. For the case in which the coupling constraints concern exclusively the aggregate variable, it can be shown that the upper limit of the summation in the right-hand side of (20) can be replaced by n+M−1, as nis the dimension of the aggregate vector. This allows to state (20) with a much higher confidence of 1 −∑n+M−1 i=0(K i)ϵi(1 −ϵ)K−i; for details, we refer the reader to Pantazis et al. (2022), where the notion of support rank is exploited (Schildbach et al.,2012). 4. Numerical example Consider a game with Nagents whose decisions are subject to deterministic local constraints and uncertain coupling constraints on the aggregate decision: min xi∈Xi xi ⊺(Cσ(x)+d) subject to bδk≤σ(x)≤bδk, k=1,...,K ⎫ ⎪ ⎪ ⎬ ⎪ ⎪ ⎭∀i∈N,(21) where C≻αIn, for some α > 0, and d∈Rn. Note that a structure similar to our numerical example has been considered in applications of aggregative games such as electric vehicle charging and traffic management under uncertainty (Fele & Margellos,2019, 2021;Paccagnan et al.,2019). We impose no knowledge of ∆ and P; we rely instead on a scenario-based approximation of the game, whereby each sample δk∈δKgives rise to bδk,bδk. Eq. (21) is an aggregative game in the form of (3). In this instance, we assume each agent’s action has negligible effect on the aggregate, and accordingly consider a GWE-seeking problem. Following the definition of FWE (Section 2.2), we get F(x)=FWE = [Cσ(x)+c]i∈N. It can be verified that Fis Lipschitz continuous and strongly monotone with respect to σ:byFacchinei and Pang (2003, Thm. 2.3.3), (21) admits a unique aggregate equilibrium σ∗=σ(x∗).5 We employ Algorithm 1to seek a WE x∗such that, by fixing M, a prespecified theoretical violation level is guaranteed for the set ΠK∩B(x∗, ρ). Due to uniqueness of σ∗, all sets B(·, ρ) — parametrised by any x∗solving (21) — are projected on the unique ball B(σ∗, ρ/N) in the aggregate space. Also note that by definition of σ, at most nnon-redundant samples will contribute to define the domain Πσ K:= 1 N(1N ⊺⊗In)X∩(⋂K k=1Cδk)in (21). For the derivation of the robustness guarantees, we can thus restrict our attention to S∗ K=Πσ K∩B(σ∗, ρ/N)⊆Rn. As remarked at the end of Section 3.3, we can apply (20) with the upper limit in the summation involved replaced by n+M−1. For the case n= 2,N=50, and different choices of M,Fig. 3 depicts the projected iterations {σ(x(κ))},κ=1,2, . . . generated by Algorithm 1on the space Πσ K. It can be observed how the region S∗ Kchanges as the value of Mis modified. It is worth noting that in this case F(x) is integrable—this can be inferred by Facchinei and Pang (2003, Thm. 1.3.1) since the Jacobian of the game is symmetric, i.e., ∇xF(x)= ∇xF(x)⊺. Therefore, a GWE x∗can also be obtained by solving min x∈Xσ(x)⊺Cσ(x)+d⊺σ(x) subject to bδk≤σ(x)≤bδk,k=1,...,K. (22) In other words, this game admits a potential function E(x):= σ(x)⊺Cσ(x)+d⊺σ(x), whose minimisers correspond to GWEs. E can be interpreted as the total cost incurred by the population of agents, and its minimisation leads to the optimum social welfare. The contour lines of Eare depicted in Fig. 3: since x∗minimises E(·), σ∗lies on the contour associated to the minimum value of Ewithin the feasible domain. Lower values of Mresult in larger regions for which guarantees are provided. Fig. 4 shows how the sequence {E(x(κ))}κ=1,2,... converges to the minimum potential within the possibly tightened feasibility region. It can be observed how in this case the efficiency of the equilibrium decreases as smaller values of Mare chosen. The three panels in Fig. 4 show the trade-off between system level efficiency and the guaranteed robustness levels. The lower the value of M, the lower the empirical constraints violation—corresponding to a better confidence bound in the right-hand side of (20). 5. Concluding remarks This work proposes a data-driven equilibrium-seeking algorithm such that probabilistic feasibility guarantees are provided for a region surrounding a game equilibrium. These guarantees are a priori and the region that is accompanied with such a probabilistic certificate is tunable. For games that admit a potential function, the proposed scheme is shown to achieve a trade-off between cost and the level of probabilistic feasibility guarantees. In fact, our scheme returns the most efficient equilibrium such that the predefined guarantees are achieved. Proving this conjecture is left for future work. Moreover, current work investigates a distributed implementation of the proposed equilibrium seeking algorithm. 5We note that this case slightly transcends the conditions in Theorem 2, as F does not comply with Assumption 3–(1). Convergence of Algorithm 1(following from the nonexpansiveness of Ton each subdomain Mj) can still be ensured here due to the affine structure of F; cf. Facchinei and Pang (2003, Sec. 12.5.1). 8 G. Pantazis, F. Fele and K. Margellos Automatica 167 (2024) 111746 Fig. 3. Iterates generated by Algorithm 1(blue diamonds) for different choices of M. In this numerical instance, N=50, ρ=10, and Xi:= {xi∈Rn:xi∈ [xi,xi]}, with xi=(0 0), xi=(3.5 3.5). The randomly generated coupling constraints form the rectangular feasibility region Πσ K(delineated by the solid black line). The red-shaded region represents the intersection between the latter and the ball B1(σ∗, ρ/N) around the aggregate equilibrium σ∗(red diamond marker). In this instance, its volume increases as larger values for Mare chosen. The value associated to the contour lines of the potential function Edecreases from top-right to bottom-left, showing that σ∗is the unique minimiser in the admissible region (shaded in green) after constraint tightening is performed by the algorithm (see Section 3.2.2). (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) Fig. 4. Potential function E(x(κ)) evaluated along the iterations of Algorithm 1. Lower values of Myield better confidence on the theoretical robustness certificates for the considered region (see Theorem 3), which results in a lower empirical probability of constraint violation. On the other hand, the system-level efficiency of the equilibrium increases for higher values of M. Acknowledgements F. Fele gratefully acknowledges support from grant RYC2021- 033960-I funded by MICIU/AEI/10.13039/501100011033 and European Union NextGenerationEU/PRTR, as well as from grant PID2022-142946NA-I00 funded by MICIU/AEI/10.13039/ 501100011033 and by ERDF/EU. The work of K. Margellos was partially supported by EPSRC under Grant EPP03277X/1. Appendix A. Proof of Lemma 2 Let µ, zbe arbitrary vectors in Mand, as in the proof of Lemma 3, define µ,  zas the vectors composed by rearranging the elements of µ, zin decreasing order. According to this arrangement, let Iµ= {i1,i2,...,im}be the ordered set of indices of µ, i.e., ik:µik= µk,k=1,...,m; as a result, i1and imwill be the indices of the largest and smallest components of µ, respectively. Applying a similar definition to z, we denote the corresponding set Iz:= {j1,j2,...,jm}. Then, the first Mindices in Iµand Iz, denoted as Lµand Lz, respectively, are relative to the constraints not tightened by the application of Q(·,M). In other words, for all ℓ∈Lµ, (Q(µ, M)ρ)ℓ=0 — and similarly for z. Vice versa, the complementary sets Lc µ=Iµ\Lµand Lc z=Iz\Lzare such that for all ℓ∈Lc µ, (Q(µ, M)ρ)ℓ=cρ, and for all ℓ∈Lc z, (Q(z,M)ρ)ℓ=cρ. Let q= [Q(µ, M)−Q(z,M)]ρ. We distinguish between the following cases: (1) ℓ∈Lc µ∩Lz: we have (Q(µ, M)ρ)ℓ=cρsince ℓ∈Lc µ, while (Q(z,M)ρ)ℓ=0 as ℓ∈Lz. Then, qℓ=cρ. (2) ℓ∈Lµ∩Lc z: from ℓ∈Lc zwe have (Q(z,M)ρ)ℓ=cρ. On the other hand, since ℓ∈Lµ, (Q(µ, M)ρ)ℓ=0. This results in qℓ= −cρ. (3) ℓ∈(Lµ∩Lz)∪(Lc µ∩Lc z). If ℓ∈Lµ∩Lzthen (Q(·,M)ρ)ℓ=0 for both µand z. Therefore, qℓ=0. Conversely, if ℓ∈ Lc µ∩Lc z, then (Q(·,M)ρ)ℓ=cρfor both µand z, which results again in qℓ=0. The sets Lc µ∩Lz,Lµ∩Lc z, (Lµ∩Lz)∪(Lc µ∩Lc z) are pairwise disjoint and exhaust the set {1,...,m}. Hence we can write U=(µ−z)⊺(Q(µ, M)−Q(z,M))ρ= m ∑ ℓ=1 (µℓ−zℓ)qℓ =∑ i∈Lc µ∩Lz (µi−zi)cρ+∑ j∈Lµ∩Lc z (µj−zj)·(−cρ) =(∑ i∈Lc µ∩Lz µi−∑ j∈Lµ∩Lc z µj    =:U1 +∑ j∈Lµ∩Lc z zj−∑ i∈Lc µ∩Lz zi    =:U2 )cρ, (A.1) Now, notice that for any i∈Lc µ∩Lz⊆Lc µand j∈Lµ∩Lc z⊆Lµ, we have by definition of Lµand Lc µthat µi≤µj(which by (13) only holds with equality if µi=µj=0). With analogous reasoning, we have zi≥zjfor any i∈Lc µ∩Lz⊆Lzand j∈Lµ∩Lc z⊆Lc z. Let h1be the cardinality of the set Lc µ∩Lz, and h2that of Lµ∩Lc z. Then, h1= |Lc µ∩Lz|(a) = |Lz\Lµ|=|Lz|−|Lµ∩Lz| (b) = |Lµ|−|Lz∩Lµ|=|Lµ\Lz|=|Lµ∩Lc z| = h2, where (a) holds since Lµ,Lz⊆ {1,...,M}, and (b) follows from |Lµ|=|Lz| = M. Therefore h1=h2=: hand 0 ≤h≤M, which implies U1≤0 and U2≤0in(A.1). We can observe that U1<0 and U2<0 if Lµ∩Lc zand Lc µ∩Lzare nonempty and the corresponding components of µand zare nonzero. In such a 9