# Dynamics of a Heterogeneous Constraint Profit Maximization Duopoly Model Based on an Isoelastic Demand

Authors: S. S. Askar; A. Ibrahim; A. A. Elsadany
Journal: Complexity
 (2021)

---

## Abstract

A Cournot duopoly game is a two-firm market where the aim is to maximize profits. It is rational for every company to maximize its profits with minimal sales constraints. As a consequence, a model of constrained profit maximization (CPM) occurs when a business needs to be increased with profit minimal sales constraints. The CPM model, in which companies maximize profits under the minimum sales constraints, is an alternative to the profit maximization model. The current study constructs a duopoly game based on an isoelastic demand and homogeneous goods with heterogeneous strategies. In the event of sales constraint and no sales constraint, the local stability conditions of the Cournot equilibrium are derived. The initial results show that the duopoly model would be easier to stabilize if firms were to impose certain minimum sales constraints. Two routes to chaos are analyzed by numerical simulation using 2D bifurcation diagram, one of which is period doubling bifurcation and the other is Neimark–Sacker bifurcation. Four forms of coexistence of attractors are demonstrated by the basin of attraction, which is the coexistence of periodic attractors and chaotic attractors, the coexistence of periodic attractors and quasiperiodic attractors, and the coexistence of several chaotic attractors. Our findings show that the effect of game parameters on stability depends on the rules of expectations and restriction of sales by firms.


---

## Body

## 1. Introduction

In 1838, Antoine Augustin Cournot introduced the concept of duopoly in his book titled Researches on the Mathematical Principles of the Theory of Wealth. He considered the situation of two firms competing in a market by producing a homogeneous product that would maximize their profits. A number of studies on Cournot games has been conducted for two (duopoly) or more players, phrased as dynamic games on discrete time domains [1–6]. Most discrete time oligopoly dynamics that have been considered in the past two decades were based on the single objective of maximizing profit of firms [7–10]. More recently, the stable, periodic, and chaotic dynamic behaviors of these markets can be analyzed on the basis of these models [11–15].From an economic point of view, multiobjective games have the advantage of being more realistic in many cases. Typically, one maximizes profits in oligopoly models as a single objective. This is not realistic since many cases have more than one objective. Baumol [16] proposed an alternative model called constraint sales maximization where a business maximizes its sales on the basis of a minimum constraint of profit. The idea for Baumol’s model was based on his experience as a business consultant, which he found that businesses maximize sales rather than profits. Fisher [17, 18] concluded that maximizing profits which subject to minimum sales constraints are more rational for a business, thus leading to a new model called constraint profit maximization. Since then, there have been a few studies on the dynamics of oligopoly models with constraints. The earliest was by Kamerschen and Smith [19], who showed that the stability of their duopoly model depends on the production cost and the coefficients in the linear demand function. In their work, by including risk minimization into their profit maximization oligopoly model, Ahmed et al. [20] showed that giving more weight to risk minimization decreases profit. In another work, Ahmed et al. [21] showed that a multiobjective oligopoly model is more stable than a single objective one. Mert [22] showed that, in a duopoly model under a hybrid of sales and profit maximization goal, Nash equilibrium occurs at sales maximization conditions rather than profit maximization conditions. More recently, Ibrahim [23] derived the local stability conditions for a Nash equilibrium of duopoly model subject to minimum sales constraint.Decision mechanism plays an important role in the output adjustment process of a firm. Some common mechanisms are naive expectation, adaptive expectation, bounded rationality, and local approximation. All the abovementioned studies consider the case of homogeneous players, where the duopolists or oligopolists adopt the same decision mechanism. However, in practical, a more realistic assumption is the players adopting heterogeneous decision mechanisms, which characterized the works by Léonard and Nishimura [24], Den Haan [25], Agiza and Elsadany [26, 27], and Tramontana and Elsadany [28]. Dubiel-Teleszynski [29] studied a heterogeneous Cournot game with nonlinear cost function and explored the nonlinear dynamics of this game. Tramontana [30] and Ding et al. [31] analyzed a heterogeneous duopoly game with isoelastic and linear demand function, respectively. Both works showed the existent of two different routes to complicated dynamics. Cavalli et al. [32] investigated a dynamic Cournot duopoly with heterogeneous players based on a local monopolistic approach versus a gradient rule with endogenous reactivity local monopolistic approach. Pecora and Sodini [33] analyzed the stability switching curves in the heterogeneous Cournot duopoly with delay.The literature review reported above has shown that the issue of profit maximization for firms has been discussed without constraints in dynamic situation. Recently, certain games models have been investigated and have given rise to complex dynamics with some constraints, such as global analysis and multistability of Cournot duopoly market game based on consumer surpluses [34, 35]. The effect of model parameters on game dynamics, coexisting attractors, and global dynamics of constrained competition games by both corporate social responsibility (CSR) and social welfare (SW) has been studied in [36, 37].The goal of this work is to reconsider the CPM Cournot duopoly game by providing an isoelastic demand function with heterogeneous players and to develop the Tramontana model [30] on the basis of a minimum sales constraint. In this paper, we analyze a similar duopoly game to Tramontana [30], but the profit maximization is subject to minimum sales constraint. The dynamic equilibrium-point structure reflects corresponding economic explanations. Jury stability criterion [38] and numerical simulation provide local stability of boundary equilibrium and Nash equilibrium for obtaining the game’s internal complexity. The theoretical stability of the equilibrium point and the derivation of critical curves of noninvertiable map are studied, and the game’s dynamic behavior is simulated numerically. We demonstrate that the sequential change in game parameters not only causes the stability of the game to collapse but also increases the complexity of consequent game behaviors such as the number of coexisting attractors, the critical bifurcation of attractors, and the global bifurcation of the attraction basin.The paper is organized as follows. In Section2, a duopoly model and the nonlinear system describing the dynamics of the productions of the firms are described. In Section 3, the local stability conditions for the Nash equilibrium are determined. In Section 4, we present the numerical simulations as well as an analysis of the local bifurcations, basin of attractions, critical curves, and route to complex dynamics. Finally, some remarks are presented in Section 5.


## 2. Model

We consider a market dominated by two firms producing perfect substitute goods. Letqi be the output of firm i=1,2 and Q denote the total output of the two firms. We recall an isoelastic demand function, which is based on the Cobb–Douglas utility function governing the market [39]:(1)p=1Q=1q1+q2,where p is the goods price. Let the cost function for each company be(2)Ciqi=ciqi.Now, each company’s profit is(3)πiq1,q2=qiQ−ciqi.i considers Si to be a constant constraint of sales. The firm i is maximizing its profit subject to its sales constraint:(4)maximizesqiQ−ciqisubject toqiQ≥Si,which is equivalent to maximizing the payoff function as the objective function of firm i:(5)Li=qiQ−ciqi−μiqiQ−Si,where μi is positive parameter and associated with the sales constraint. By differentiating the payoff function Li with respect to qi, we obtain(6)∂Li∂qi=q−i1−μiQ2−ci=0,where q−i represents the output of the other duopolist. The solution of (6) gives the firm’s i reaction function as follows:(7)qi=q−i1−μici−q−i.The response of the firmi to the output of the other firm is the response function.We consider the situation where each firm adopts different mechanisms in adjusting its output in each time period. Following Tramontana and Elsadany [28] and Tramontana [30], we assume firm 1 to have bounded rationality, where an increase or decrease of its output is dependent on the marginal payoff function of the last period. This adjustment process is given by(8)q1t+1=q1t+αq1t∂L1∂q1,where α>0 represents the speed of adjustment.For firm 2, we assume that it naively expects that the output of its rival is the same as the last period, i.e.,q1et+1=q1t. Given this assumption, using the best reaction function in (7), firm 2 will determine its output in period t+1 using the following the adjustment process:(9)q2t+1=q1t1−μc2−q1t.The dynamics of the outputs of the two firms are expressed by the discrete dynamical system as follows:(10)Tq1,q2:q1t+1=q1t+αq1tq2t1−μ1q1t+q2t2−c1,q2t+1=q1t1−μ2c2−q1t,where 0<μ2<1. When μ1=μ2=0, system (10) is reduced to the game considered by Tramontana [30].System (10) is a noninvertible two-dimensional map whose iteration determines competing firms’ trajectories. We investigate the effect of parameters on this system’s dynamics. To discuss the local qualitative behaviors of the solutions of system (10), first, we need to investigate the local stability and bifurcation of equilibrium points of the game (10) between these two firms.


## 3. Nash Equilibrium and Local Stability

The equilibrium points of system (10) are obtained through the fixed point conditions q1t+1=q1t=q1∗ and q2t+1=q2t=q2∗ in the discrete system (10), and we obtain(11)αq1∗q2∗1−μ1q1∗+q2∗2−c1=0,(12)q1∗1−μ2c2−q1∗−q2∗=0.Solving equations (11) and (12) simultaneously yields the following Nash equilibrium point:(13)E=q1∗,q2∗=c21−μ21−μ12c11−μ2+c21−μ12,c11−μ11−μ22c11−μ2+c21−μ12,where E is a positive equilibrium provided that μi<1,i=1,2. The Nash equilibrium E of game (10) is a unique fixed point and is the same equilibrium in Tramontana [30] when μ1=μ2=0. The Jacobian matrix of system (10) is(14)J=1−αc1+α1−μ1q22−q1q2q1+q23α1−μ1q12−q1q2q1+q23121−μ2c2q1−10,with characteristic polynomial(15)fλ=λ2−trJEλ+detJE,where trJ and detJ are the trace and determinant of the Jacobian matrix in (14), respectively.The required condition of the local asymptotic stability of Nash equilibrium point is that the eigenvalues of the corresponding Jacobian matrix are inside the unit circle. According to the Jury criterion [34], the local stability conditions of the Nash equilibrium point are as follows:(16)f−1=1+trJE+detJE>0,f1=1−trJE+detJE>0,detJ<1.If only two of the three inequalities in (16) hold, one of the following three behaviors occurs [38]:(i)
A period doubling bifurcation is obtained whenf−1=0(ii)
A transcritical or fold bifurcation is obtained whenf1=0(iii)
A Neimark–Sacker bifurcation is obtained whendetJ<1Proposition 1.
The Nash equilibrium point (13) is stable provided that α<minαf,αns, where(17)αf=41−μ2c11−μ2+c21−μ11−μ22c12+61−μ11−μ2c1c2−31−μ12c22,αns=21−μ2c11−μ2+c21−μ231−μ12c22−21−μ11−μ2c1c2−1−μ22c12.Proof.
The proof depends on the Jacobian matrix of map (10) given in (14), which when evaluated at the Nash equilibrium point (13) becomes(18)JM=1−α2c1c21−μ1c11−μ2+c21−μ1αc21−μ1c21−μ1−c11−μ21−μ2c11−μ2+c21−μ1−12c11−μ2+3c21−μ1c21−μ10,and the trace and determinant of matrix (18) are(19)trJE=1−α2c1c21−μ1c11−μ2+c21−μ1,detJE=c21−μ1−c11−μ2c11−μ2+3c21−μ121−μ2c11−μ2+c21−μ1.
Substituting (19) in (16) and after some algebraic calculations, the first condition of (16) is satisfied if 3c21−μ1>c11−μ2, while the other two conditions are reduced to α<αf and α<αns. This completes the proof.Proposition 2.
The Nash equilibrium point (13) can be unstable by either period-doubling or Neimark–Sacker bifurcation if(20)αns≤αf⟺131−μ11−μ2c2≤c1≤31−μ11−μ2c2,αns≥αf⟺c1≤131−μ11−μ2c2∪c1≥31−μ11−μ2c2.Proof.
The proof is straightforward and is obtained by reducing the expressionαns−αf.
In Section4, the features of bifurcation and chaos are demonstrated in depth through local analysis numerical simulation such as chaotic attractors, initial sensitivity, intermittent chaos, and multistationary properties. The hidden complexity of the game is studied using chaos and bifurcation theories to reveal the complexity of competition between two heterogeneous firms.


## 4. Numerical Simulation of Dynamic Game Behaviors

The numerical simulations in this section show some insights about the local stability of the Nash equilibrium point (13) and confirm our results in Section 2. In fact, we will see that the dynamics of map (10) becomes more complex due to the map’s trajectories behaving differently when there is a slight change in the initial datum. For example, by fixing the parameters c1=0.35,c2=0.2,μ1=0.1, and μ2=0.5, we get αf≈5.8 and αns≈99.3, which means that αf<αns, and hence, the Nash equilibrium becomes unstable due to period doubling bifurcation. Interestingly, at the initial data x0,1,x0,2=0.20,0.40 and x0,1x0,2=0.40,0.11, we report two different bifurcation diagrams corresponding to c1=0.35,c2=0.2,μ1=0.1, and μ2=0.5 on varying the parameter α. Figure 1(a) shows these bifurcation diagrams, where the blue one is for x0,1,x0,2=0.20,0.40 while the red one is for x0,1,x0,2=0.40,0.11. As we can see, a period 2-cycle arises even before αf≈5.8, particularly at α≈5.69. This indicates that different attractors may coexisted on varying the parameter α. Such attractors are important and we give detailed investigations about them later. Now, we focus and report some attractors of map (10) at the initial datum x0,1,x0,2=0.40,0.11 by keeping the other parameter values fixed. At the bifurcation threshold αf≈5.8, we obtain(21)JM=−1.05860.029408−1.98610,whose eigenvalues are λ1=−1.0002 and λ2=−0.058395, both of which are real with one having an absolute value greater than 1. This means that the Nash equilibrium becomes unstable due to period doubling bifurcation.Figure 1
The fixed parameters arec1=0.35,c2=0.2,μ1=0.1, and μ2=0.5. (a) Bifurcation diagram of q1 with respect to α at the initial data x0,1,x0,2=0.40,0.11 (red color) and x0,1,x0,2=0.20,0.40 (blue color). (b) Phase plane for period 2-cycle at α=5.8. (c) Phase plane for period 4-cycle at α=6.58. (d) Phase plane for period 8-cycle and 16-cycle at α=6.738 and α=6.77, respectively. (e) Phase plane of the four pieces chaotic areas at α=6.78 and α=6.82, respectively. (f) Phase plane of the one-piece and two-piece chaotic areas.
(a)(b)(c)(d)(e)(f)Figure1(b) shows that the map’s dynamic begins with the period 2-cycle, which is represented by squares while the Nash equilibrium point is represented by a circle. In order to investigate the evolution of the map’s dynamic, we have carried out further numerical simulations, as shown in Figures 1(c)–1(f). In Figure 1(c), we present a stable period 4-cycle at α=6.58. This dynamic is followed by a stable period 8-cycle and a stable period 16-cycle which occurred at α=6.738 and α=6.77, respectively (Figure 1(d)). After that and particularly at α=6.773, a stable period 32-cycle arises. We should highlight that the basin of attraction of those stable period cycles becomes more complicated as we approach period 32-cycle. Further increase in the parameter α to 6.78 gives rise to a dynamic of four unconnected areas that turns into four connected areas as α goes to 6.782 (Figure 1(e)). Moreover, the dynamic of the map then goes to an unstable period 20-cycle α=6.787 followed by four pieces of chaotic areas α=6.788 which turn into an unstable period 12-cycle α=6.795. Eventually, the map’s behavior changes into a two-piece chaotic at α=6.85 which then finally evolved into a one-piece chaotic attractor (Figure 1(f)).All numerical simulations performed so far was at parameter values that satisfied the second condition of (16). Conversely, we assume now parameter values that satisfy the first condition. Let us assume the following parameter values, c1=0.1,c2=0.22,α=2,μ1=0.3, and μ2=0.7 with the same initial datum x0,1x0,2=0.40,0.11. Substituting in (18), one obtains(22)JM=0.665220.69188−1.59740,whose eigenvalues are λ1,2=0.33261±0.99729i which are complex and have an absolute value greater than 1. This means the stable Nash equilibrium point will become unstable due to Neimark–Sacker bifurcation on varying the parameter α above the critical value αns. In Figure 2(b), at α=7.10, the attractor becomes a stable spiral point (red color), and as the parameter α increases, the stable spiral enlarges in size (yellow color). The Nash equilibrium point is represented by circle in the previous figure. As α continues to increase, the spiral becomes larger before changing into an attracting invariant closed curve with rough selvedge due to Neimark–Sacker bifurcation. These spiral and selvedge are plotted in Figure 2(c) at α=7.17 (red color) and α=7.18 (blue color), respectively. At α=7.19 and α=7.195, the attracting invariant closed curves becomes larger with some rough edges as shown in Figure 2(d). It is clear from the close values of the parameter α how quasi-periodic motions are produced around Nash point due to this type of bifurcation. Increasing α furthermore gives rise to continued invariant closed curves around the Nash point, but their selvedges begin to vanish as plotted in Figure 2(e) (blue color). After that and particularly when α approaches 7.847, the closed invariant curve is turned into a stable period 7-cycle, as illustrated by the squares in Figure 2(f). This figure analyzes the basin of attraction of this cycle where the cyan color refers to the basin of Nash point, while the yellow color represents the basin of the cycle. When we increase the parameter α until it is equal to 7.894, the previous cycle becomes a closed invariant curve.Figure 2
The fixed parameter values arec1=0.1,c2=0.22,α=2,μ1=0.3, and μ2=0.7. (a) Bifurcation diagrams of q1and q2 on varying the parameter α at the initial datum (the red color) x0,1,x0,2=0.40,0.11. (b) Phase plane of a stable spiral point. (c) Phase plane of a stable spiral point at α=7.17 and α=7.18. (d) Phase plane of an invariant closed curve at α=7.19 and α=7.195. (e) Phase plane of an invariant closed curve at α=7.23 and α=7.30. (f) Stable period 7-cycle and its basin of attraction at α=7.847.
(a)(b)(c)(d)(e)(f)Other invariant closed curves that are different from the previous ones are constructed for the map atα=8.10 and α=8.12, as shown in Figure 3(a). Following this dynamic situation, in Figure 3(b), we obtain an interesting behavior of the map after increasing αto 8.143, that is, the coexistence of multiple stable period 10-cycle. Increasing the parameter α again changes the period 10-cycle into a chaotic attractor which then becomes a stable period 23-cycle at α=8.149. In Figure 3(c), we plot this period 23-cycle with its complicated basin of attraction. The gray color in this figure indicates divergent trajectories. Moving above α=8.149 and, particularly, in the interval 8.150,8.154, the dynamic evolution of map (10) will be mutual between period 23-cycle and chaotic attractors. Consequently, a stable period 26-cycle is formed around the Nash point in Figure 3(d) at α=8.1565, and above this value, the dynamic evolution of map (10) becomes a chaotic attractor, as evidenced in Figures 3(e) and 3(f).Figure 3
The fixed parameter values arec1=0.1,c2=0.22,α=2,μ1=0.3, and μ2=0.7. (a) Phase plane of an invariant closed curve at α=8.10 and α=8.12. (b) Basin of attraction of multiple stable period 10-cycle at α=8.143. (c) Stable period 23-cycle and its basin of attraction at α=8.149. (d) Stable period 23-cycle and its basin of attraction at α=8.1565. (e) Phase plane of chaotic attractor at α=8.15651. (f) Phase plane of chaotic attractor at α=8.158.
(a)(b)(c)(d)(e)(f)So far, we have discussed the complex dynamics of map (10) by varying the parameter α and keeping the other parameter values fixed. Now, we are going to analyze the behavior of the map when either μ1 or μ2 varies and the other parameters become fixed including α. Let us start our investigation by assuming the following parameter values, c1=0.35,c2=0.2,μ1=0.1, and α=5, while μ2 is considered as the bifurcation parameter.In Figure4(a), a period doubling bifurcation diagram is reported on varying μ2. Increasing μ2 further gives rise to a cascade of period doubling bifurcations ensued by periods of higher periodicity appeared and routes to chaos are formed. A confirmation of existence of chaos is the largest Lyapunov exponent corresponding to the bifurcation parameter μ2. As can be seen in Figure 4(b), the largest Lyapunov exponent is negative for μ2<0.7247. Actually, the first period doubling occurs at μ2=0.648. It is plotted in Figure 4(c) with its basin of attraction at the same parameter values and μ2=0.648. The gray color denotes the divergent and unfeasible trajectories B∞, while the white color indicates nonconvergent points. As μ2 increases, the period 2-cycle continues to appear till μ2 approaches 0.721 at which a period 6-cycle arises. Increasing the parameter μ2 further to 0.7245 gives rise to a period 12-cycle. The basin of attraction of those two cycles is given in Figure 4(d). Shifting this parameter to 0.725 causes six-piece chaotic attractors to emerge. The phase plane of these attractors is given in Figure 4(e) with an interesting behavior, that is, the dynamic of map (10) changes from a chaotic attractor to a stable period 2-cycle whose basin is given in the same figure.Figure 4
The fixed parameter values arec1=0.35,c2=0.2,μ1=0.1, and α=5. (a) Bifurcation diagrams of q1and q2 on varying the parameter μ2. (b) Largest Lyapunov exponent with respect to μ2. (c) Basin of attraction of the period 2-cycle at μ=0.648. (d) Basin of attraction of the periods 6-cycle and 12-cycle at μ2=0.721 and μ2=0.7245, respectively. (e) The phase plane of the six pieces chaotic attractors at μ2=0.725 and the basin of attraction of the chaotic behavior converted to a period 2-cycle at μ2=0.7257. (f) The basin of attraction of period 8-cycle at μ2=0.7259 and the chaotic behavior converted to two spirals at μ2=0.7269.
(a)(b)(c)(d)(e)(f)We end this part of numerical investigation by the following dynamic behaviors of the map given by Figures4(f) and 5(a). They give information about a stable period 8-cycle, a chaotic attractor that changes to two spiral points, and two closed invariant curves with chaotic attractor due to Neimark–Sacker bifurcation. Conversely, when c1<c2, a Neimark–Sacker bifurcation with respect to the parameter μ2 takes place. We give in Figure 5 some numerical simulations showing different behaviors of the map for different values of the bifurcation parameter.Figure 5
The fixed parameter values areq0,1,q0,2=0.40,0.11,c1=0.1,c2=0.2,μ1=0.3, and α=8. (a) Two closed invariant curves at μ2=0.73. The time series of q1 and q2 corresponding to these two invariant curves are also shown. (b) Bifurcation diagrams of q1and q2 on varying the parameter μ2. (c–f) Different behaviors of the map at different values of the bifurcation parameter μ2. Time series of q1 and q2 corresponding to the chaotic attractor are also shown.
(a)(b)(c)(d)(e)(f)


## 5. Noninvertible Map and Critical Curves

Due to the complex structure of the basin of attraction, we discuss here some important characteristics of map (10). By setting q1,t+1=q´1 and q2,t+1=q´2 in the two-dimensional map (10), the time evolution of the two competing firms is obtained by the iteration of the map T:q1,q2⟶q´1,q´2 given by(23)T:q´1=q1+αq11−μ1q2q1+q22−c1,q´2=1−μ2q1c2−q1.To understand the structure of basin and its qualitative changes, the inverse of the above map becomes important. Map (23) is a noninvertible map. This can be obtained from the fact that given a point q´1,q´2∈ℝ, its rank-1 preimages may be up to four; those four points are computed by solving the algebraic system given in (23) with respect to q1 and q2. The phase space can be divided into several regions denoted by Zi, where i refers to the number of rank-1 preimages corresponding to each region in the phase plane. Any two contiguous regions are different from each other by two real preimages. The critical curve LC which generalizes the critical value (local maximum or minimum) of one dimension to the two-dimensional framework is the locus which has two or more coincide preimages. Similarly, LC−1 denotes the set of critical points (local extreme points). LC−1 contains all the points at which the following Jacobian determinant of (23) vanishes:(24)α1−μ1q1q1−q22c2c21−μ2q1+μ2−1q1+q23c21−μ2q1=0.We should highlight that map (23) is not defined at q1,q2=0,0. This means that LC−1 becomes LC−1⊆q1,q2∈ℝ+2:q1=q2,q1≠0≠q2. The critical curve can be obtained easily from LC=TLC−1. One can see from Figure 6(b) that LC separates the phase plane into two regions: Z0 and Z2. The region Z0 has no preimages of rank-1, while the region Z2 has two preimages or rank-1. Even if we are unable to get analytical expressions for these preimages, we can still use numerical calculation to obtain those points. At the parameter values c1=0.35,c2=0.2,μ1=0.1,μ2=0.5, and α=5.8, we get only two rank-1 preimages’ points located in region Z2. One of which is located below LC−1 and the other point located above it, as shown in Figure 6(b). Therefore, map (23) is a Z0−Z1 noninvertible map. Figure 6(a) represents the phase plane corresponding to the same parameter values. The gray color refers to the divergent and unfeasible trajectories denoted by B∞. The other two colors indicate the basin of attraction of the Nash point and the basin of attraction of the period 2-cycle coexisting with the Nash point, denoted by BE and BC, respectively. In general, divergent trajectories along the q2−axis can be obtained starting from an initial condition out of this axis. This means that the preimages on this axis can be calculated as follows. Let the point q´1,0∈ℝ+2. Then, its preimages are the real solutions of the following algebraic system obtained by substituting this point in (23):(25)q´1=q1+αq11−μ1q2q1+q22−c1,0=1−μ2q1c2−q1.Figure 6
(a and b) Phase plane corresponding to the parameter valuesc1=0.35,c2=0.2,μ1=0.1,μ2=0.5, and α=5.8. The gray color refers to the divergent and infeasible trajectories, while the other two colors are for the basin of attraction of Nash point and the period 2-cycle coexisting with it.
(a)(b)From the second equation of (25), one can see that the preimages of that point are either located on the same invariant axis q1=0 or on the straight line q1=1−μ2/c2.Proposition 3.
The real preimages of the point0,q´2∈ℝ+2 belong to the same invariant axis q1=0 or lie on the curve:(26)q2=α1−μ12αc1−1−q1+α1−μ1αc1−1α1−μ14αc1−1−q1,q1≤α1−μ14αc1−1.Proof.
Substituting0,q´2 in (23), we obtain(27)0=q1+αq11−μ1q2q1+q22−c1,q´2=1−μ2q1c2−q1.
From the first equation of (27), the preimages of that point belong to the same invariant axis q1=0 or lie on the curve defined by the following equation:(28)0=1−αc1+α1−μ1q1+q22q2.
From (28), the term q1+q22=α1−μ1/αc1−1q2 must be positive under the condition αc1≥1. Solving (28) with respect to q2 completes the proof.Proposition 4.
Map (23) has two real preimages O−11 and O−12, where(29)O−11=0,α1−μ1αc1−1,O−12=1−μ2c2,α1−μ12αc1−1−q1+α1−μ1αc1−1α1−μ14αc1−1−q1,q1≤α1−μ14αc1−1.Proof.
The proof is straightforward and is obtained by substitutingq1=0 or q1=1−μ2/c2 in (27).
Therefore, the two preimagesO−11 and O−12 are corner points on the hyperplane that describes the basin of attraction. All the points outside this hyperplane cannot generate any feasible trajectories.


## 6. Conclusion and Future Work

This paper features complex behaviors of heterogeneous Cournot duopoly game in which players optimize their profits under minimal sales constraints. In cases of sales constraints and no sales constraints, the steady states points and the local stability conditions of the Nash equilibrium have been obtained. We have proven that the proposed game inherently produces two routes to chaos: the flip and Neimark–Sacker bifurcation, even though different values of the parameters of the minimum sales constraints are present in the firms. The plots of 1D bifurcation diagrams are used to further illustrate those two routes to chaos. We also use numerical simulations to describe the dynamic phenomenon of the game through 2D bifurcation diagram, 1D bifurcation diagram, Lyapunov maximum exponent, time series map, phase diagram, basins of attractions, and so forth.Our numerical results show the Nash equilibrium point is very sensitive to any change in the adjustment speed and minimal sales constraints parameters, the latter parameter being the main difference between this study and Tramontana [30]. Furthermore, the two routes to chaos give rise to complex dynamic behaviors when the changes in the parameters are large. Numerical simulations also show that the basin of attractions of the feasible and unfeasible trajectories of map (10) depend on the cost, adjustment speed, and minimum sales constraint parameters.Based on the result in this study, economical firms have better chance of reaching the Nash equilibrium by setting their minimum sales constraints to a certain value while varying their adjustment speed from time to time according to the stability condition. However, if this is not possible, another approach that can control their output towards a stable Nash is by employing a controlling parameter method such as delayed feedback control, which can be explored in future study.



---

*Source: 6687544-2021-03-28.xml*