A membrane computing approach to the generalized Nash equilibrium

A membrane computing approach to the generalized Nash equilibrium Alejandro Luque-Cerpa 1 •Miguel A ´. Gutie ´rrez-Naranjo 2 Accepted: 17 March 2025 The Author(s) 2025 Abstract Generalized Nash Equilibrium is an extended version of the standard Nash Equilibrium with important implications in reallife problems such as economics, wireless communication, the electricity market, or engineering among other areas. In this paper, we propose a first approach to computing Generalized Nash Equilibria using Membrane Computing techniques. We model an efficient P system that, based on Euler’s method, computes approximations of Generalized Nash Equilibria of population games under Brown–von Neumann–Nash dynamics, bridging both areas and opening a door for a flow of problems and solutions in both directions. Keywords Membrane computing Generalized Nash equilibrium Evolutionary game theory 1 Introduction Evolutionary Game Theory (EGT) studies the evolution of a population of agents that interact with each other and get a payoff in each interaction Hofbauer and Sigmund (2000). The obtained payoff depends on the chosen strategies of the agents which participate in the interaction. Each agent selects only one strategy at a time, but this choice can be modified over time. The driving principle in this situation is that individuals tend to be selfish, choosing strategies that result in higher payoffs for themselves. In this context, a Nash equilibrium is reached when no agent can increase its payoff by changing its strategy while other agents maintain their current ones Nash (1951). In a Nash equilibrium problem, all the agents compete among them to maximize their payoffs, and each agent can freely choose its strategy. The generalized Nash equilibrium problem (GNEP) is a variant of the Nash problem introduced in 1952 by G. Debreu Debreu (1952). In a GNEP, the strategy set of each player may also depend on the other players’ strategies. This GNEP models a large number of real-life situations, such as power allocation in a telecommunication system, environmental pollution control, or energy market model (for a detailed survey, see, e.g., Facchinei and Kanzow (2007)). In this paper, we propose to study the GNEP in the framework of Membrane Computing Pa ˘un (2002); Pa ˘un et al. (2010). Membrane Computing is a well-known area of Computer Science that takes inspiration from the biochemical reactions inside the vesicles of living cells. P systems Pa ˘un (2000), the so-called Membrane Computing devices, have been successfully considered to model many dynamic processes in real-life problems Colomer et al. (2011,2010); Garcı ´a-Quismondo et al. (2017). From the initial definition of P systems, many variants have been explored by adding new features to the initial model (see, e.g., Song (2021) for a recent survey). Recently, Probabilistic P systems Cardona et al. (2011), a kind of P system designed to deal with probability distributions in the application of rules, was considered to model the spread of behaviors in structured populations in the framework of EGT Garcı ´a-Victoria et al. (2022). In this paper, we study the GNEP by considering transition P systems with active membranes Pa ˘un (2001), also called membrane polarization. &Alejandro Luque-Cerpa [email protected] Miguel A ´. Gutie ´rrez-Naranjo [email protected] 1 Department of Computer Science and Engineering, Chalmers University of Technology, Gothenburg, Sweden 2 Department of Computer Science and Artificial Intelligence, Universidad de Sevilla, Seville, Spain 123 Natural Computing https://doi.org/10.1007/s11047-025-10014-z(0123456789().,-volV)(0123456789().,-volV) The paper is organized as follows: Sect. 2establishes some background on P systems and specifies the type of P system we use: transition P systems with membrane polarization. Section 3introduces population games under the Brown-von Neumann-Nash (BNN) dynamics, that will be used as the framework to define our P system. In Sect. 4, we describe the design of our P system and analyze its time complexity, showing that it does not depend on the number of players or strategies involved. We also present an experiment to illustrate its functioning. Finally, some conclusions and hints for future work are presented. 2 Transition P systems with membrane polarization In this section, we define the variant of P systems that we used to solve our problem: transition P systems with membrane polarization. Then, in section 2.1, we provide an example of such a P system. After the development of the first model of P system by Gh. Pa ˇun in 1998 Pa ˘un (2000), many variations have been presented. In this work, we use a combination of two proposed variants. The P system designed is a transition P system Pa ˘un (2000) with active membranes Pa ˘un (2001) without division rules, i.e., a transition P system with membrane polarization. A transition P system with membrane polarization of degree q1 is a construct: P¼hC;l;w1;...;wq;ðR1;q1Þ;...;ðRq;qqÞ;iouti where: 1. Cis the alphabet of objects; 2. lis a hierarchical tree-like membrane structure of q membranes that have a polarization among 0;þ;; 3. w1;...;wqare multisets of objects over C; 4. R1;...;Rqare finite sets of evolution rules of the form: •u½va h!u0½v0b hwhere u;u0;v;v0are multisets over C,h2f1;...;qg,his not the label of the root membrane in l, and a;b2f0;þ;g. •½va h!u0½v0b hwhere u;u0;v;v0are multisets over C, h2f1;...;qg,his the label of the root membrane in l, and a;b2f0;þ;g. The difference between the expressions resides in that no objects should be able to enter the skin membrane (the root of l) from the environment. The meaning of these rules can be easily understood as combinations of the following examples: •½ua h!½va h, also expressed as ½u!va h, is an object evolution rule, that transforms the multiset uinto the multiset v. •½ua h!v½ b his a send-out communication rule, that ejects the multiset u, and transforms it into the multiset v. •u½ a h!½vb his a send-in communication rule, that absorbs the multiset u, and transforms it into the multiset v. The general expression considers combinations of these cases, where some multisets can be absorbed into membrane hat the same time as others are transformed or ejected. 5. q1;...;qqare partial order relations over R1;...;Rq, called priority relations. Given two rules r;r0,we represent that rhas higher priority than r0by qr[qr0. Priority indicates what rule should be applied if both are applicable. 6. iout 2f0;1;...;qgis the output region, where 0 represents the environment. A configuration of Pis defined by Ct¼ ððw1;t;a1;tÞ;...;ðw1;t;a1;tÞ;w0;tÞfor an instant t, where wh;t is the multiset of objects in membrane hat instant t,ah;tis the membrane polarization of membrane h, and w0;tis the multiset of objects of the environment. The initial configuration of Pis C0¼ððw1;0Þ;...;ðwq;0Þ;;Þ. We use the notation Ct¼l0to denote specific parts from the configuration Ctwhere only the membranes in the subtree l0from lare considered. For each configuration, the rules are applied in a parallel and maximal way. By maximal, we indicate that no more rules can be applied at the same time. Formally, a multiset Uof rules is maximal if there is no multiset of applicable rules U0such that UU0. If two applicable rules with the same priority are exclusive, this is, triggering one would prevent the other one from triggering, then only one of them is selected at random and applied. As in Garcı ´a-Victoria et al. (2022), the semantics of the P system follow the next principles: (I1) When an object crosses a membrane, its polarization may change. Rules can only be applied if the polarization is appropriate. (I2) If two rules that affect the same membrane can be applied at the same time, and one of the rules changes the polarization of the membrane, both rules are applied. This means that the change of polarization is performed after all other evolution rules are applied. Notice that principle (I2) ensures that rules are applied in a parallel and maximal way. If this principle is not assumed, and a rule can change the membrane polarization, then the order of application of the rules would be important during a single transition step. In that case, multisets of rules that A. Luque-Cerpa, M. Á. Gutiérrez-Naranjo 123 can be applied would not be well-defined, breaking the parallelism and maximality of the system. Notice also that, while the membranes in Phave labels in f1;...;qg, we can always define a set Hof labels with jHj¼qsuch that there is a bijection between the elements of Hand the membranes of l. The same applies to the rules, that we can express as the tuple ðR;qÞ. We use this fact in Sect. 4.1 to provide a better indexing. 2.1 Example of a P system Let P¼hC;l;w1;w2;ðR1;fqr1[qr2gÞ;ðR2;;Þ;iouti be a transition P system with membrane polarization of degree 2 where: •C¼fa;b;cg; •l¼½½ 0 20 1; •w1¼fk0g; •w2¼fa3cg; •R1¼fr1k0½ 0 2!½kþ 2;r2½k0!k0 1g; •R2¼fr3½a2!b0 2;r4½a!c0 2;r5½c!b0 2g; •iout ¼2. We have that U1¼fr1;r3;r4;r5g,U2¼fr3g,U3¼ fr1;r3 4;r5gare multisets of applicable rules. U2is not maximal because U2U1.U1and U3are the only maximal multisets of applicable rules, and they could both be applied in this configuration because there is no priority between the rules involved in each set. r2can not be part of a maximal multiset of applicable rules because rule r1has priority over rule r2(qr1[qr2). The multiset U1would lead to the configuration C1¼½½b2ckþ 20 1, and the multiset U3would lead to C1¼½½bc3kþ 20 1. Because the polarization of membrane 2 changed to þin both configurations, none of the rules in R2can now be applied. Because there are no objects k0in membrane 1, none of the rules in R1can be applied. The computation of the system is then finished for both cases after one transition step. 3 Population games under BNN dynamics The purpose of this section is to introduce population games. Specifically, we introduce population games under BNN dynamics, which are central to this paper. In section 3.1, we give an example of such a game: the Energy Market Game, where players decide when to buy energy and modify their strategies depending on the decisions of the rest until an equilibrium is reached. We use the Energy Market Game as a framework to define our P system, and we explain how the P system can be modified to adapt it to other population games under BNN dynamics. In a population game, we have a society of decisionmaking agents divided into disjoint populations that receive different payoffs depending both on the decisions they make and the decisions the rest of the agents make. The goal of each population is to maximize the payoff received. The decisions that agents can make depend on the population they form part of. Each agent is endowed with a revision protocol, which provides conditional switch rates between strategies according to their associated payoffs Sandholm (2010). These rates allow the agents to change their strategies over time. When the number of agents is large enough, this process can be described by differential equations, referred to as the evolutionary dynamics model (EDM). In EDMs, the agents can be modeled as real numbers, the mass of agents, instead of being modeled as discrete independent entities. There are multiple EDMs, but we focus on a specific EDM known as BNN dynamics Brown and von Neumann (1951), which are described next (see Martinez-Piazuelo et al. (2022) for details). Let us consider a society of agents divided into N2 Z1disjoint populations indexed by P¼f1;2; :::; Ng. Each population k2Pis comprised of a constant mass of decision-making agents mk2R[0. The set of strategies of each agent in population k2Pis SkZ1with 2nk¼jSkj\1. The amount of agents selecting strategy i2Skat population kis denoted as xk i2R0. Notice that agents from different populations k1and k2can select the same strategy iif i2Sk1and i2Sk2. Similarly, the proportion of agents selecting strategy i2Skat population kis denoted as zk i¼xk i=mk. Furthermore, xk¼ðxk i1; :::; xk inkÞ> and zk¼ðzk i1; :::; zk inkÞ>denote the strategic distributions of population k,x¼ðx1>;x2>; :::; xN>Þ>, and z¼ðz1>;z2>; :::; zN>Þ>.Lett2Z0be the discrete-time index; x(t) the value of xat time t;z(t) the value of zat time tand pk iðtÞ2Rthe payoff received by the agents selecting strategy i2Skat population k2P. Following the revision protocol introduced in MartinezPiazuelo et al. (2022), the equations that define the EDM describing the evolution of x(t) over time are: _ xk iðtÞ¼mk½^ pk iðtÞþxk iðtÞX j2Sk½^ pk jðtÞþð1Þ ^ pk jðtÞ¼pk jðtÞ 1 mkX l2Sk xk lðtÞpk lðtÞð2Þ where ½þ¼maxð;0Þ, and _ xdenotes the derivative of x. This EDM is known as the BNN dynamics Brown and von Neumann (1951). A membrane computing approach to the generalized Nash equilibrium 123 Intuitively, Eq. (2) computes the benefit of having more agents following strategy jin population k. The reason is that Equation (2) computes the difference between the payoff obtained by agents xk jand the average payoff obtained at time step t. A positive value of ^ pk jðtÞwould indicate that it would be better for population kto have agents switch to strategy j. Equation (1) can then be used to decide how many agents should switch to other strategies at the next time step, given by the derivatives _ xk iðtÞ. A payoff dynamics model (PDM) that describes the evolution of p(t) is also introduced, defined by: _ lðtÞ¼AxðtÞbð3Þ pðtÞ¼fðxðtÞÞA>lðtÞð4Þ where fis a fitness function that provides the payoff for the strategies chosen at a given population state, lrepresents some constraints over such decisions, and such constraints are given by a matrix Aand a vector bthat depends on the specific problem. In this context, Aand bjust introduce penalizations to the payoff, instead of introducing hard constraints. Since the importance of this system lies in updating the payoff signal p(t) and having a closed-loop configuration between p(t) and x(t), a simplified version of this system, where we remove the constraints over the strategies chosen, is considered: pðtÞ¼fðxðtÞÞ ð5Þ An EDM whose payoff function follows a PDM is then called an EDM-PDM. It is not hard to modify the P system proposed later in Sect. 4.1 to compute the effect of the constraints introduced by Aand bin Eq. (5). Because they are linear, the extra computation time required is constant per iteration t. However, because our goal is to show how can P systems be used for computing Generalized Nash Equilibria (GNE), we limit ourselves to the case without constraints. 3.1 Energy market game In the previous section, we expressed a population game under BNN dynamics through Eqs. (1)(2) and (5). To solve a specific population game, we need to define the payoff function (Eq. 5) considering a specific function fthat is different for each game. Taking this into account, a specific EDM must be selected as a framework to define our P system. Because of this, we consider an example of the Energy Market Game Martinez-Piazuelo et al. (2022)as the framework. In the Energy Market Game, N2Z1players compete to purchase energy over a time horizon of T2Z1time slots. Players who try to purchase energy in the same time slots will end up paying more for the energy, and the base price for energy is higher for some time slots than for others. The goal is to buy energy as cheaply as possible, considering that other players have the same goal. For this problem, we can consider each player a population k, and the agents xk irepresent the amount of energy purchased by player kin the time slot i. Following the notation described at the beginning of Sect. 3, let Ck2RTnkbe a matrix such that each column of Ckhas exactly one element equal to 1 and the rest equal to 0, each row of Ckhas at most one element equal to 1, and the j-th element of the i-th column of Ckis 1 iff player kcompetes in time slot jT. Let C¼½C1;C2;...;CN2 RTnbe the concatenation of the Ckmatrices of all players, where n¼Pk2Pnk. Then Cx corresponds to the collective energy demand for all time slots. Let J:Rn!RTbe the pricing function given by JðxÞ¼DCx þ J, where D2 RTT 0is diagonal and  J2RT 0, and let Qk:Rnk 0!Rbe the individual cost of each player k2P, given by QkðxkÞ¼P i2Skðak i=2Þðxk iÞ2þbk ixk i, where ak i2R0and bk i2R0. Following the results from Martinez-Piazuelo et al. (2022), the payoff function pðtÞ¼fðxðtÞÞ for the Energy Market Game can be expressed by fðxðtÞÞ ¼ SxðtÞ C> JaxðtÞbwhere •M¼diagðm1In1;m2In2;...;mNInNÞ, •S¼diagðC1>DC1;...;CN>DCNÞþR>R, and •R¼½ ffiffiffiffi D pC1;ffiffiffiffi D pC2;...;ffiffiffiffi D pCN. •diag is the operation that constructs a matrix using the input elements as the diagonal, and where the rest of the elements are null. To define the payoff over zk iðtÞ, the following transformation is performed over Eq. (5): pðtÞ¼fðxðtÞÞ ¼ fðMzðtÞÞ ¼SMzðtÞC> JaðMzðtÞÞbð6Þ Equations (1) and (2) also change for zk iðtÞ: _zk iðtÞ¼_ xk iðtÞ mk¼½^ pk iðtÞþzk iðtÞX j2Sk½^ pk jðtÞþð7Þ ^ pk jðtÞ¼pk jðtÞX l2Sk zk lðtÞpk lðtÞð8Þ A. Luque-Cerpa, M. Á. Gutiérrez-Naranjo 123 4 Design and functioning of the P system In this section, we introduce first the general idea behind the design of a P system proposed to compute approximations of GNE for population games under BNN dynamics. Then, we define the P system in Sect. 4.1. After that, we perform a computation analysis in Sect. 4.2, where we indicate the evolution rules defined in Sect. 4.1 that are applied to the configurations. Finally, we include our experimental results in Sect. 4.3. Let us consider the EDM-PDM system introduced in Sect. 3.1 by Eqs. (6), (7), and (8). In this section, a P system that computes approximations of GNE under the BNN dynamics for this system is described. The computation can be summarized in a loop of five stages, represented in Algorithm 1. Stages 1 and 2 are used to compute ^ pðtÞusing Eq. (8), Stages 3 and 4 are used to compute _zðtÞ using Eq. (7), and Stage 5 is used to update the value of z(t). To solve any other EDM-PDM system, only the first stage of the P system has to be modified, while the rest remains unchanged. The fundamental idea behind the system is to compute approximations and discretize the values involved in the EDM-PDM system by rounding to ndecimals and multiplying by 10n. To show the functioning of our P system, we fix n¼2 from now on. However, the system can easily be modified for other values of n, providing better approximations with the cost of a longer runtime. After discretizing, a P system can evolve objects representing z(t)to compute GNE. For n¼2, a single object that represents zk iðtÞ, represents 1% of the agents of population kthat follow the strategy i. For example, if we have 16 objects that represent zk iðtÞ, then zk iðtÞ¼0:16. Formally, to approximate and discretize, we perform roundðx;nÞ¼b10nxc. To obtain the next value of the variables in the next instant tþtstep using the values of instant t, we use Euler’s method Butcher (2016), this is, zðtþtstepÞ¼zðtÞþ_zðtÞtstep. In Algorithm 1, other stop conditions can be easily defined, for example, comparing the z(t) values of one iteration with those of the previous one (in constant time) and stopping if no difference is found, but more rules would be necessary. For the sake of simplicity, our stop condition is to limit the number of iterations in the loop. Because performing multiplications using P systems is not trivial, we define a P system that replicates the Russian peasant multiplication algorithm Cameron (1994)in Appendix A. The reason for using this specific algorithm is that the number of time steps required to compute a multiplication is upper bound by a constant for all multiplications of our P system. We use this multiplication P system as a module for our P system. Algorithm 1 General overview of the P system computation 4.1 Definition of the P system The P system to compute approximations of GNE under the BNN dynamics is defined as the construct: P¼hC;H;l;ðwhÞh2H;ðR;qÞ;iouti where the alphabet of objects is given by: C¼fhProd;k;i;li;hk;i;lijk2P;i2Sk;l¼X j\kjSjjþig [fhProd2;k;i;li;ha;k;i;lijk2P;i2Sk;l¼X j\kjSjjþig [fc;rem;mult0;mult1;prod;prod0;e;e0;pos;q;s0;s1;zneg;zvarp;zvarng [fa;b;d;m;f;y0;p;n;compg[fplj1lX k2PjSkjg [fmi;ki;ai;bi;fi;yij1i6g[fover;p1;err;vg [fy0;0;y0;1;y0;2;y2;0;y2;1;y2;2;y3;0;y4;0;y4;1;y4;2;y4;3;y5;0g [fy7;kjk2Pg[fhAUX;ni;hAUX1;ni;hCLK;nijn0g [fy3;j;i;multzi;0;multzi;1jk2P;i2Sk;1j7g [fCk;hq;ii;di;qi;negi;posijk2P;i2Sk;1j7g [fy0;k;y0;i;y5;j;y5;3;i;y5;11;ijk2P;i2Sk;1j10g [fy5;12;k;wi;compwi;zi;jk2P;i2Sk;1j10g [fEXITi;hEXIT;k;i;l;nijk2P;i2Sk;1j10; n1g[fhINIT;k;i;lijk2P;i2Sk;1j10;n1g; the set of membrane labels is given by H¼f0g[P[fSi;kg8i2Sk8k2P[fRESi;kg8i2Sk8k2P [fMULTi;kg8i2Sk8k2P[fM1;M2g [fMULT2i;kg8i2Sk8k2P[fUPDi;kg8i2Sk8k2P [fACUMkg8k2P; the membrane structure lis represented in Fig. 1, and is defined as follows: •Membrane skin with label 0, inside of which we find: 1. One membrane with label P. 2. Nmembranes with labels P. Inside of each membrane k2P: 2:1. Skmembranes with labels Si;k8i2Sk. Inside each Si;k: A membrane computing approach to the generalized Nash equilibrium 123 – One membrane with label RESi;k 2:2. Skmembranes with labels MULTi;k8i2Sk. Inside each membrane MULTi;k: – One membrane with label M1 – One membrane with label M2 2:3. Skmembranes with labels MULT2i;k8i2Sk. Inside each membrane MULT2i;k: – One membrane with label M10 – One membrane with label M20 2:4. Skmembranes with labels UPDi;k8i2Sk 2:5. One membrane with label ACUMk the output region iout is the skin (label 0); the initial multisets are wP¼y0,wSi;k¼hk;i;lizk i8i2Sk 8k2Pwith zk i:¼b 100 nkcfor 1 i\maxfSkgand zk maxfSkg:¼100 ðjSkj1Þb100 nkc,wRESi;k¼hAUX;0i8i2 Sk8k2P, and for any other membrane m, the initial multiset is wm¼;; and the set of rules Ris given by the following rules, separated by the corresponding stage of Algorithm 1, where krepresents the empty multiset, the rule RSm;nrepresents the n-th rule of the m-th stage, the rules are defined 8k2P,8i2Sk, and l¼P j\kjSjjþi, and qm;nrepresents the priority of the rule RSm;n: Stage 1 (Computate payoff p(t)) To use f(z(t)) in the P system, let j:¼b100 ðC> J bÞcbe the constant part of Eq. (6), and let aj;l:¼bðSMÞjlc and bl:¼bðSMÞll þðaMÞlcfor 1 j;lPk2PjSkjbe the coefficients that will multiply zk iin Eq. 6. Notice that because aj;land blare used to compute products by multiplying them by zk i, and zk iis already multiplied by 100, neither aj;lnor blare multiplied by 100 when rounded. The intuition behind the rules of this stage is that objects hk;i;lirepresent zk iðtÞ, and they will be used to compute plðtÞby generating objects pl(see Eq. (6)). Rules RS1;1,RS1;3, and RS1;4move objects hk;i;li, which represent zk iðtÞ, until they are in membrane P. At the same time, they produce objects that will be used in later stages (c,hProd;k;i;li), and objects that will be used for Fig. 1 Membrane structure of our P system. All membranes start with polarization 0 A. Luque-Cerpa, M. Á. Gutiérrez-Naranjo 123 coordination (Ck). Rule RS1;2generates objects plrepresenting the constant part of Eq. 6. RS1;1½hk;i;li0 Si;k!½c0 Si;khk;i;li RS1;2½y0!pj1 1pj2 2...pjn n0 P RS1;3½hk;i;li0 k!½hProd;k;i;li0 khk;i;liCk RS1;4hk;i;li½ 0 P!½hk;i;li0 P Rules RS1;5and RS1;6are used to coordinate the system. Once the polarization of the membrane Pchanges to þ, rule RS1;7is applied, and the objects plgenerated are mixed with the objects plthat were already in membrane P because of rule RS1;2, computing indeed p(t) as expressed in Eq. 6. RS1;5½C100 1C100 2...C100 N!y10 0 RS1;6y1½ 0 P!½y2þ P RS1;7½hk;i;li!pa1;l 1pa2;l 2...pal1;l l1pbl lpalþ1;l lþ1...pan;l nþ P Rules RS1;8,RS1;9, and RS1;11 are used for coordinating. The coordination is achieved by changing the polarization of membrane P. Rule RS1;10 extracts from Pobjects pl, that now represent plðtÞ¼pk iðtÞ, as objects ha;k;i;li. RS1;8½y2!y3þ P RS1;9½y3þ P!½y4 Prem RS1;10 ½pl P!½  Pha;k;i;li RS1;11 ½y4!y5 P Rule RS1;12 moves each object ha;k;i;lito its corresponding membrane k. These objects will later be used in Stage 2. Rules RS1;13 to RS1;15 are used to coordinate the beginning of Stage 2. Rule RS1;16 is a cleaning rule used to remove objects rem that are now useless. RS1;12 ha;k;i;li½ 0 k!½ha;k;i;li0 k RS1;13 ½y5!y6 P RS1;14 ½y6 P!½ 0 Py7;1y7;2:::y7;N RS1;15 y7;k½ 0 k!½multjSkj 0 k,8k2P. RS1;16 ½rem !ks m,8m2H,8s2fþ;;0g(Cleaning rule). Stage 2 (Compute sums P j2Sk zk jðtÞpk jðtÞ) To compute the products of zk jðtÞand pk jðtÞ, we use the membranes MULTi;k, which work as multiplication modules. These modules compute the product of two numbers and return the result in a maximum of 43 transition steps. Each membrane MULTi;kcontains two membranes, M1 and M2. When objects aare placed in M1, objects bare placed in MULTi;k, an object k1is placed in M1, and the polarization of the three membranes is 0, the multiplication module will compute the product of the numbers represented by objects aand b. Membranes MULTi;kexpel objects d, representing the product, and an object frepresenting that the multiplication process is finished. For the sake of simplicity, we left the details about the multiplication process in Appendix A, including a computation analysis. The intuition behind the rules of this stage is that objects hProd;k;i;liwere produced in Stage 1 to represent zk iðtÞ. These objects are transformed into objects prod. Together with objects ha;k;i;li, that represent pk iðtÞ, and the multiplication modules MULTi;k, we obtain the result of multiplying zk iðtÞby pk iðtÞ. The sum Pj2Skzk jðtÞpk jðtÞis obtained by grouping all the resulting objects in the membranes ACUMi;k. The rest of the rules are for coordination, rounding, or producing objects necessary for later stages. Rules RS2;1to RS2;9are used to coordinate the generation of objects a,b, and k1in membranes M1, MULTi;k, and M1, respectively. Once this is done, the multiplication process will begin, and the product of zk iðtÞ(represented by objects a) and pk iðtÞ(represented by objects b) will be computed. Objects hProd2;k;i;liand posiare also generated for later stages. RS2;1½hProd;k;i;li!prod hProd2;k;i;li k RS2;2½ha;k;i;li!e posi k RS2;3mult0½ 0 MULTi;k!½mult0þ MULTi;k RS2;4prod½ þ MULTi;k!½prodþ MULTi;k RS2;5e½ þ MULTi;k!½eþ MULTi;k RS2;6½mult0þ MULTi;k!½k00 MULTi;krem RS2;7prod !½a0 M1 RS2;8½e!b0 MULTi;k RS2;9k0½ 0 M1!½k10 M1 Rule RS2;10 is used to send objects posito membranes UPDi;k, which will be used in later stages. RS2;10 posi½ 0 UPDi;k!½posi0 UPDi;k Rule RS2;12 takes the multiplication results, represented by objects d, and sends them to membrane ACUMktransformed into objects neg. When all multiplications are finished, the sum Pj2Skzk jðtÞpk jðtÞis given inside ACUMk, represented by objects neg. Then, rule RS2;11 changes the polarization of ACUMkto indicate that the sum is computed. Because zk jðtÞand pk jðtÞwere rounded by multiplying by 100, it is necessary to use rules RS2;13,RS2;14, and RS2;15 to round the product again and eject them from membrane ACUMk. A membrane computing approach to the generalized Nash equilibrium 123 RS2;11 fjSkj½ 0 ACUMk!½y2;0þ ACUMk RS2;12 d½ 0 ACUMk!½neg0 ACUMk RS2;13 ½neg100þ ACUMk!½ þ ACUMkneg RS2;14 ½neg51þ ACUMk!½ þ ACUMkneg, with q2;13 [q2;14. RS2;15 ½neg !kþ ACUMk, with q2;14 [q2;15. Rules RS2;16,RS2;17, and RS2;18 coordinate the beginning of the next stage. RS2;16 ½y2;0!y2;1þ ACUMk RS2;17 ½y2;1þ ACUMk!½ 0 ACUMky2;2 RS2;18 ½y2;2 k!½y3;00 krem Stage 3 (Compute ½^ pk iþand P j2Sk½^ pk jþ) The goal of this rules is, because of Eq. (8), to compute ^ pk iby computing the difference of two numbers: pk iðtÞand Pl2Skzk lðtÞpk lðtÞIn this stage, we compute ½^ pk iþand Pj2Sk½^ pk jþ. From now on, in this stage, let Sk¼fi1;...;ijSkjg. Rules RS3;1and RS3;2are used for coordination. RS3;1½y3;0!y3;1;i1... y3;1;ijSkj0 k RS3;2y3;1;i½ 0 UPDi;k!½remþ UPDi;ky3;2;i Because we need to compute ^ pk iðtÞfor each i, we use rule RS3;3to create jSkjcopies of Pl2Skzk lðtÞpk lðtÞ, each represented by negi. These copies are then sent into membranes UPDi;kby using rule RS3;4. RS3;3½neg !negi1... negijSkj0 k RS3;4negi½ þ UPDi;k!½negiþ UPDi;k Rules RS3;5and RS3;6are used to coordinate the computation of pk iðtÞPl2Skzk lðtÞpk lðtÞby changing the polarization of membranes UPDi;kto -. RS3;5½y3;2;i!y3;3;i0 k RS3;6½y3;3;i½ þ UPDi;k!½y3;4;i UPDi;k0 k Rules RS3;7,RS3;8, and RS3;9compute pk iðtÞPl2Skzk lðtÞpk lðtÞ. If the difference is positive, objects qiremain. Otherwise, because we are interested in using this difference to compute ½^ pk iþ, objects negiare eliminated. ½^ pk iþis then represented by objects qi. RS3;7½negiposi!k UPDi;k RS3;8½negi!k UPDi;k, with q3;7[q3;8. RS3;9½posi!qi UPDi;k, with q3;7[q3;9. Rules RS3;10 and RS3;11 are used for coordination. Rule RS3;12 generates objects qand ejects them from membranes UPDi;k. By doing this, we compute the sum Pj2Sk½^ pk jþ, while ½^ pk iþis still represented by objects qi. Rule RS3;13 coordinates the beginning of the next stage. RS3;10 ½y3;4;i!y3;5;i UPDi;k RS3;11 ½y3;5;i UPDi;k!½y3;6;i0 UPDi;krem RS3;12 ½qi0 UPDi;k!½ 0 UPDi;kqq i RS3;13 ½y3;6;i0 UPDi;k!½ 0 UPDi;ky3;7;i Stage 4 (Compute _zk iðtÞ) Because of Eq. (7), to compute _zk iðtÞwe need to compute first zk iPj2Sk½^ pk jþ. From Stage 3, we have objects qthat represent Pj2Sk½^ pk jþ, and from Stage 2 we have objects hProd2;k;i;lithat represent zk i. As in Stage 2, we can use membranes MULT2i;kas multiplication modules. These modules work the same as MULTi;kfrom Stage 2, computing products in a maximum of 43 transition steps, and with the only difference that MULT2i;kreturns objects di instead of dand objects f1instead of f(see Appendix A for more details). When the multiplication process is finished, we will have objects qifrom Stage 3 that represent ½^ pk iþ, and we will have obtained objects dithat represent zk iPj2Sk½^ pk jþ. We can then compute _zk iðtÞ(see Eq. (7)), resulting in objects zvarp if it is positive, or zvarn if it is negative. From now on, in this stage, let Sk¼fi1;...;ijSkjg. Rules from RS4;1to RS4;12 coordinate the beginning of the multiplication process inside membranes MULT2i;k, multiplying zk i, represented by objects hProd2;k;i;li, and Pj2Sk½^ pk jþ, represented by objects q. RS4;1½y3;7;i1...y3;7;ijSkj!multzi1;0...multzijSkj;00 k RS4;2½multzi;0!multzi;10 k RS4;3½q!hq;i1i...hq;ijSkji0 k RS4;4hProd2;k;i;li½ 0 MULT2i;k!½prod0 MULT2i;k RS4;5hq;ii½ 0 MULT2i;k!½e0 MULT2i;k RS4;6multzi;1½ 0 MULT2i;k!½mult0þ MULT2i;k RS4;7½prod !prod0þ MULT2i;k RS4;8½e!e0þ MULT2i;k RS4;9½mult0þ MULT2i;k!½mult10 MULT2i;krem RS4;10 prod0!½a0 M10 RS4;11 ½e0!b0 MULT2i;k RS4;12 mult1½ 0 M10!½k10 M10 Rules RS4;13 and RS4;14 coordinate the beginning of the computation of ½^ pk iðtÞþzk iðtÞP j2Sk½^ pk jðtÞþin membranes Si;k. RS4;13 fjSkj 1!yjSkj 4;0 A. Luque-Cerpa, M. Á. Gutiérrez-Naranjo 123 RS4;14 y4;0½ 0 Si;k!½y4;1þ Si;k Rules RS4;15 to RS4;23 compute such difference. If the difference is positive (negative), then objects zvarp (zvarn) are produced. RS4;15 qi½ þ Si;k!½s0þ Si;k RS4;16 di½ þ Si;k!½diþ Si;k RS4;17 ½s0!s1þ Si;k RS4;18 ½d100 i!znegþ Si;k RS4;19 ½d51 i!znegþ Si;k, with q4;18 [q4;19. RS4;20 ½di!kþ Si;k, with q4;19 [q4;20. RS4;21 ½s1zneg !kþ Si;k RS4;22 ½s1!zvarpþ Si;k, with q4;21 [q4;22. RS4;23 ½zneg !zvarnþ Si;k, with q4;21 [q4;23. Rules RS4;24,RS4;25, and RS4;26 are just for coordination, preparing the system for the next and final stage. RS4;24 ½y4;1!y4;2þ Si;k RS4;25 ½y4;2!y4;3þ Si;k RS4;26 ½y4;3þ Si;k!½y5;0 Si;krem Stage 5 (Update z(t) and output results) Since Stage 1, we have objects cthat represent zk i. In this stage, we combine them with objects zvarp or zvarn, representing _zðtÞ, to update z(t) using Euler’s method: zðtþtstepÞ¼zðtÞþtstep _zðtÞ. Rule RS5;1is used for coordination. Because we take tstep ¼0:01, rules RS5;2to RS5;10 are used to compute zðtÞþtstep _zðtÞ. Because nothing guarantees that the new values satisfy 0 zk iðtÞ100, there are some special cases to consider. Most of the rules from RS5;11 to RS5;38 are used to deal with these cases, and the rest are for coordination. These cases are further developed in section 4.2, Stage 5. RS5;1½y5;0!y5;1 Si;k RS5;2½zvarn100 c!k Si;k RS5;3½zvarn51 c!k Si;k, with q5;2[q5;3. RS5;4½zvarp100 !p Si;k RS5;5½c!p Si;k,withq5;3[q5;5. RS5;6½zvarn100 !n Si;k, with q5;3[q5;6. RS5;7½zvarn51 !n Si;k, with q5;6[q5;7. RS5;8½zvarn !k Si;k, with q5;7[q5;8. RS5;9½zvarp51 !p Si;k, with q5;4[q5;9. RS5;10 ½zvarp !k Si;k, with q5;9[q5;10. RS5;11 ½y5;1!y5;2comp100 Si;k RS5;12 ½p100 Si;k!½over Si;kw100 i RS5;13 ½y5;2 Si;k!½y5;30 Si;ky5;3;i RS5;14 ½p0 Si;k!½ 0 Si;kp RS5;15 ½over comp100 !k Si;k RS5;16 ½n0 Si;k!½ 0 Si;kn RS5;17 ½comp0 Si;k!½ 0 Si;kcompwi RS5;18 ½p comp !p1 Si;k, with q15 [q18. RS5;19 ½p10 Si;k!½ 0 Si;kwi RS5;20 ½pn!k0 k RS5;21 ½p compwi!wi0 k, with q20 [q21. RS5;22 ½nw i!compwi0 k, with q20 [q22. RS5;23 ½p!err0 k, with q21 [q23. RS5;24 ½n!err0 k, with q22 [q24. RS5;25 ½y5;3!y5;40 Si;k RS5;26 ½y5;4!y5;50 Si;k RS5;27 ½y5;50 Si;k!½ þ Si;ky5;6 RS5;28 ½y5;3;i1...y5;3;ijSkj!y5;40 k RS5;29 ½y5;40 k!½y5;5v100þ krem RS5;30 ½wiv!ziþ k RS5;31 ½wi!kþ k, with q5;30 [q5;31. RS5;32 ½compwi!kþ k RS5;33 ½v!ziþ k, with q5;30 [q5;33. RS5;34 ½y5;5þ k!½ 0 krem RS5;35 zi½ þ Si;k!½ziþ Si;k RS5;36 y5;6½ þ Si;k!½y5;70 Si;k RS5;37 y5;7½ 0 RESi;k!½y5;8þ RESi;k RS5;38 zi½ þ RESi;k!½EXITiþ RESi;k Rules from RS5;39 to RS5;47 are used for coordination and for preparing the output of the P system. Objects hAUX;ni are used to count how many iterations of Algorithm 1 have been completed. Objects hEXIT;k;i;l;nirepresent the values zk iðtÞat iteration n. Objects hINIT;k;i;liare used to reset the P system, preparing Stage 1 for a new iteration of the algorithm. RS5;39 ½hAUX;ni!hCLK;nþ1i100hAUX1;nþ1iþ RESi;k RS5;40 ½y5;8þ RESi;k!½y5;9 RESi;krem RS5;41 ½EXITihCLK;ni RESi;k!½  RESi;khEXIT;k;i;l;ni, with n1 RS5;42 ½hCLK;ni!k RESi;k, with n1 and q5;41 [q5;42. RS5;43 ½y5;9 RESi;k!½ 0 RESi;ky5;10 RS5;44 ½hAUX1;ni!hAUX;ni0 RESi;k, with n1 A membrane computing approach to the generalized Nash equilibrium 123 Hofbauer J, Sigmund K (2000) Evolutionary games and population dynamics. J Amer Stat Assoc. https://doi.org/10.2307/2669431 Martinez-Piazuelo J, Ocampo-Martinez C, Quijano N (2022) Generalized Nash equilibrium seeking in population games under the Brown-von Neumann-Nash dynamics, 2022 European Control Conference 2161–2166. https://doi.org/10.23919/ECC55457. 2022.9838437 Martı ´nez-del-Amor MA, Garcı ´a-Quismondo M, Macı ´as-Ramos LF, Valencia-Cabrera L, Riscos-Nu ´n ˜ez A, Pe ´rez-Jime ´nez MJ (2015) Simulating P systems on GPU devices: a survey. Fundam Informaticae 136(3):269–284. https://doi.org/10.3233/FI-20151157 Nash J (1951) Non-cooperative games. Annals Math 54(2):286–295 Pe ´rez-Hurtado I, Valencia-Cabrera L, Pe ´rez-Jime ´nez M, Colomer M, Riscos-Nu ´n ˜ez A (2010) Mecosim: A general purpose software tool for simulating biological phenomena by means of P systems, IEEE Fifth international conference on bio-inspired computing: theories and applications (BIC-TA) 637–643. https://doi.org/10. 1109/BICTA.2010.5645199 Pa ˘un Gh (2000) Computing with membranes. J Computer Syst Sci 61(1):108–143. https://doi.org/10.1006/jcss.1999.1693 Pa ˘un Gh (2002) Membrane Computing. Springer-Verlag, Germany Pa ˘un Gh (2001) P systems with active membranes: attacking NPcomplete problems, J Automata, LangCombin 6(1): 75-90. https://doi.org/10.25596/jalc-2001-075 Pa ˘un Gh, Rozenberg G, Salomaa A (eds) (2010) The Oxford Handbook of Membrane Computing. Oxford University Press, Oxford, England Sandholm WH (2010) Population games and evolutionary dynamics. MIT Press Song B, Li K, Orellana-Martı ´nD,Pe ´rez-Jime ´nez MJ, Pe ´rez-Hurtado I (2021) A survey of nature-inspired computing: membrane computing. Assoc Comput Mach Comput Surv 54(1):1–31. https://doi.org/10.1145/3431234 Zhang G, Pe ´rez-Jime ´nez MJ, Riscos-Nu ´n ˜ez A, Verlan S, Konur S, Hinze T, Gheorghe M (2021) P systems implementation on GPUs. In: Membrane Computing models: implementations. Springer, Singapore, 163-215. https://doi.org/10.1007/978-98116-1566-56 Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. A. Luque-Cerpa, M. Á. Gutiérrez-Naranjo 123