An application of membrane computing to humanitarian relief via generalized Nash equilibrium

Vol.:(0123456789) Journal of Membrane Computing https://doi.org/10.1007/s41965-025-00187-y RESEARCH PAPER An application ofmembrane computing tohumanitarian relief viageneralized Nash equilibrium AlejandroLuque‑Cerpa1· DavidOrellana‑Martín2,3· MiguelÁ.Gutiérrez‑Naranjo3 Received: 1 October 2024 / Accepted: 13 March 2025 © The Author(s) 2025 Abstract Natural and political disasters, including earthquakes, hurricanes, and tsunamis, but also migration and refugees crisis, need quick and coordinated responses in order to support vulnerable populations. In such disasters, nongovernmental organizations compete with each other for financial donations, while people who need assistance suffer a lack of coordination, congestion in terms of logistics, and duplication of services. From a theoretical point of view, this problem can be formalized as a generalized Nash equilibrium (GNE) problem. This is a generalization of the Nash equilibrium problem, where the agents’ strategies are not fixed but depend on the other agents’ strategies. In this paper, we show that membrane computing can model humanitarian relief as a GNE problem. We propose a family of P systems that compute GNE in this context, and we illustrate their capabilities with Hurricane Katrina in 2005 as a case study. Keywords Membrane computing· Game theory· Nash equilibrium 1 Introduction According to [1], the economic cost of damages as a result of global natural disasters in 2023 was 203.35 billion of USA dollars. They include droughts, floods, extreme weather, extreme temperatures, landslides, dry mass movements, wildfires, volcanic activity, and earthquakes. Not only natural but also political disasters, such as migration and refugee crises, need quick and coordinated responses to support vulnerable populations. In such disasters, nongovernmental organizations (NGOs) are the main actors to reduce suffering and mortality and support quality of life. However, several studies (e.g., [2, 3]) point out that humanitarian aid has not been successful due to a lack of coordination, congestion in terms of logistics, and duplication of services. Although there are huge differences between humanitarian logistics and commercial logistics [4], both problems can be formalized in a similar way. In [5], the authors compare fundraising with and without an earmarking option using optimization models. It seems to be the first study where game theory is considered in order to model the interaction between donors and humanitarian organizations. In [3], the authors develop a generalized Nash equilibrium network model for post-disaster humanitarian relief by nongovernmental organizations which is the starting point for our study by using Membrane Computing techniques. In such a paper, the authors consider Hurricane Katrina as a case study, the costliest disaster in the history of the United States. This hurricane caused huge damage to property and infrastructure, left 450,000 people homeless, and took 1833 lives in Florida, Texas, Mississippi, Alabama, and Louisiana. For this work, we take the data of the disaster provided by Nagurney etal. [3]. From a theoretical point of view, this problem can be formalized as a generalized Nash equilibrium (GNE) problem. This is a generalization of the Nash equilibrium problem, where the agents’ strategies are not fixed but they depend on the other agents’ strategies and it can be considered as a problem in the area of Evolutionary Game Theory (EGT) * Alejandro Luque-Cerpa [email protected] David Orellana-Martín [email protected] Miguel Á. Gutiérrez-Naranjo [email protected] 1 Department ofComputer Science andEngineering, Chalmers University ofTechnology, Gothenburg, Sweden 2 Research Group ofNatural Computing, Universidad de Sevilla, Seville, Spain 3 Department ofComputer Science andArtificial Intelligence, Universidad de Sevilla, Seville, Spain A.Luque-Cerpa et al. [6]. Beyond the intrinsic interest of GNE as a theoretical problem, it can be applied to a wide range of applications, such as the control of interacting vehicles [7], the intersection management problem [8] or the energy market [9]. GNE (or its simplified version, Nash equilibrium) has also provided a theoretical model for other helping actions such as blood donations [3, 10], competition for medical supplies [11], or logistics for humanitarian services [12, 13] among many others. Membrane computing is a well-known bio-inspired computing paradigm [14, 15], whose devices, called P systems, have been used to successfully model many reallife problems, such as the dynamics of the population of giant panda in captivity [16], fault propagation paths in power systems [17] or the ecosystem of some scavenger birds [18] among many others. To the best of our knowledge, the first time that a GNE problem was simulated by membrane computing devices was in Luque-Cerpa etal. [19]. This paper follows the research line started in that paper by showing that the Membrane Computing paradigm can be useful in the simulation of humanitarian relief. In opposition to that paper, we propose a family of P systems that computes GNE on a game out of the framework of Evolutionary Game Theory. The paper is organized as follows: Sect.2 recalls some theoretical aspects related to how GNE can be used in the framework of humanitarian relief after a disaster. Section3 gives some basic information about the P system model used for the simulation. In Sect.4, the design of our Membrane Computing device for dealing with humanitarian relief based on the GNE problem is presented. Next, Sect.5 shows, as a case study, the use of our device for simulation of humanitarian relief with the data of Hurricane Katrina, and finally, the paper ends with some conclusions and some future research lines. 2 GNE model forpost‑disaster humanitarian relief The generalized Nash equilibrium problem (GNEP), first presented by G. Debreu [20] in 1952, is a generalization of the Nash Equilibrium problem in which the players’ strategies depend on their rivals’ strategies. Constraints in the game define these dependencies, and the payoffs obtained by each player depend on the strategies selected. Each player can choose only one strategy at a time, although this choice can be changed for subsequent actions. The guiding idea here is that players typically select courses of action that provide them with greater rewards, but the selected strategy may vary over time to adapt to the behavior of the rest of the players. When no player in this situation may change its strategy and whereas other agents keep using their current ones, a Nash equilibrium is reached [21]. Generalized Nash equilibrium problem can be concretized into problems from many different areas. Among other applications, in [3] authors model the humanitarian relief post-disaster problem as a GNEP. This problem is then transformed into an optimization problem. The formulation is as follows. Let us have m NGOs that provide disaster relief to n different locations. NGOs compete to provide relief items to demand points. Let qij ≥0 be the amount of items that the NGO i provides to the demand point j. Let si be the maximum amount of relief items that the NGO i can provide. Then, the following equation must hold for each NGO i: n ∑ j=1 qij ≤s i ; that is, the total of the items provided by each NGO cannot be higher than the maximum relief items that such an NGO can provide. Furthermore, for each NGO i and each location j, we can consider the mapping cij ∶ℝ+ → ℝ+ , which maps qij to cij(qij) ; that is, the mapping cij represents the cost of sending the relief to the location j by the NGO i. Similarly, each NGO i gets utility associated with the relief items it sent to the location j. The utility over all demand points is measured by n ∑ j=1 𝛾ijq ij , where 𝛾ij is a positive factor. Each NGO i has a positive weight 𝜔i associated with such a utility measure that represents the monetization of this utility, and the monetization of an NGO i is measured by 𝜔 i n ∑ j=1 𝛾ijqij . Finally, each NGO i receives funding from donors based on media attention and the visibility of NGOs at location j. These funds are represented by 𝛽 i n ∑ j=1 Pj(q ) , where q is the vector of all the item flows of all the NGOs to all the demand points, Pj(q) represents the funds in donation dollars due to the visibility of all NGOs at location j, and 𝛽i represents the proportion of total donations collected that is received by NGO i. Other constraints are usually imposed by an authority. For example, there are lower ( dj ) and higher bounds ( dj ) for the number of relief items needed at a location j. These constraints are expressed as m ∑ i=1 qij ≥d j and m ∑ i=1 qij ≤d j , which limit the amount of objects distributed by all NGOs to a specific location j. Another important assumption is that NGOs have enough resources to comply with the An application ofmembrane computing tohumanitarian relief viageneralized Nash equilibrium lower bounds of relief items required by all locations, that is, m ∑ i=1 si≥ n ∑ j=1 d j . Finally, the optimization problem faced by the NGO i can be expressed as which, subject to the constraints indicated, is equivalent (Theorem1, [3]) to the following optimization problem: Note (Existence and uniqueness) A solution q∗ to the optimization problem (1) is guaranteed to exist when the objective function consists of continuous functions and the feasible set defined by the constraints is compact. In [3], the functions employed satisfy these conditions, and because the objective function is strictly convex, the uniqueness of the solution is also guaranteed. 2.1 Euler method applied tothegame theory model To solve the minimization problem (1), we can use a variational inequality formulation of the problem1 and then apply the Euler method [22, 23] to obtain closed-form expressions in the qij and Lagrange multipliers 𝜆 i,𝜆 1 j ,𝜆 2 j associated with the constraints. Such Euler method approximates the solution iteratively. We have the following expressions of qt+1 kl at each iteration t for k=1, ..., m ; l=1, ..., n : Minimize −𝛽i n ∑ j=1 Pj(q)−𝜔i n ∑ j=1 𝛾ijqij + n ∑ j=1 cij(qij ) (1) Minimize − n ∑ j=1 Pj(q)− m ∑ i=1 n ∑ j=1 𝜔 i 𝛾 ij 𝛽i qij + m ∑ i=1 n ∑ j=1 1 𝛽i cij(qij ) (2) q t+1 kl =max {0, {qt kl +at( n ∑ l=1(𝜕Pj(q t ) 𝜕qkl )+𝜔k𝛾 kl 𝛽k −1 𝛽 k 𝜕ckl(qt kl) 𝜕q kl −𝜆t k+𝜆1 l t−𝜆2 l t)}} (3) 𝜆 t k=max { 0, 𝜆t k+at ( −sk+ n ∑ l=1 qt kl )} (4) 𝜆 1 l t+1=max { 0, 𝜆1 l t+at ( − m ∑ k=1 qt kl +dl )} (5) 𝜆 2 l t+1=max { 0, 𝜆2 l t+at ( −dl+ m ∑ k=1 qt kl )} where the sequence {at} must satisfy at>0, at → 0 y ∞ ∑ t=0 at =∞ to guarantee the convergence of the iterative scheme [22]. When designing the P system that will model these equations, the main problem to solve is dealing with the two partial derivatives in Eq.2. The discrete encoding of the information of P systems as multisets of objects is an added difficulty in order to deal with derivatives, therefore, we have made the following decisions in the design of the P system: • The functions Pj(q) , which represents the funds in donation dollars due to the visibility of all NGOs at location j, in many cases can be expressed as with kj>0 , thus and This usually gives the term various orders of magnitude lower than the other terms on Eq.2. Because of this, we have decided to ignore it. In fact, in Sect.5, we experimentally show that the impact derived from this decision is low. • Similarly, the cost functions are of the form with akl,bkl >0 , so 𝜕 ckl(q t kl) 𝜕 q kl =2a2 klqkl +2aklb kl . Finally, we can rewrite then Eq.2 as: We note that the conditions that guarantee the existence and uniqueness of the solution have not been compromised. In Sect.4 a P system that simulates the algorithm defined by Eqs.3, 4, 5 and 6 is designed. P j(q)=kj √ √ √ √ m ∑ i=1 qij 𝜕 Pj ( q t) 𝜕qkl =0 if j≠ l 𝜕 Pj(qt) 𝜕qkl =kj 2 ( m ∑ i=1 qij )− 1 ∕ 2 when j= l 𝜕P j (qt) 𝜕qkl cij (q ij )=(a kl q ij +b kl ) 2 (6) q t+1 kl =max {0, {qt kl +at( 𝜔 k 𝛾 kl 𝛽k − 1 𝛽k (2a2 klqkl +2aklbkl ) −𝜆t k+𝜆1 l t−𝜆2 l t )}} 1 An interested reader can find a detailed description in [3]. A.Luque-Cerpa et al. 3 Transition P systems withmembrane polarization In this section, we introduce the chosen P system model for designing our solution to the GNE problems. We have chosen the model of transition P systems [14] endowed with membrane polarization, i.e., the membranes have a charge that acts as a state of the membrane. That state controls the evolution of the objects within the membrane. For some definitions about membrane systems and formal languages, we refer the reader to [15, 24]. A transition P system with membrane polarization of degree q≥1 is a tuple where: 1. Γ is an alphabet whose elements are called objects; 2. 𝜇 is a hierarchical tree-like structure; 3. M1,…,Mq are multisets of objects over Γ ; 4. R1,…,Rq are sets of rules of the following form: • [u → v]𝛼 h,u,v∈M(Γ),h∈{1, …,q},𝛼∈{0, +,−} is an object evolution rule; • [ u]𝛼 h →v[] 𝛽 h ,u,v∈M(Γ),h∈{1, …,q},𝛼,𝛽∈{0, +, −} is an send-out communication rule; • u [] 𝛼 h →[v] 𝛽 h ,u,v∈M(Γ),h∈{1, …,q}, h is not the label of the skin membrane ,𝛼,𝛽∈{0, +,−} is an send-in communication rule; 5. 𝜌1,…,𝜌q are partial weak relations between rules from R1,…,Rq , respectively; 6. iout ∈{0, 1, …,q} is the output region (if iout =0 , the output region is the environment of the system). A transition P system with membrane polarization can be seen as a set of q membranes organized in a rooted tree-like structure whose root node is called the skin membrane, each of them having a polarization among 0, + , or −. A configuration of Π in an instant t can be described by the multisets of objects in each membrane in such a moment and the polarization of each membrane and the multiset of objects in the environment, denoted by M0,t ; that is, Ct= ((M1,t,𝛼1,t),…,(Mq,t,𝛼q,t),M0,t) . It can also be described in a more graphical way such that [u]𝛼 h represents that membrane h contains the multiset of objects u and its charge is 𝛼 . If the graphical description contains a membrane h′ such that it is in the membrane h, like [[ ] h � ]h , then h′ is a child membrane of h. The initial configuration of Π is C0= ((M1 ,0 ) , … , (Mq ,0 ) , �) . A configuration is called a halting configuration if no more rules are applicable to it. For the sake of simplicity, we will use the notation Π=(Γ ,𝜇, M1 , … , Mq , (R1 ,𝜌 1) , … , (Rq ,𝜌 q) ,i out), Π=(Γ ,𝜇, M1 , … , Mq , (R1 ,𝜌 1) , … , (Rq ,𝜌 q) ,i out) ℂt to denote specific parts of the P system in such a way that ℂt=𝜇� , where 𝜇′ is a subtree of 𝜇 , and the multisets of objects associated to each membrane h of 𝜇′ at an instant t is denoted by [u]h , where u is a string that represents the multiset Mh,t . An object evolution rule [u → v]𝛼 h , u , v∈M(Γ), h∈{1, …,q}, 𝛼 ∈{0, +,−} is applicable to a configuration Ct if there exists a membrane labeled by h in Π such that it contains the multiset of objects u and its polarization is 𝛼 . Applying such a rule leads to the removal of u from the membrane h and the generation of the multiset of objects v in the membrane h. An send-out communication rule [ u]𝛼 h →v[] 𝛽 h ,u,v∈M(Γ) , h∈{1, …,q},𝛼,𝛽∈{0, +,−} is applicable to a configuration Ct if there exists a membrane labeled by h in Π such that it contains the multiset of objects u and its polarization is 𝛼 . The application of such a rule leads to the removal of u from such a membrane h, the generation of the multiset of objects v in the parent region of h and the change of polarization of such a membrane h from 𝛼 to 𝛽 . An send-in communication rule u [] 𝛼 h →[v] 𝛽 h ,u , v∈M(Γ),h∈{1, …,q},h is not the label of the skin membrane ,𝛼,𝛽∈{0, +,−} is applicable to a configuration Ct if there exists a membrane labeled by h in Π such that its parent membrane contains the multiset of objects u and the polarization of such membrane h is 𝛼 . The application of such a rule leads to the removal of u from the parent membrane, the generation of the multiset of objects v in such a membrane h, and the change of polarization of such a membrane h from 𝛼 to 𝛽 . In addition, if (r1,r2)∈𝜌h , it means that r1 has priority over r2 in the following sense: Rule r2 is applicable only if the remaining objects from membrane h cannot fire rule r1 any more times. If one object can fire more than one rule, then it will be selected in a non-deterministic way. Apart from that, object evolution rules are applied in a maximal parallel way; that is, any object that can fire an applicable rule will fire it. A multiset of rules U is maximal if there is no other multiset of rules U′ in the set of multisets of applicable rules such that U⊂U′ . Send-in and send-out rules are applied in a maximal parallel way as follows: Let r1 and r2 two rules where either r1 and r2 can be a send-in or a send-out rule. Let h1 ( h2 , respectively) be the label of the membrane affected by the rule r1 ( r2 , resp.). Let 𝛼1 ( 𝛼2 , respectively) be the polarization of the left-hand side of the rule (that is, the part of the rule that is to the left of the arrow) of rule r1 ( r2 , resp.), and let 𝛽1 ( 𝛽2 , resp.) be the polarization of the right-hand side of the rule r1 ( r2 , resp.). We say that r1 and r2 are compatible if the following holds: h1=h2∧𝛼1=𝛼2 → 𝛽1=𝛽2 ; that is, two rules that affect the same membrane and can be applied in the same moment t are compatible if and only if the resulting An application ofmembrane computing tohumanitarian relief viageneralized Nash equilibrium polarization by the application of each rule is the same. Each multiset of applicable rules can contain only compatible rules. A transition or a computational step of Π is made by applying the rules in the aforementioned way in all membranes at the same time, and it is denoted by Ct ⇒ ΠCt+1 . A computation of Π is a sequence of configurations C=(C0,C1,…,Cn) such that C0 is the initial configuration of Π and for each t, Ct ⇒ Ct+1 . A computation is halting if n∈ℕ , and Cn is a halting configuration. 3.1 An example ofaP system For the sake of simplicity, we give a simple example to explain the behavior of transition P systems with membrane polarization. Let be a transition P system with membrane polarization of degree 1 where: 1. Γ={a,b,c,d,e} ; 2. 𝜇=[ ] 1 ; 3. M1={a3,d} 4. R1={r1 ≡ [a2 → b]1,r2 ≡ [a → c]1,r3 ≡ [d → e]1} ; 5. 𝜌1= {(r1,r2)} ; 6. iout =0 . The initial configuration can be represented as ℂ0 =[a 3 ,d] 0 1 . Then, the following multisets of rules are applicable: • U1={ r 1, r 2} ; • U2={r1} ; • U3={ r 1, r 3} ; • U4={r1,r2,r3} ; Since we look for maximal multisets of applicable rules, we can observe that U1,U2,U3⊂U4 . Therefore, the multiset U4 will be applied in the first computational step, leading to the following configuration: ℂ1=[b,c,e]1 . 4 Design andfunctioning oftheP system Let us consider the Eqs.3, 4, 5 and 6 defined in Sect.2.1. In this section, a P system that computes the solutions of these equations is analyzed.2 We first need to define a sequence {at} that satisfies Π=(Γ,𝜇,M1,(R1,𝜌1),iout) as indicated in Sect.2.1. The sequence at=0.1 ⋅ bt is defined by {bt}= taking r copies of the element 0.1 ⋅ 1∕2w with r=max{1024, 2w} and w=0, 1, … . This sequence guarantees convergence and can be easily updated and employed using the evolution rules of a P system. In opposition, the implementation of the sequence used in [3] can raise the complexity of the model. The terms at are represented through objects s in the P system. The computation of the P system can be summarized in a loop of three stages, represented in Algorithm1. Algorithm 1 General overview of the P system computation Next, an analysis of the computation of each stage is provided. The result of the analysis is that, for each time step, the number of transition steps in the Initialization and Comparison stages is bounded by a constant. Meanwhile, the number of transition steps in the Update stage is constant for the first 10,240 iterations (1024 identical elements for each of the first 10 values in the sequence at ) and then grows logarithmically with time. This means that the time complexity of the global computation only depends on the number of iterations required by the original algorithm. In the computation analysis, ℂ𝜏 represents the 𝜏 -th configuration of each iteration; that is, ℂ0 is the initial configuration in each iteration. The main difference between iterations is given by the existence of objects s, and this fact is considered along the computation. Stage 1: initialization In this stage, all the necessary objects for updating the values q t + 1 kl ,𝜆t + 1 k ,𝜆1 t+1 l and 𝜆 2 t+1 l are created and distributed to the corresponding membranes. We assume that this stage starts with some objects xk,l,lak,la1,l,la2,l inside the membrane INIT which represent the values qt kl ,𝜆 t k ,𝜆1 t l and 𝜆 2 t l at a time step t. For simplicity, we consider only the objects xk,l , but the behavior of the other objects is analogous. If we only consider a t>0, at→0 and ∞ ∑ t=0 at =∞ { 1, 1, … ,1 ⏟⏞ ⏞⏟⏞ ⏞⏟ 1024 times ,1 ∕ 2, 1 ∕ 2, … ,1 ∕ 2 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 1024 times , … ,1 ∕ 1024, 1 ∕ 1024, … ,1 ∕ 1024 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 1024 times , 1∕2048, 1∕2048, …,1∕2048 ⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏞⏟ 2048 times ,…} 2 A detailed description of the P system can be found in the Appendix A. A.Luque-Cerpa et al. the membranes involved in each step, we have the following configuration: In the next transition step, different copies of each object in the membrane INIT are ejected into the skin membrane. The objects y are applied to coordinate the computation: In one transition step, objects representing the constant part of equations (3) to (6) (before multiplying by at ) are created in the corresponding membranes. For simplicity, only the membranes Qk,l , which represent the Eq.6, are considered. The objects p0 represent originally the term 𝜔 k 𝛾kl 𝛽k , while objects ct0 represent the term −2a kl bkl 𝛽k (see rule RS1,6 in the Appendix A). In one transition step: Now that the charge of the membranes is negative, objects representing qt kl ,𝜆 t k ,𝜆 1t l and 𝜆2t l are inserted into the corresponding membranes. The terms qt kl are represented by the objects p and c0 , while the 𝜆 -terms are represented by objects p0 or n0 . These objects will be used to update the variables using Equations (3) to (6) (see rules RS1,10 to RS1,19 in the Appendix A). In one transition step: In this step, all membranes Qk,l , LAMBk , LAMB1,l and LAMB2,l have the necessary objects to start the Update stage, and all of them are coordinated through the usage of the object y1 . The Initialization stage is then completed in three transition steps. Stage 2: update In this stage, the values q t+1 kl ,𝜆t+1 k ,𝜆1 t+1 l and 𝜆 2 t+1 l are computed following Equations (3) to (6). The computation has two phases: one of computing the parentheses in each equation and multiplying it by at , and one of computing the max function. For simplicity, only membranes Qk,l are considered, but the behavior for membranes LAMBk , LAMB1,l and LAMB2,l is analogous. This stage starts with the configuration ℂ3 . In two transition steps, all terms in the parentheses of Equations (3) to (6) will be represented by objects p0 or n0 depending on their sign: In one transition step, objects y3 will be introduced in the membranes REDUCE, changing their charge to negative. In one more transition step, objects p0 and n0 will be introduced ℂ 0= [[x q k,l k,l ] 0 INIT y0] 0 1 ℂ 1=[x q k,l k,l ,xt q k,l k,l ,xl q k,l 0,k,l ,xl q k,l 1,k,l ,xl q k,l 2,k,l ,y0,k,l,ylk,yl1,l,yl2,l] 0 1 ℂ 2=[[y0,p 𝜅 0 0,ct 𝜅 1 0] − Q k,l x q k,l k,l,xt q k,l k,l,xl q k,l 0,k,l,xl q k,l 1,k,l,xl q k,l 2,k,l] 0 1 ℂ 3=[[y1,p 𝜅 0 0,ct 𝜅 1 1,p q k,l,c q k,l 0] − Q k,l ] 0 1 ℂ 5=[[y3,p𝜅0 0,n𝜅1 +⌊ 2⋅a 2 kl ∕ 𝛽k ⌋ ⋅qk,l 0,pqk,l]− Q k,l ] 0 1 in the membranes REDUCE. For the first 1024 iterations, in the next step objects p0 and n0 will be added, generating objects p or n, and the object y6 will be generated, changing the charge of membranes REDUCE to positive. In the rest, there will be objects s present that will cause the reduction in the number of objects p and n by multiplying them by 10 ⋅ at . There is one more transition step for each object s present, and there are n objects s present for the sequence term at=0.1 ⋅ 1∕2w . Assuming we end with objects p (analogous otherwise), in at least three transition steps we have the configuration: In two more transition steps, the objects y6 will evolve and change the charge of membranes REDUCE back to neutral and the charge of membranes Qk,l,LAMBk , LAMB1,l and LAMB2,l to positive. During these steps, the number of objects p or n that were previously in membranes REDUCE will be multiplied by 0.1 and added to the rest in membranes Qk,l,LAMBk , LAMB1,l and LAMB2,l , updating the terms q t+1 kl ,𝜆t+1 k ,𝜆1 t+1 l and 𝜆 2 t+1 l : In the two next transition steps, if objects p and n are present they cancel each other. If objects p are left, they are ejected to the skin as objects okl , ikl and xt1,k,l , and if objects n are left they are deleted. The Update stage is then completed in a minimum of nine transition steps, and the number of transition steps grows linearly to nineteen for the 10240 first iterations and logarithmically with t for the rest. In practice, a solution with convergence tolerance 10−10 is found in the first 2048 iterations, so we have an effective upper bound of ten transition steps. Stage 3: comparison In this stage, the new values qt+1 kl are compared with the previous values qt kl . If a difference is found, the execution continues. Otherwise, the new values are sent to the membrane OUTPUT, and the computation finishes. This stage starts with the following configuration: In two transition steps, objects xtk,l and xt1,k,l are in membrane COMP to be compared: ℂ ≥8=[[y6[p10⋅at⋅ ( 𝜅0 − 𝜅1 −⌊ 2⋅a 2 kl ∕ 𝛽k ⌋ ⋅qk,l ) ] + REDUCE 0,k,l pqk,l] − Q k,l ] 0 1 ℂ ≥10 =[y8,k,l[[] 0 REDUCE 0,k,l pq t+1 kl ]+ Q k,l ] 0 1 ℂ ≥12 =[ym⋅n 10 ,o(q t+1 k,l) k,l ,i(qt+1 kl ) k,l ,xt(qt+1 kl ) 1,k,l [[] 0 REDUCE 0,k,l ]0 Q k,l ] 0 1 ℂ ≥12 =[ym⋅n 10 ,o(q t+1 k,l) k,l ,i(q t+1 kl ) k,l ,xt(q t+1 kl ) 1,k,l ,xt(q t kl) k,l [] 0 COMP [] 0 INIT [] 0 OUTPUT ] 0 1 ℂ ≥ 14 =[i(q t+1 kl ) k,l [y 12 ,o(q t+1 k,l) 1,k,l ,xt(q t+1 kl ) 1,k,l ,xt(q t kl) k,l ]0 COMP [] 0 INIT [] 0 OUTPUT ] 0 1 An application ofmembrane computing tohumanitarian relief viageneralized Nash equilibrium In the next two transition steps, objects xtk,l and xt1,k,l cancel each other. If at least one of them is present, then y16 evolves to y13 and is ejected as y14 . Otherwise, an object stop is created and ejected. The rest of the objects xtk,l and xt1,k,l are deleted. At this moment, two things can happen. The first one is that in the next three transition steps, the stop object changes the membrane OUTPUT charge to negative, the objects o3,k,l go into it, and the computation finishes after changing again the charge of COMP to neutral: The other possibility is that, in the following two transition steps, the y14 object changes the charge of membrane INIT, the objects ik,l get into membrane INIT, and an object y0 is created, leaving the P system in a state where a new iteration can start again: In this last step, a count0 object is also created. This object will contribute to the update of at . We can see then that this stage is completed in six transition steps except for the last iteration, which is completed in seven steps. The result of this analysis is that the whole computation of the system only depends on the number of iterations required by the original algorithm. The algorithmic time complexity of an iteration of the algorithm is constant for Stage 1, logarithmic with t for Stage 2, and constant for Stage 3. The algorithmic time complexity of the whole computation is then O(t log(t)) . In practice, it is O(t) , as explained in Stage 2. In contrast, the original algorithm has to compute Eqs.3, 4, 5, and 6 at every iteration, so the algorithmic time complexity of the original algorithm is O((m ⋅ n+m+2n) ⋅ t) . Consequently, we present a reduction of a factor of m ⋅ n+m+2n . 5 Experiments Different P system simulators are available, such as P-lingua (MeCoSim) [25–27] or UPSimulator [28, 29]. To perform experiments, MeCoSim has been updated and chosen to implement our P system. The last simulator update allows for the design of transition P systems with membrane polarization. To test the correct behavior of our P system, we have considered several examples taken from [3], where the numerical data of some case studies are provided. Namely, the ℂ ≥ 16 =[o(q t+1 k,l) 3,k,l ,i(q t+1 kl ) k,l ,y 14|| stop [] 0 COMP [] 0 INIT [] 0 OUTPUT ] 0 1 ℂ ≥19 =[i(q t+1 kl ) k,l [] 0 COMP [] 0 INIT [o(q t+1 k,l) k,l ]0 OUTPUT ] 0 1 ℂ ≥18 =[y0[] 0 COMP [x(q t+1 kl ) k,l ,count0]0 INIT [] 0 OUTPUT ] 0 1 authors provide four toy examples and the real-world case study of Hurricane Katrina. All of these examples have been simulated. For the five numerical examples, the convergence tolerance chosen is 10−5 , while for the Hurricane Katrina case study, the tolerance is 10−10 . A solution was found for all the numerical examples in less than 1024 iterations. The Update stage took then nine transition steps for every iteration. For the Hurricane Katrina case study, a solution was found in less than 1100 iterations. The Update stage took then nine transition steps for the first 1024 iterations, and ten transition steps for the rest (see Sect.4). The results of the simulations of the four toy numerical examples can be found in Tables1, 2, 3 and 4, where the obtained values for qij (i.e., the number of items that each NGO i provides to a demand point j) are compared. On the right, the equilibrium values obtained in [3] are shown, and, on the left, the values obtained with our P system. Considering that the true values are rounded to one decimal, it can be observed that the results returned by our P system Table 1 Results of the simulation for Example 1 of [3] qij P system Solution q11 352.50012 352.5 q21 247.50012 247.5 Table 2 Results of the simulation for Example 2 of [3] qij P system Solution q11 352.50012 352.5 q12 452.50004 452.5 q21 247.50012 247.5 q22 347.50004 347.5 Table 3 Results of the simulation for Example 3 of [3] qij P system Solution q11 423.75003 423.8 q12 471.25002 471.3 q13 436.87498 436.9 q21 176.25011 176.3 q22 328.75007 328.8 q23 563.12492 563.1 Table 4 Results of the simulation for Example 4 of [3] qij P system Solution q11 411.25004 411.3 q12 458.75002 458.8 q13 499.37496 499.4 q21 138.75012 138.8 q22 291.25007 291.3 q23 750.62488 750.6 A.Luque-Cerpa et al. coincide perfectly with the true values, even after replacing Eq.2 with Eq.6 as discussed in Sect.2.1. Regarding the Hurricane Katrina case study, the replacement of Eq.2 with Eq.6 has impacted the solutions returned, as can be observed in Table5. However, the solution returned by the P system is a good approximation: the average error between the correct solution and the solution returned by our P system is 1.98%, the median is 0.82%, the maximum error is 7.89%, and only in four cases out of thirty is the error higher than 5%. These experiments support our claim that the contribution of the functions Pj(q) can be ignored without incurring significant errors. This implies that the financial funds due to the visibility of all NGOs at each location have little impact on the final solution of the problem of how to distribute humanitarian relief subject to the constraints stated if they follow the assumptions in [3]. 6 Technical conclusions The simulation of continuous processes by intrinsically discrete models is always a complex task, and each realworld problem needs to be deeply studied to adopt the best possible solution. In this paper, we have considered the minimization question expressed in Eq.2. Such an equation involves two terms with derivatives. After a deep study of the problem, we can conclude that the term with the first derivative is several orders of magnitude lower than the remaining terms. By removing this term, we lose accuracy; however, from a practical point of view, this loss is not significant, as the experiments show. Concerning the second term with a derivative, a solution for a general function cij is not possible, but it can be reached in this case bearing in mind that it can be expressed as a second-grade polynomial equation. From a membrane computing perspective, we would like to emphasize that our design shows that the number of transition steps in the Initialization and Comparison stages is bounded by a constant and the number of transition steps in the Update stage grows linearly for the first 10240 iterations and then grows logarithmically with time. This can be considered as the main contribution of this paper from the theoretical side. Due to the intrinsic massive parallelism of the Membrane Computing devices, the time complexity of the global computation, considered on the basis of the P system steps, only depends now on the number of iterations required by the original algorithm. Besides, the modularity of the design allows us to extend the process with new agents by introducing new modules and only minimal changes in the other modules. Additionally, the object-based approach is similar to agent-based models in the sense that the behaviors of these individuals can be tracked, having a one-toone relationship between objects and real-life resources, and leading to a fine-grained model that can be calibrated for different scenarios. As a final remark, it is worth stressing that, although the use of transition P systems with polarizations to generalized Nash equilibria was introduced in [19], a first detailed reallife use case, specifically applied to model humanitarian relief, is presented here. 7 Final remarks In this paper, we have shown that Membrane Computing is a useful computational paradigm to model humanitarian relief. The main contribution of this paper is twofold. On the one hand, we contribute to highlighting that it is unacceptable that efforts to distribute humanitarian aid are not enough successful due to a lack of coordination, congestion in terms of logistics, and duplication of services. This situation demands an answer from the scientific community, and many theoretical efforts must be made to optimize the resources. One possible solution is to model the problem in terms of GNE, but other optimization methods can be Table 5 Results of the simulation for the use case of Hurricane Katrina in [3] Others Red cross Salvation army Location P system Solution Error P system Solution Error P system Solution Error St.Charles 16.10 17.48 7.89% 29.47 28.89 2.01% 4.35 4.192 3.77% Terrebonne 267.93 267.02 0.34% 410.31 411.67 0.33% 73.38 73.57 0.26% Assumption 47.92 49.02 2.24% 77.59 77.26 0.43% 13.09 12.97 0.93% Jefferson 264.57 263.69 0.33% 405.34 406.68 0.33% 72.30 72.45 0.21% Lafourche 186.55 186.39 0.09% 287.22 287.96 0.26% 51.11 51.18 0.14% Orleans 466.04 463.33 0.58% 710.66 713.56 0.41% 126.63 127.1 0.37% Plaquemines 20.53 21.89 6.21% 37.02 36.54 1.31% 4.37 4.23 3.31% St.Barnard 74.52 72.31 3.06% 120.12 115.39 4.10% 17.13 16.22 5.61% St.James 57.66 58.67 1.72% 92.31 92.06 0.27% 15.77 15.66 0.70% St.John Baptist 16.83 18.2 7.52% 30.60 29.99 2.03% 4.51 4.40 2.50% An application ofmembrane computing tohumanitarian relief viageneralized Nash equilibrium explored. On the other hand, to the best of our knowledge, this is the first paper bridging Membrane Computing with the distribution of humanitarian relief and we have shown that P systems can be a useful tool to model complex situations in this area. As future research, many applications of Membrane Computing to logistics in humanitarian aid can be explored, not only with transition P systems but with other P systems models. In parallel, other optimization methods can be considered to optimize humanitarian aid distribution. As a consequence of the good behavior of the model, it seems reasonable to study specific situations where the general case is not enough to model real-life events, such as abnormal donation patterns, changes in the political landscape, or any other fact that could impact the distribution of humanitarian relief. Appendix AAppendix: Definition oftheP system Let N={1, ..., m} and D={1, ..., n} with m and n the number of NGOs and demand locations respectively. Let P−1=10−p with p∈ℕ be the tolerance for convergence of the algorithm. Let us consider the following transition P system with membrane polarization where the alphabet of objects is given by: The set of membrane labels is given by: Π=⟨Γ,H,EC,𝜇,{wh}h∈H,(R,𝜌)⟩ Γ={y 0 ,y 1 ,y 2 ,y 3 ,y 4 ,y 5 ,y 6 ,y 7 ,y 10 ,y 11 ,y 12 ,y 13 ,y 14 ,y 15 }∪ ∪{y0,k,l,yl0,k,yl1,l,yl2,l|k∈N,l∈D}∪ ∪{y8,k,l,y9,k,l,yla8,k,yla9,k,yla1,8,l,yla1,9,l,yla2,8,l,yla2,9,l|k∈N,l∈D}∪ ∪{xk,l,xtk,l,xt1,k,l,xl0,k,l,xl1,k,l,xl2,k,l|k∈N,l∈D}∪ ∪{lak,laq0,k,l,la0,k,la1,l,laq1,k,l,la1,l,la2,l,laq2,k,l,la2,l|k∈N,l∈D}∪ ∪{p0,n0,ct0,ct1,ct2,c0,c1,p,n,o}∪ ∪{ik,l,lao0,k,lao1,l,lao2,l,ok,l,o0,k,l,o1,k,l,o2,k,l,o3,k,l,o4,k,l,|k∈N,l∈D }∪ ∪{rem,stop,s,s0,countn,u0,un,maxn | n≥1}∪ ∪{sq k,l ,sl k ,sl 1,l ,sl 2,l| k∈N,l∈D} H ={1, INIT,COMP,OUTPUT}∪{Qk,l|k∈ N ,l∈ D }∪{LAMBk| k∈N}∪{LAMB1l|l∈D}∪{LAMB2l|l∈D}∪{REDUCE0,k,l | k∈N,l∈D}∪{REDUCE1,k|k∈N} ∪{REDUCE 2,l |l∈D}∪{REDUCE 3,l |l∈D}. The membrane structure 𝜇 can be defined as follows: • The skin membrane with label 1, inside of which there exist: 1. One membrane with label INIT. 2. One membrane with label OUTPUT. 3. One membrane with label COMP. 4. One membrane with label Qk,l∀ k ∈N ; ∀ l ∈D . Inside of each membrane Qk,l : – One membrane with label REDUCE0, k , l∀k∈ ; ∀l∈ 5. One membrane with label LAMBk ∀k∈ N . Inside of each membrane LAMBk : – One membrane with label REDUCE1,k ∀k∈N 6. One membrane with label LAMB1l∀l∈D . Inside of each membrane LAMB1l : – One membrane with label REDUCE2,l∀l∈D 7. One membrane with label LAMB2l∀l∈D . Inside of each membrane LAMB2l : – One membrane with label REDUCE3,l∀l∈D