Analyzing reciprocity dynamics in supply chains of public goods: a stochastic evolutionary game approach

Introduction

Public goods are a sort of commodity. In The Pure Theory of Public Expenditure, published in 1954, (Samuelson, 1954) an American economist, indicated that public goods are products and services that several people utilize in the same period. Public goods do not have rivals when they are consumed or used and are not excludable in advantage.

A large amount of data on public goods in experimental games, as outlined by Cardenas & Carpenter (2008), confirms that complete rationality in human behavior is not a strong assumption. Chen, Ye & Wang (2012) indicated that individuals do not show complete rationality and choose one of the options when preferences are a concern. Individuals cooperate ubiquitously when, for example, environmental protection grows as an issue. Generally, residents living in rural areas raise money to fix water canals and buy cleaning services when public resources or public goods are not on time. These actions helped form repeated or multiple-stage games of participants since individuals are aware of the result of collective benefits greater than the total of individual ones when cooperation is in use. Hence a synergistic effect occurs. Friedman showed in 1971 that the result of any Nash equilibrium can be constructed in a perfect equilibrium of repeated games (Friedman, 1971) when the Pareto condition dominates the original game. In 1976, Aumann and Shapley suggested the replacement of the “Nash Equilibrium” with “Subgame Perfect Equilibrium”. The reasons explained by theorists why cooperations emerged and maintained based on varied perspectives. Robert (1981) believes that the repeated game of complete information is pertinent to the evolution of the fundamental structure of interactions between agents and proves that cooperation, altruism, revenge, and threat are the outcomes of bounded rationality in real situations. The most proper means to promote cooperation is to employ the robust reciprocity theory (Rand, Dreber & Elingsen, 2009) to expect future interests in a long-term repeated game. Cooperation can be formed and maintained by building a reasonably stout reciprocity mechanism for games on any public goods where cooperative equilibrium cannot be formed.

Evolutionary games present a robust theoretical framework to explain how to promote and keep cooperation. However, players may not always maintain the same rationality when strategies are adjusted in real-life situations. Besides, conventional evolutionary game models cannot reflect the effects of uncertain factors such as information opacity and income volatility. Foster & Young (1990) proposed the Stochastic Stable Strategy (SSS) concept in 1990 when studying stochastic issues in evolutionary games. Then, many scholars proposed different stochastic evolutionary game models (Michihiro, George & Rafael, 1993; Christine et al., 2004; Lorens & Martin, 2006; Wallace & Young, 2015). So far, stochastic games have become a hot research topic for game theorists. In 2018, Ji et al. (2018) studied the game model of the random public goods in finite groups. Li, Zichun & Hui (2020) considered the impact of random perturbation payment on the equilibrium point in the evolutionary game model.

However, few studies investigate the long-term benefits and stochastic issues in the repeated game of public goods. More up-to-date research regarding reciprocity dynamics in supply chains and evolutionary game theory in the literature (Ahmad, Shah & Al-Fagih, 2023; Hao et al., 2022). Therefore, the work and innovation of the article introduce a strong reciprocity mechanism in the long-run benefits of the game based on the repeated games of the manufacturer and seller groups in the supply chains of public goods, proposing a model called the evolutionary game to cooperate between these groups in the supply chain of public goods based on the classical evolutionary game models by adding WGN to reflect random interference in the game process (Sun, Wang & Xue, 2016; Li & Xin, 2017; Sheng & Chen, 2019; Wang & Xia, 2019). The cooperative evolution process of manufacturer and seller group are analyzed to interpret changes in system stability with random disturbances, and finally examining the effects of different situations on stability through numerical simulations. More up-to-date research regarding the latest research can be found (Saeedi et al., 2024; Ghanei, Contreras & Cordeau, 2023; Bilancini et al., 2024; Okada, 2023).

Repeated Game Evolution Model of Public Goods

Repeated Random Evolution Model and Stability of Public Goods Supply

According to the analysis in Section ‘Dynamic System’, when the parameters satisfy A₁>C₁($1-\overline{\omega }$)−δ, A₂>B₂($1-\overline{\omega }$)−δ the conditions, Eq. (1) has a stable equilibrium point at (0,0), but when the system is subject to random disturbance, the question “Is the strategy (0,0) still stable?” needs to be answered. Random disturbance influences stability, so the zero solution stability of stochastic differential equation is conceptualized.

Definition 1 (Cobbl, 1985; Zhang, Xue & Zhou, 2019; Xu, Liu & Qian, 2011; Hu, Huang & Wu, 2008): Let the stochastic process X ={X(t), t ≥0} be the solution to the initial value problem of the following It $\hat{\mathrm{o}}$ differential equation, $\left\{\begin{array}{c}\mathrm{dX}\left(\mathrm{t}\right)=\mathrm{f}\left(\mathrm{t},\mathrm{X}\left(\mathrm{t}\right)\right)\mathrm{dt}+\mathrm{g}\left(\mathrm{t},\mathrm{X}\left(\mathrm{t}\right)\right)\mathrm{d}{\mathrm{B}}_{\mathrm{t}}\forall \mathrm{t}\ge 0\hfill \\ \mathrm{X}\left({\mathrm{t}}_0\right)={\mathrm{x}}_0\hfill \end{array}\right.$

If there is a negative k-moment Liapunov exponent in $X\left(t\right)$, i.e., $\underset{\mathrm{t}\to \mathrm{\infty }}{\overline{lim}}{\mathrm{t}}^{-1}\mathrm{lnE}|\mathrm{X}\left(\mathrm{t}\right){|}^{\mathrm{k}}$<0, then the equation’s solution is said to be stable in terms of the k-expected moment exponent; if lim_t→∞t⁻¹lnE|X(t)|^k>0, then the solution is unstable.

To use Lemma 1 to find the stability of the solution,

Lemma 1 (Cobbl, 1985; Zhang, Xue & Zhou, 2019; Xu, Liu & Qian, 2011; Hu, Huang & Wu, 2008): Let the random process X ={X(t),t ≥0} be the solution of the initial value problem in Eq. (5), (5)${\displaystyle \left\{\begin{array}{c}\mathrm{dX}\left(\mathrm{t}\right)=\mathrm{f}\left(\mathrm{t},\mathrm{X}\left(\mathrm{t}\right)\right)\mathrm{dt}+\mathrm{g}\left(\mathrm{t},\mathrm{X}\left(\mathrm{t}\right)\right)\mathrm{d}{\mathrm{B}}_{\mathrm{t}}\forall \mathrm{t}\ge 0\hfill \\ \mathrm{X}\left({\mathrm{t}}_0\right)={\mathrm{x}}_0\hfill \end{array}\right.}$

where V(t,x) is a continuously differentiable mapping with c₁, c₂ positive real numbers that lead to (6)${\displaystyle {\mathrm{c}}_1|\mathrm{x}{|}^{\mathrm{k}}\le V\left(t,x\right)\le {\mathrm{c}}_2|\mathrm{x}{|}^{\mathrm{k}}.}$ (1) If there exists a positive constant γ such that LV(t,x) ≤ - γV(t,x), the solution of the initial value problem of the It $\hat{\mathrm{o}}$ differential equation is called exponentially stable in terms of the k-expected moment, i.e., $\underset{\mathrm{t}\to \mathrm{\infty }}{\overline{lim}}{\mathrm{t}}^{-1}\mathrm{lnE}|\mathrm{X}\left(\mathrm{t}\right){|}^{\mathrm{k}}$<0.

(2) If there is a positive constant γ that makes LV(t,x) ≥γV(t,x), the solution of the initial value problem of It $\hat{\mathrm{o}}$ differential equation is said to be exponentially unstable to the k-expected moment, that is, $\underset{\mathrm{t}\to \mathrm{\infty }}{\overline{lim}}{\mathrm{t}}^{-1}\mathrm{lnE}|\mathrm{X}\left(\mathrm{t}\right){|}^{\mathrm{k}}>$0.

where LV(t,x) = V _t(t,x)+ V _x(t,x) f(t,x)+ $\frac{1}{2}{\mathrm{g}}^2$(t,x) V _(x,x)(t,x).

According to Lemma 1, sufficient conditions for the exponential stability of the k-expected moment of the zero solution in the stochastic system are attained.

Proposition 1 For the stochastic differential Eq. (3), ${\displaystyle dx=x\left(t\right)\left(\left[A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)\right]y-\left(D_1-B_1-\delta \right)\left(1-\overline{\omega }\right)\right)dt+\sigma x\left(t\right)dw\left(t\right)}$

There exist c₁ = c₂ =1, k =1, γ =1 and Liapunov equation V (t,x) = x(t),

When (1) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$ =0, $\delta \le D_1-B_1-\frac{1}{1-\overline{\omega }}$, or

(2) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$>0, δ<min{$D_1-B_1-\frac{1}{1-\overline{\omega }}$, $\frac{C_1-A_1-1}{1-A_1\overline{\omega }+\overline{\omega }-\overline{\omega }\overline{\omega }}$}or

(3) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$<0, $\delta \ge D_1-B_1-\frac{1}{1-\overline{\omega }}$, the expected moment of the equation’s zero solution is exponentially stable.

When (4) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$ =0, $\delta \ge D_1-B_1+\frac{1}{1-\overline{\omega }}$, or

(5) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$>0, $\delta \ge D_1-B_1+\frac{1}{1-\overline{\omega }}$, or

(6) When $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$<0, $\delta >max\left\{D_1-B_1+\frac{1}{1-\overline{\omega }},1-A_1+C_1\left(1-\overline{\omega }\right)\right\}$, the k-expected moment of the equation’s zero solution is exponentially unstable.

Prove: Considering the stochastic differential Eq. (3)

Assume that the Liapunov function V(t,x) = x(t), c₁ = c₂ =1, k =1 exists which satisfies Eq. (5) in lemma 1, and there is LV(t,x) = f(t,x) = $\mathrm{x}\left(\mathrm{t}\right)\left(\left[{\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)\right]\mathrm{y}\left(\mathrm{t}\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)\right)$. According to Lemma 1, if γ =1 makes Lv(t,x) ≤− V(t,x), then, the zero-solution k-expected moment of the equation is exponentially stable, so

LV(t,x) = f(t,x) = $x\left(t\right)\left(\left[A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)\right]y\left(t\right)-\left(D_1-B_1-\delta \right)\left(1-\overline{\omega }\right)\right)\le -x\left(t\right)$,

Therefore, $\mathrm{x}\left(\mathrm{t}\right)\left(\left[{\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)\right]\mathrm{y}\left(\mathrm{t}\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)+1\right)$ ≤0

As $\mathrm{x}\left(\mathrm{t}\right)$ ∈ [0,1], so long as

$\left[{\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)\right]\mathrm{y}\left(\mathrm{t}\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1\right)-\delta$) $\left(1-\overline{\omega }\right)+1$ ≤0.

Therefore, the following situations are discussed:

(1) When ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)$ =0, there is (D₁ − B₁ − δ)) $\left(1-\overline{\omega }\right)\ge 1$, namely ${\mathrm{D}}_1-{\mathrm{B}}_1-\frac{1}{1-\overline{\omega }}\ge \delta$;

(2) When ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)$>0,

and when $\mathrm{y}\left(\mathrm{t}\right)$ =0 时, there is $\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)$) $\left(1-\overline{\omega }\right)\ge 1$, 即$\delta \le {\mathrm{D}}_1-{\mathrm{B}}_1-\frac{1}{1-\overline{\omega }}.$

and if $\mathrm{y}\left(\mathrm{t}\right)$ =1, ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)$) $\left(1-\overline{\omega }\right)+1$ ≤0, then $\delta \left(1-{\mathrm{A}}_1\overline{\omega }+\overline{\omega }-\overline{\omega }\overline{\omega }\right)+{\mathrm{A}}_1-{\mathrm{C}}_1$+1 ≤0, namely δ ≤ $\frac{{\mathrm{C}}_1-{\mathrm{A}}_1-1}{1-{\mathrm{A}}_1\overline{\omega }+\overline{\omega }-\overline{\omega }\overline{\omega }}$.

Therefore: δ<min{${\mathrm{D}}_1-{\mathrm{B}}_1-\frac{1}{1-\overline{\omega }}$, $\frac{{\mathrm{C}}_1-{\mathrm{A}}_1-1}{1-{\mathrm{A}}_1\overline{\omega }+\overline{\omega }-\overline{\omega }\overline{\omega }}$}, ${\mathrm{D}}_1-{\mathrm{B}}_1-\frac{1}{1-\overline{\omega }}\ge \delta$.

(3) When ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)$<0, then, there is (D₁ − B₁ − δ) $\left(1-\overline{\omega }\right)-1$ ≤0, (D₁ − B₁ − δ)) ≤ $\frac{1}{1-\overline{\omega }}$, $\delta \ge {\mathrm{D}}_1-{\mathrm{B}}_1-\frac{1}{1-\overline{\omega }}$, (D₁ − B₁ − δ) $\left(1-\overline{\omega }\right)-1$>0, that is, ${\mathrm{D}}_1-{\mathrm{B}}_1-\delta >\left(1-\overline{\omega }\right)$.

On the other hand, from lemma 1 , if there is γ =1, such that LV(t,x) = f(t,x) ≥x(t), the k-expected exponential moment of the initial value problem’s solution of It $\hat{\mathrm{o}}$ differential equation is unstable. It can be obtained from

LV(t,x) = f(t,x) = $\mathrm{x}\left(\mathrm{t}\right)\left[{\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)\right]\mathrm{y}\left(\mathrm{t}\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)\ge$x(t)

that

$\mathrm{x}\left(\mathrm{t}\right)\left(\left[{\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)\right]\mathrm{y}\left(\mathrm{t}\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)-1\right)\ge$0

As $\mathrm{x}\left(\mathrm{t}\right)$ ∈ (0,1), so long as

$\left[{\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)\right]\mathrm{y}\left(\mathrm{t}\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0;

Therefore, the following situations are discussed:

(4) When ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)$ = 0, $-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0, (D₁ − B₁ − δ) ≤ $\frac{-1}{1-\overline{\omega }}$, $\delta \ge {\mathrm{D}}_1-{\mathrm{B}}_1+\frac{1}{1-\overline{\omega }}$.

(5) When 当 ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)$>0, $\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)+1$ ≤0, $\delta \ge {\mathrm{D}}_1-{\mathrm{B}}_1+\frac{1}{1-\overline{\omega }}$ ;

(6) When ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)$<0, and when $\mathrm{y}\left(\mathrm{t}\right)$ = 0, there is $\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)+1$ ≤ 0, that is: $\delta \ge {\mathrm{D}}_1-{\mathrm{B}}_1+\frac{1}{1-\overline{\omega }}$.

and when $\mathrm{y}\left(\mathrm{t}\right)$ =1, there is ${\mathrm{A}}_1+\delta \overline{\omega }+\left({\mathrm{D}}_1-{\mathrm{B}}_1-{\mathrm{C}}_1\right)\left(1-\overline{\omega }\right)-\left({\mathrm{D}}_1-{\mathrm{B}}_1-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0, and $\delta +{\mathrm{A}}_1-{\mathrm{C}}_1\left(1-\overline{\omega }\right)-1\ge$0, that is: $\delta \ge 1-{\mathrm{A}}_1+{\mathrm{C}}_1\left(1-\overline{\omega }\right)$;

So: $\delta >max\left\{{\mathrm{D}}_1-{\mathrm{B}}_1+\frac{1}{1-\overline{\omega }},1-{\mathrm{A}}_1+{\mathrm{C}}_1\left(1-\overline{\omega }\right)\right\}$

Proposition 2 is obtained similarly.

Proposition 2 For the stochastic differential Eq. (4)

For $dy=y\left(t\right)\left(\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)\right)dt+\sigma y\left(t\right)dw\left(t\right)$, there are c₁ = c₂ =1, k =1, γ =1 and the Liapunov function V(t,y) = y(t), when

(1) $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$ =0, δ ≤ $D_2-C_2-\frac{1}{1-\overline{\omega }}$, or

(2) $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$>0, δ<min{$D_2-C_2-\frac{1}{1-\overline{\omega }}$, $A_2-B_2\left(1-\overline{\omega }\right)-1$}, or

(3) when $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$<0, δ ≤ $D_2-C_2-\frac{1}{1-\overline{\omega }}$, the expected moment of the zero solution of the equation is exponentially stable when

(4) $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$ =0, $\delta \ge D_2-C_2+\frac{1}{1-\overline{\omega }}$, or

(5) $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$>0, $\delta \ge D_2-C_2+\frac{1}{1-\overline{\omega }}$, or

(6) $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$<0, $\delta >max\left\{D_2-C_2+\frac{1}{1-\overline{\omega }},1-A_2+B_2\left(1-\overline{\omega }\right)\right\}$, the zero solution of the equation is exponentially unstable with the k moment.

To prove, consider the stochastic differential Eq. (4)

Assume that the Liapunov function V(t,y) = y(t), then there is c₁ = c₂ =1, k =1, which satisfies Eq. (6) in Lemma 2, and there is LV(t,y) = f(t,y) = $\mathrm{y}\left(\mathrm{t}\right)\left(\left[{\mathrm{A}}_2+\delta \overline{\omega }+\left({\mathrm{D}}_2-{\mathrm{B}}_2-{\mathrm{C}}_2\right)\left(1-\overline{\omega }\right)\right]\mathrm{x}\left(\mathrm{t}\right)-\left({\mathrm{D}}_2-{\mathrm{C}}_2-\delta \right)\left(1-\overline{\omega }\right)\right)$. From Lemma 2, if there is γ =1, such that Lv(t,y) ≤− V(t,y), then the zero-solution k-expected moment of the equation is exponentially stable, and therefore, there is

LV(t,y) = f(t,y) = $y\left(t\right)\left(\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)\right)\le -y\left(t\right)$,

and therefore, $y\left(t\right)\left(\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)+1\right)$ ≤0

As $\mathrm{y}\left(\mathrm{t}\right)$ ∈ (0,1), so long as

$\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)+1$ ≤0

Therefore, the following situations are discussed.

(1) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$ = 0, there is $-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)+1$ ≤0, namely: $D_2-C_2-\delta \ge \frac{1}{1-\overline{\omega }}$, then ${\mathrm{D}}_2-{\mathrm{C}}_2-\frac{1}{1-\overline{\omega }}\ge \delta$ ;

(2) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$>0,

Suppose that x(t) = 0, if there is $-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)+1$ ≤ 0, then $D_2-C_2-\frac{1}{1-\overline{\omega }}$ ≤0.

When x(t) =1, if there is $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)+1$ ≤0, then there is $A_2-B_2\left(1-\overline{\omega }\right)-1-\delta \ge$0, that is: δ ≤ $A_2-B_2\left(1-\overline{\omega }\right)-1$.

So δ ≤min{$D_2-C_2-\frac{1}{1-\overline{\omega }}$, $A_2-B_2\left(1-\overline{\omega }\right)-1$}

(3) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$<0, $\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)\ge$1, namely: $D_2-C_2-\frac{1}{1-\overline{\omega }}\ge \delta$

On the other hand, according to lemma 1 , if there is γ =1, such that LV(t,x) = f(t,x) ≥x(t), the k-expected moment exponent of the initial value problem solution of It $\hat{\mathrm{o}}$ differential equation is said to be unstable. It can be obtained from

LV(t,y) = f(t,y) = $y\left(t\right)\left(\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)\right)\ge$y(t).

That $y\left(t\right)\left(\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)\right)-1\ge$0

As $\mathrm{y}\left(\mathrm{t}\right)$ ∈ (0,1), so long as

$\left[A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)\right]x\left(t\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0

Therefore, the following situations are discussed.

(4) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$ =0, $-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0, namely, $D_2-C_2+\frac{1}{1-\overline{\omega }}\le \delta$.

(5) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$>0, $-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0, namely, $D_2-C_2+\frac{1}{1-\overline{\omega }}\le \delta$.

(6) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$<0,

Let x(t) =0, and there is $-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)-1\ge$0, that is: $D_2-C_2+\frac{1}{1-\overline{\omega }}\le \delta$.

When x(t) =1, and there is $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)-\left(D_2-C_2-\delta \right)\left(1-\overline{\omega }\right)-1\ge 0$, so ${\mathrm{A}}_2-{\mathrm{B}}_2\left(1-\overline{\omega }\right)+\delta -1\ge 0$, that is: $\delta \ge 1-A_2+B_2\left(1-\overline{\omega }\right)$

and therefore, δ ≥max{ $D_2-C_2-\frac{1}{1-\overline{\omega }}$, $1-A_2+B_2\left(1-\overline{\omega }\right)$}.

From Propositions 1 and 2, if the parameters satisfy any one of the exponential stability conditions 1 through 3 of the zero solution k-expected moment, the stochastic dynamic system’s equilibrium point (0,0) is constituted of Eqs. (3) and (4), respectively, and is exponentially stable in terms of the k-expected moment, that is, (0,0) is still an evolutionarily stable strategy under random disturbances. Likewise, if the parameters satisfy any one of the conditions of Propositions 4 through 6, the equilibrium point of the system (0,0) for the k-expected moment is exponentially unstable.

Numerical Simulation

To describe intuitively the differential equations’ evolution process, the dynamic evolution process is simulated for distinct parameters. From the income matrix, assuming that the incomes of the manufacturer and seller groups in the repeated game of public goods take different parameter values. Then, the numerical simulation of the evolutionary game process is executed using the combination of these different parameters. The following two situations are mainly considered for numerical simulation:

Scenario 1 In the stable evolution process of cooperative interaction, Set A₁ = 2, δ = 0.4, $\overline{\omega }$ = 0.5, D₁ = 3, B₁ = 0, C₁ = 4, A₂ = 4.5, D₂ = 3.8, B₂ = 2.5, C₂ = 1, then the parameters satisfies

Proposition 1 and Condition 2) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$>0, δ = 0.4<min{$D_1-B_1-\frac{1}{1-\overline{\omega }}$, $\frac{C_1-A_1-1}{1-A_1\overline{\omega }+\overline{\omega }-\overline{\omega }\overline{\omega }}$}=1 and

Proposition 2 and Condition 2) $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$>0, δ = 0.4<min{$D_2-C_2-\frac{1}{1-\overline{\omega }}$, $A_2-B_2\left(1-\overline{\omega }\right)-1$}= 0.8, 且 A₁>C₁($1-\overline{\omega }$) −δ, A₂>B₂($1-\overline{\omega }$) −δ.

Then, the equation’s zero-solution k-expected moment is stable exponentially, and the stochastic system will evolve from (1,1) to the stable strategy (0,0).

Numerical simulation 1: Numerical simulation is carried out according to the parameter values in case 1, and the random interference intensity is taken as σ =0, 0.25, 0.75, and 2, respectively. The cooperative strategy’s evolution process is simulated for the manufacturer and seller groups and is presented in Fig. 1.

1. Figures 1A shows that when there is no random disturbance in the system, the strategies of the manufacturer and seller groups evolve to a stable equilibrium state of cooperative strategy after some time passes. This is mainly because their subjective cognition and information mastery rationally limit the players. Still, manufacturers and sellers have gradually mastered commodity and shopping characteristics and become more rational through market operation. Rational players can only benefit from cooperation.

2. Figures 1B depict that as the system is subject to less random disturbance, the cooperation strategy of the manufacturer and seller groups will fluctuate and then evolve to a stable equilibrium state of cooperation. Due to the players’ rational limitations in subjective cognition and information mastery, there will be inconsistent strategic judgments when the system is affected by uncertain and random events, such as emergencies, thus resulting in inconsistent decision-making and volatile cooperation strategy. This also aligns with objective facts since rational players can benefit only through cooperation.

3. Figures 1C and 1D show that when the random interference is intense, the cooperation strategy of the two groups fluctuates wildly and will deviate from the cooperative strategy (1,1) and evolve to the equilibrium steady state of betrayal policy (0,0). This is because the entirely rational players are subject to greater random interference, making them unable to control their loss of income from random interference such as emergencies. The greater the strength of random interference, the faster the rational players choose a betrayal strategy, and only in this way can the players reduce the losses caused by random interference and obtain the maximum profit.

Numerical simulation 2: In case 1, only the strong reciprocity coefficient, δ = 0.6, in the repeated game is changed, and the values of other parameters and the strength of random interference still take σ = 0, 0.25, 0.75, and 2, respectively. The cooperative strategy’s evolution process is simulated for the manufacturer and the seller groups and can be seen in Fig. 2.

(4) When Figs. 1 and 2 are compared, the larger the strong reciprocity coefficient δ is, the more the evolution of the cooperation strategy fluctuates with the same random interference intensity, which indicates that the manufacturer group and the seller group pay more attention to the impact from the strong reciprocity coefficient on their earnings when making decisions.

Scenario 2 The evolution process of cooperative interaction is unstable, set A₁ = 3, δ = 0.4, $\overline{\omega }$ = 0.5, D₁ = 1, B₁ = 4, C₁ = 4.4, A₂ = 3, D₂ = 1.5, B₂ = 4.5, C₂ = 5, then the parameters satisfies

Proposition1 and Condition 6) $A_1+\delta \overline{\omega }+\left(D_1-B_1-C_1\right)\left(1-\overline{\omega }\right)$<0, $\delta >max\left\{D_1-B_1+\frac{1}{1-\overline{\omega }}\right\}$, $1-A_1+C_1\left(1-\overline{\omega }\right)$} and

Proposition 2 and Condition 6) When $A_2+\delta \overline{\omega }+\left(D_2-B_2-C_2\right)\left(1-\overline{\omega }\right)$<0, $\delta >max\left\{D_2-C_2+\frac{1}{1-\overline{\omega }},1-A_2+B_2\left(1-\overline{\omega }\right)\right\}$, 且 A₁>C₁($1-\overline{\omega }$) −δ, A₂>B₂($1-\overline{\omega }$) −δ,

Where the system’s zero-solution k-expected moment is unstable exponentially, the random system will not evolve from (1,1) to the stable policy (0,0).

Numerical simulation 3: Numerical simulation is carried out according to the parameter values in case 2, and the random interference intensity is taken as σ = 0, 0.25, 0.75, and 2, respectively. The cooperative strategy’s evolution process is simulated for the manufacturer and seller groups and is presented in Fig. 3.

(5) Figures 3A and 3B indicate that when there is no or a little random interference in the system, the strategies of the manufacturer and seller groups will soon evolve to a stable equilibrium state of cooperative strategy. This is mainly because rational game players can quickly grasp the market operating rules, and cooperation will result in the maximum benefit.

(6) Figures 3C and 3D show that the larger the random interference intensity, the larger the fluctuation of the cooperation strategy of the two groups would be, and it is unable to evolve from state (1,1) to (0,0). As the parties of the game cannot predict the market environment that randomly occurring factors will impact, the decision-making processes of the manufacturers and seller groups are always swaying.

Numerical simulation 4: In case 3, only the strong reciprocity coefficient δ = 0.6 in the repeated game is changed, and the values of other parameters and the strength of random interference still take σ = 0, 0.25, 0.75, 2, respectively. The cooperative strategy’s evolution process is simulated for the two groups and is depicted in Fig. 4.

(7) When Fig. 3 and 4 are compared, the larger δ, the more significant fluctuation of the cooperation evolution strategy in the case of the same random interference intensity, which means that the two groups pay more attention to the effect of the strong reciprocity coefficient on their earnings when decisions are made, so their hesitation is more prominent.

Conclusion

This article investigates the evolution process of manufacturers’ and sellers’ behavior strategies in the public goods supply chain. Due to the large period and space in public goods projects and the many factors involved, various uncertainties exist in the operation of public goods. To fully reflect the real state of the system, the manuscript begins with the repeated game between the manufacturer group and the seller group in supply chains of public goods. It introduces a strong reciprocity mechanism for in-game revenue. It also builds a model called the evolutionary game for the cooperative strategy of the manufacturer and the seller groups by implementing the conventional evolutionary game method. Then, WGN is added to establish a stochastic dynamic system, and the k-moment of the cooperative strategy is simulated to investigate the evolution process.

Research shows that the cooperation strategy of manufacturers and sellers will fluctuate and then stabilize after some time passes when random interference occurs. The players will adopt a cooperative strategy when the strong reciprocity coefficient of betrayal is small. On the other hand, when it is large, the players will adopt a betrayal strategy. Thus, the greater the random interference intensity, the faster it evolves to equilibrium. Namely, the stronger the reciprocity coefficient and the strength of random interference to the system, the greater the changes in selecting strategies when the decision-making processes of players in the repeated game of public goods are under consideration.

According to the research results, the following suggestions are put forward to manage the supply chains of public goods:

(1) To overcome the influence of random factors not directly related to the system but external factors, manufacturers and sellers should pay close attention to the impact of market policies and emergencies and reduce the influence of uncertainties stemming from external factors. Currently, the exchanges and cooperation between manufacturers and sellers need to be strengthened for common information sharing and give benefits to the people to maximize the supply chain’s overall revenue in line with the principle of cooperating for a win-win situation.

(2) To make the supply chain of public goods run stably, it is necessary to strengthen the behavior supervision of manufacturers and sellers through the effect of a strong reciprocity policy, to avoid one party choosing a betrayal strategy to maximize individual benefits.

The limitation of the research is that WGN is used to represent randomness. However, other types of randomness could be more realistic in modeling real cases. Instead of using real datasets, we implement simulation data. Even though it is out of the scope of the research, real datasets can provide more reliable results.

Supplemental Information