Urban traffic resilience control - An ecological resilience perspective

Shengling Gao School of Mathematical Sciences, Beihang University, Beijing,
China<EMAIL_ADDRESS>
Zhikun She Corresponding author. School of Mathematical Sciences, Beihang
University, Beijing, China<EMAIL_ADDRESS>
Quanyi Liang Co-Corresponding author. School of Mathematical Sciences, Beihang
University, Beijing, China<EMAIL_ADDRESS>
Nan Zheng Co-Corresponding author. Institute of Transport Studies, Department
of Civil Engineering, Monash University, Clayton, Australia.
<EMAIL_ADDRESS>
Daqing Li School of Reliability and Systems Engineering, Beihang University,
Beijing, China<EMAIL_ADDRESS>
Urban traffic resilience has gained increased attention, with most studies
adopting an engineering perspective that assumes a single optimal equilibrium
and prioritizes local recovery. On the other hand, systems may possess
multiple metastable states, and ecological resilience is the ability to switch
between these states according to perturbations. Control strategies from these
two different resilience perspectives yield distinct outcomes. In fact,
ecological resilience oriented control has rarely been viewed in urban
traffic, despite the fact that traffic system is a complex system in highly
uncertain environment with possible multiple metastable states. This absence
highlights the necessity for urban traffic ecological resilience definition.
To bridge this gap, this paper defines urban traffic ecological resilience as
the ability to absorb uncertain perturbations by shifting to alternative
states. The goal is to generate a system with greater adaptability, without
necessarily returning to the original equilibrium. Our control framework
comprises three aspects: portraying the recoverable scopes; designing
alternative steady states; and controlling system to shift to alternative
steady states for adapting large disturbances. Among them, the recoverable
scopes are portrayed by attraction region; the alternative steady states are
set close to the optimal state and outside the attraction region of the
original equilibrium; the controller needs to ensure the local stability of
the alternative steady states, without changing the trajectories inside the
attraction region of the original equilibrium as much as possible. Note that,
this paper gives inner and outer estimations of attraction region with
explicit algebraic expressions, as the attraction region for nonlinear systems
are usually very complex and difficult to obtain. Comparisons with classical
engineering resilience oriented urban traffic resilience control schemes show
that, proposed ecological resilience oriented control schemes have better
adaptability and can generate greater resilience. These results will
contribute to the fundamental theory of future resilient intelligent
transportation system.

Urban traffic resilience; resilience control; ecological resilience;
macroscopic fundamental diagram; attraction region; stability; alternative
steady state

##  1 Introduction

Traffic congestion is plaguing most of megacities (Huang et al. (2020a), Zeng
et al. (2019)). Congested transportation systems are highly vulnerable, and
often fall into collapse under extreme events such as poor weather and
accidents. Such vulnerability of transportation network operations is a
"persistent urban disease". Therefore, how to sustain functional road
transportation systems has been a long-standing research topic, in particular
the challenging objective is to improve the resilience of operations towards
uncertain perturbations.

Resilience, first proposed by Holling (1973) in the ecology field, denotes the
system ability to absorb unexpected disturbances and to keep persistence. Over
the past few decades, resilience has been greatly developed in various
disciplines including ecology, engineering, and psychology (Martin-Breen and
Anderies (2011), Hosseini, Barker, and Ramirez-Marquez (2016a)). In general,
there are two perspectives on resilience: engineering resilience and
ecological resilience (Holling (1996)). Noteworthy, this categorization of
resilience is discipline-independent. One can use ball-and-cup heuristic to
illustrate these two perspectives (Scheffer et al. (1993), Rölfer, Celliers,
and Abson (2022), Liao (2012), Walker et al. (2004)), as shown in Fig. 1: when
the ball is at the cup bottom, it symbolizes a steady state, known as an
“attractor”; the shape of the cup portrays the recoverable phase space of the
corresponding attractor, usually described as the shape of potential function
(landscape) or attraction basin; the yellow arrows represent perturbations and
system responses. Fig. 1(a) illustrates the engineering resilience, which
usually assumes there exists a single steady state and characterizes the
system ability to absorb disturbances and recover to this single steady state.
Studies related to engineering resilience typically aim at rapid local or
global recovery of a single steady-state system (Tilman and Downing (1994),
O’Neill et al. (1986), Pimm (1984)). In contrast, Fig. 1(b) shows the
ecological resilience, which features multiple metastable states and the
perturbations are absorbed through a multi-state shifting, without requiring
the returning to a single steady state. Clearly, ecological resilience
emphasizes persistence with uncertainty (Holling (1973), Folke (2006), Zeng et
al. (2022)). As such, ecological resilience can transform a system under
instability from one basin of attraction to another (Holling (1973)), even for
cases where they are far from the stable equilibrium as defined otherwise in
engineering resilience.

(a) Engineering resilience

(b) Ecological resilience

Figure 1: Schematic diagram of (a) engineering resilience and (b) ecological
resilience (Liao (2012), Rölfer, Celliers, and Abson (2022))

Road network resilience is an emerging research subject. Most of the developed
works are towards engineering resilience, aiming at a rapid recovery to a non-
congested equilibrium or a single optimal equilibrium. For example, the
classic perimeter control schemes developed based on the macroscopic
fundamental diagram (short for MFD) characterization of congestion dynamics,
have received much attention in recent years (Geroliminis and Daganzo (2008),
Haddad and Geroliminis (2012), Aghamohammadi and Laval (2020), Laval,
Leclercq, and Chiabaut (2017), Johari et al. (2021), Haddad (2015), Khwais and
Haddad (2022), Haddad and Zheng (2020), Mariotte, Leclercq, and Laval (2017),
Haddad and Mirkin (2020), Ramezani, Haddad, and Geroliminis (2015), Li,
Yildirimoglu, and Ramezani (2021), Mohajerpoor, Cai, and Ramezani (2023),
Kouvelas, Saeedmanesh, and Geroliminis (2023), Tsitsokas, Kouvelas, and
Geroliminis (2023), Yang, Menendez, and Zheng (2019), Dantsuji, Fukuda, and
Zheng (2021), Ampountolas, Zheng, and Geroliminis (2017), Yang, Zheng, and
Menendez (2018), Su et al. (2020)). The perimeter control schemes regulate the
proportion of entering traffic at the perimeter of the networks via traffic
signal control. In this direction, Haddad and Geroliminis (2012) calculated
the "region of attraction" (short for RA) and designed a state feedback
control to increase the RA for enlarging the recoverable area towards a single
non-congested stable equilibrium. One pre-condition though, is that the system
is within the states as defined by RA. In face of the so-called “hyper-
congestion” which is far beyond the stable equilibrium state (e.g. a
congestion density value twice as much as the critical density) and beyond RA,
recovering to the desired system equilibrium cannot be guaranteed, and the
system will inevitably fail, e.g. becoming unrecoverable grid lock. On the
other hand, Zhong et al. (Zhong et al. (2020a, 2018a, 2018b)) and Huang et al.
(Huang et al. (2020b)) developed control schemes where the control would aim
at a designated point, for example the point near the critical density or
critical accumulation where the network flow reaches its maximum. However, the
sufficient conditions of these controllers require that the disturbance to the
systems, i.e. the incoming traffic demand to the road networks, cannot be high
and fast-varying over time, which limits the controllability towards real
world traffic situations. Most recently, Gao et al. (2022) proposed a novel
control scheme based on RA, denoted as RCS-single, focusing specifically
towards the recovery to the single optimal equilibrium from hyper-congested
states. This is one of the first exploratory works looking at the “resilience”
of system control. Nevertheless, this work was still within the framework of a
single equilibrium system, and the recovery speed is very slow. For systems
where instability is coupled with complex disturbances, in the case of hyper-
congestions accompanied by persistently high travel demand, the existing
controls are not only time-consuming, but also the objective of returning to
single equilibrium is impractical. In literature, the impact of large changes
in travel demand causing the instability of MFD was pointed out in multiple
studies, such as Daganzo, Gayah, and Gonzales (2011), Zhong et al. (2020b),
Haddad and Zheng (2020).

Urban traffic is a complex dynamic system, which may exist multiple metastable
states (Zeng et al. (2020)). Given the above discussion on the limitation of
the single-objective control schemes, a more realistic way appears to be
guiding the system to switch between multiple metastable states under
perturbations. While the engineering resilience-oriented system controls have
its values, there is a need to address the limitation towards system
perturbations. Future resilient traffic system, that is capable of
withstanding large and rapidly changing environments and contingencies,
absorbing environmental uncertainties and developing different adaptive
landscapes, would be by nature an ecosystem-like dynamic system. In other
words, controlling such system should embrace an ecological perspective, where
the ecological operation essentially enables the system to shift among
alternative stable states (Holling (1996), Folke et al. (2010), Walker and
Salt (2012), Gunderson (2000), Dakos and Kéfi (2022), Scoones et al. (2020)).
Therefore, this work develops an ecological resilience control for managing
congested urban road systems. The control aim is to enhance road traffic flow
systems with greater adaptability towards different perturbations, such as
congestion, unusual travel demand, sudden infrastructure closure, and
accidents. Our work will, for the first time, provide an analytical modeling
framework and the proves on system shifting between alternative steady states.
Our work will demonstrate that such system is superior and much more
resilient, comparing to controlled systems that aiming at returning to a
single state, such as an original equilibrium or a pre-defined optimal state.

Built upon the MFD-based traffic flow representation system, this paper
proposes a methodology for deriving ecological resilience control under large
disturbances. Our control methodology comprises three aspects: portraying the
recoverable system states, designing the target steady system states, and
developing control mechanisms to regulate system towards steady states.
Specifically, the RA framework will be adopted but extended, where the
recoverable scopes are portrayed by attraction region; the alternative steady
states are set close the optimal state and outside the attraction region of
the original equilibrium; the controller is devised to guarantee the local
stability of the alternative steady states, without changing the dynamics
around the original equilibrium as much as possible. We showcase the proposed
methodology in two-region road systems (i.e. two congested road networks whose
traffic flow and congestion interact with each other, and each road network
has its own MFD to describe its traffic dynamics). We first derive the
explicit algebraic expressions of the boundaries of inner and outer
estimations of the attraction regions. Building upon these theoretical
boundary delineations, we then design two distinct resilience control schemes,
denoted as RCS-1 and RCS-2, respectively. Considering a four-equilibria system
as an illustration, we test two resilient control schemes. The control
performance can be seen in Fig. 13 and Fig. 15 (more details and discussions
will be provided in section 5). CPC (denotes constant perimeter control) and
RCS-single (proposed by Gao et al. (2022)) are two benchmark methods. As a
general remark, the proposed resilience control schemes and RCS-single can
prevent the evolution of vehicle density (also called vehicle accumulation)
$n_{1}$ or $n_{2}$ to jam vehicle density (corresponding to zero completion
flow), while CPC is unable to impede such occurrences. Furthermore, the
resilience measures summarized in Table 2 indicate that, RCS-1 and RCS-2
outperform RCS-single as they exhibit greater resilience measures in most
cases, corresponding to smaller resilience triangle (resilience loss). This
superiority stems from the ability of RCS-1 and RCS-2 to guide the
trajectories towards alternative stable states that are more easily
attainable, thereby enhancing adaptability. Moreover, the local landscape
(local Lyapunov function (Wang et al. (2021), Willems (1971), Blanchini
(1995))) near the target point of controlled system with RCS-2 clearly shows
the effect of our control scheme.

The main contributions of this study encompass three key aspects:

  1. (1)

We establish a theoretical foundation for the approximate depiction of
recoverable states for multi-equilibrium systems, by deriving the inner and
outer estimations of attraction regions with explicit algebraic expressions
for multi-equilibrium MFD dynamics, utilizing invariant sets. Note that the
computation of RA for multiple equilibria nonlinear systems is especially
intricate and hard to obtain.

  2. (2)

We innovatively propose a leap from single-steady state control (from an
engineering resilience perspective) to multi-steady states control (from an
ecological resilience perspective), by leveraging the developed theory of
approximate depiction of recoverable states and the stability theory of switch
systems. This effort yields explicit control methodological, providing a clear
pathway for the implementation of a control system that navigates across
different stability landscapes.

  3. (3)

We innovatively design a resilience measure utilizing the classical resilience
triangle and completion flow. And it intuitively depicts the loss of completed
trips compared to the maximum potential number of trips that could have been
completed during the recovery period.

##  2 The control system dynamics and its ecological resilience

This paper showcases the advancement and performance of an ecological
resilience-oriented control in a two-region system (known also as two-
reservoir system in literature) traffic flow control problem, where the
traffic dynamics of the two regions follow an MFD-based representation. Sec.
2.1 firstly gives the mathematical description of the MFD-based modeling
framework, Sec. 2.2 illustrates the ecological resilience control framework
for MFD dynamics, and then Sec. 2.3 defines ecological resilience measure.

###  2.1 Two-region MFD dynamics

To start with, we consider a two-region MFD system, illustrated in Fig. 2. The
urban traffic network is partitioned into three parts: Region 1, Region 2, and
the Outside. Our primary focus lies on Regions 1 and 2, and each endowed with
its own MFD. The MFD delineates the correlation between the travel completion
flow $G_{i}(n_{i}(t))$ and vehicle accumulation (density) $n_{i}$ of region
$i$ at time $t$, manifesting a single peak, characterized by an initial ascent
succeeded by a descent. Classical perimeter control methodologies regulate the
traffic flow ratio entering Regions 2 and 1 through signal controls, denoted
as allowed pass rates $u_{1}$ and $u_{2}$. As illustrated in Fig. 2, for each
region at time $t$, there are two inflows (marked as gray arrows directed
toward region i) and one completion flow $G_{i}(n_{i}(t))$ (marked as gray
circle). The first inflow, ${u_{j}}{G_{j}}{(n_{j}(t))}$ ($j=1,2$, $j$ $\neq$
$i$) signifies the allowed traffic inflow transferred from $j$ to i with the
pass rate $u_{j}$. The second inflow, denoted as $d_{i}$, represents the fixed
net inflow from the outside to $i$, encapsulating uncontrollable demand,
commonly known as disturbance.

Figure 2: Two region MFD dynamics (Aboudolas and Geroliminis (2013), Gao et
al. (2022)).

The model can be described by MFD dynamic (Aboudolas and Geroliminis (2013),
Gao et al. (2022)):

$\begin{split}\frac{\mathrm{d}{n}_{1}(t)}{\mathrm{d}t}=F_{1}({n}_{1}(t),{n}_{2}(t)):=-{G_{1}}{(n_{1}(t))}+{u_{2}}{G_{2}}{(n_{2}(t))}+d_{1},\\\
\frac{\mathrm{d}{n}_{2}(t)}{\mathrm{d}t}=F_{2}({n}_{1}(t),{n}_{2}(t)):=-{G_{2}}{(n_{2}(t))}+{u_{1}}{G_{1}}{(n_{1}(t))}+d_{2},\end{split}$
(1)

with the boundary condition:

$\begin{split}&0\leq n_{1}(t)\leq p_{1},\,0\leq n_{2}(t)\leq p_{2},\\\ &0\leq
u_{1}\leq
u_{1,max},\,u_{1,max}=\min\left\\{1,\frac{G_{2,max}}{G_{1,max}}\right\\},\\\
&0\leq u_{2}\leq
u_{2,max},\,u_{2,max}=\min\left\\{1,\frac{G_{1,max}}{G_{2,max}}\right\\}.\end{split}$
(2)

Let $\bm{n}=(n_{1}(t),n_{2}(t))$,
$\bm{F}(\bm{n})=(\left(F_{1}(n_{1}(t),n_{2}(t),F_{2}(n_{1}(t),n_{2}(t))\right)$,
and $\dot{\bm{n}}=\frac{\mathrm{d}{\bm{n}}}{\mathrm{d}t}$, system (1) can be
simplified to $\dot{\bm{n}}=\bm{F}(\bm{n})$. Note that perimeter controllers
$u_{1}$ and $u_{2}$ are assumed to be constants. Therefore, they can be
denoted as constant perimeter control, short for CPC. The fixed net inflow
$d_{i}$ is also assumed to be constants. Moreover, as a realistic scene, we
here use the parabolic-MFD: $G_{i}(n_{i})=-{a_{i}}{n_{i}}({n_{i}-{p_{i}}})$
with opening size $a_{i}>0$ and jam accumulation $p_{i}$ of region $i$, as
shown in Fig. 2 and Fig. 3. Then we have the maximum capacity
$G_{i,max}=\frac{{a_{i}}{{p_{i}}^{2}}}{4}$ at the critical vehicle
accumulation $\frac{p_{i}}{2}$. Note that $d_{j}\leq{G_{j,max}}$,
$u_{i}{G_{i,max}}\leq{G_{j,max}}$ and $u_{1}\times u_{2}\neq 1$ are satisfied
to avoid overflow. Without loss of generality, we can assume $G_{1,max}\leq
G_{2,max}$, then the condition of $u_{1}$ and $u_{2}$ in Eq. (2) can be
simplified as Condition $(\mathbb{H})$:

$0\leq u_{1}\leq 1,\ 0\leq
u_{2}\leq\frac{G_{1,max}}{G_{2,max}},~{}\text{but}~{}u_{1}\times u_{2}\neq 1.$
(3)

For the two-region MFD dynamic (1) with boundary condition (2) under traffic
demand $(d_{1},d_{2})$, given the MFD function $(G_{1}(n_{1}),G_{2}(n_{2}))$,
we aim to generate a transportation system from ecological resilience
perspective by combining the stability characteristics of the MFD dynamics.

To achieve this objective, we first introduce concepts related to stability
characteristics. The stability characteristics can be elucidated by phase
portrait, where the phase portrait is a geometric representation of the
dynamic system trajectories on the phase plane. As shown in Fig. 3(a), the
phase plane is $n_{1}-n_{2}$, and each point (state) in $n_{1}-n_{2}$ plane
denotes the vehicle accumulation in two regions. Progressing to the right or
upward along the $n_{1}-n_{2}$ plane signifies a corresponding increase in the
vehicle accumulation within Region 1 (Region 2). Moreover, the axis $n_{2}$ in
Fig. 3(a) corresponds to the ordinate of the right MFD and axis $n_{1}$
corresponds to the abscissa of the upper MFD. The completion flow
$G_{1}(n_{1})$($G_{2}(n_{2})$) increases with the vehicle accumulation before
the critical vehicle accumulation $\frac{p_{1}}{2}$($\frac{p_{2}}{2}$),
reaching a peak $G_{1,max}$ ($G_{2,max}$), and subsequently decreases until it
diminishes to zero at the right (upper) boundary, signifying a traffic jam. As
shown in Fig. 3(a), the attraction region (light pink region) of the
equilibrium point $\bm{n}^{0}=(n^{0}_{1},n^{0}_{2})$ (green circle) is a set
of the initial points. The trajectories start from these initial points will
converge towards $n^{0}$, while trajectories start from outside attraction
region will move away from $n^{0}$ and towards either the upper or right
boundary. Note that equilibrium point $\bm{n}^{0}=(n^{0}_{1},n^{0}_{2})$
satisfies $F_{1}(n^{0}_{1},n^{0}_{2})=0,~{}F_{2}(n^{0}_{1},n^{0}_{2})=0$. We
can introduce the definition of attraction region for the general equilibrium
$\bm{n}^{\ast}=(n^{\ast}_{1},n^{\ast}_{2})$ as follows:

###### Definition 2.1

(Gao et al. (2022)) For system (1) with boundary restrictions (2), denoting
$n(t)=(n_{1}(t),n_{2}(t))$ and $\mathbb{D}=\left\\{(n_{1}(t),n_{2}(t))|0\leq
n_{i}(t)\leq p_{i},\,i=1,2\right\\},$ define

$\mathcal{R}(n^{\ast})=\left\\{n_{init}\in\mathbb{D}|n(0)=n_{init}\
and\lim_{t\to+\infty}n(t)=n^{\ast}\right\\},$ (4)

as the attraction region for $n^{\ast}=(n^{\ast}_{1},n^{\ast}_{2})$.

Note that, the classical concept of the "attraction domain" was originally
defined for asymptotically stable equilibria (Wang, She, and Ge (2020b), Zheng
et al. (2018), Ratschan and She (2010), Wang, She, and Ge (2020a)).
Especially, for MFD dynamic (1), Gao et al. (2022) defined the "attraction
region" for the unstable single equilibrium point. Differently, we here in
this paper define for the general equilibrium $(n^{\ast}_{1},n^{\ast}_{2})$,
including both unstable equilibrium and stable equilibrium.

Figure 3: Schematic representation of the phase portraits under (a) CPC, (b)
RCS-single, and (c) desired control. Note that desired phase portraits exhibit
multiple metastable states. The attraction region of each metastable state is
demarcated by different colors. The green area signifies the selectable range
of the alternative steady states. the green circle represents the steady
states; the blue lines with arrows depict state trajectories, indicating the
evolution of a point (state) over time. The phase portraits and the
Macroscopic Fundamental Diagram (MFD) exhibit a correlation. Each point in the
phase portraits corresponds to the vehicle density (accumulation) in two
regions, and the corresponding completion flows can be identified on the MFD.

###  2.2 Ecological resilience control framework for MFD dynamics

In the traffic flow theory and traffic control community, how to recover the
traffic system from congestion has been a challenging research question. The
classical perimeter control schemes, which utilizes MFD as references for
optimal control, target on recovering the traffic system to a specified
equilibrium. Such schemes are known to be constrained and uneconomical for
hyper-congestion conditions, because these conditions are far from the single
optimal control objective and there exist complex dynamics in the system that
prevent the system from even getting closer to the targeted equilibrium. For
this reason, our approach defines a new control objective, built upon which a
corresponding control methodology is developed.

In essence, ecological resilience considers the existence of multiple
metastable states in a dynamic system, and it enables the system to shift
between states as a way to absorb or adapt to rapid external changes. Inspired
by this nature, we aim to construct an ecological resilience traffic system
resembling Fig. 3(c), with multiple metastable states (marked as green
circles), where multiple metastable states includes the original uncongested
equilibrium $\bm{n}^{0}$ and alternative steady states
$\bm{n}^{i}=({n}^{i}_{1},{n}^{i}_{2})$ ($i=1,\cdots,r$ with number of
alternative steady states $r$).

To generate the ecological resilience traffic system resembling Fig. 3(c), our
first critical issue is how to select multiple metastable states. Initially,
we aim to leverage the inherent landscape of the system as much as possible.
Therefore, we preserve the original non-congested equilibrium point of the
system, denoted as $\bm{n}^{0}$. This point coincidentally aligns with the
single equilibrium point set by the classical CPC control system (as shown in
Fig. 3(a)) or the single equilibrium point of the RCS-single control system
proposed by Gao et al. (2022) (as depicted in Fig. 3(b)), satisfying
$F_{1}(\bm{n}^{0})=0,~{}F_{2}(\bm{n}^{0})=0$. After determining $\bm{n}^{0}$,
the next step involves selecting alternative stable states
$\bm{n}^{i}=({n}^{i}_{1},{n}^{i}_{2})$. There are two key criteria for
selecting alternative stable states: one to position them as close to the
optimal and critical vehicle accumulation
$\left(\frac{p_{1}}{2},\frac{p_{2}}{2}\right)$ as possible; the other is to
position them outside the attraction region $\mathcal{R}(\bm{n}^{0})$
(depicted in the light pink region) of the original non-congested equilibrium
point $\bm{n}^{0}$. The first criterion seeks to increase the completion flow,
given that such completion flow tends to decline as the vehicle accumulation
gets further from the critical vehicle accumulation (see Fig. 3(c)). The
second criterion is crucial because states within the attraction region of
$\bm{n}^{0}$ will spontaneously recover to $\bm{n}^{0}$. Here, we provide a
quantitative description of the selectable range. For the first criterion, we
aim to select $\bm{n}^{i}=(x,y)$, such that the absolute distance between
completion flow $(G_{1}(x),G_{2}(y))$ and maximum completion flow
$(G_{1,max},G_{2,max})$ is less than $\Delta_{0}$, where $\Delta_{0}$ denotes
the absolute distance between completion flow of the original equilibrium
$\bm{n}^{0}$ and maximum completion flow. That is,

$\left|G_{1}(x)-G_{1,max}\right|+\left|G_{2}(y)-G_{2,max}\right|\leq\Delta_{0},$
(5)

substituting $G_{i}(n_{i})=-{a_{i}}{n_{i}}({n_{i}-{p_{i}}})$ and
$G_{i,max}=\frac{{a_{i}}{{p_{i}}^{2}}}{4}$, yields

$\mathcal{E}=\left\\{(x,y)|\frac{\left(x-\frac{p_{1}}{2}\right)^{2}}{a_{2}}+\frac{\left(y-\frac{p_{2}}{2}\right)^{2}}{a_{1}}\leq\frac{\Delta_{0}}{a_{1}a_{2}}\right\\},$
(6)

$\mathcal{E}$ corresponds to the interior of an ellipse. Combining the second
criterion, the selectable range of alternative stable states $\mathcal{S}$ is

$\mathcal{S}=\left\\{(x,y)|(x,y)\in\mathcal{E}\bigwedge(x,y)\notin\mathcal{R}\right\\},$
(7)

as the green region shown in Fig. 3(c).

Once multiple metastable states have been chosen, the next consideration is
how to devise a control scheme that enables the system to recover to
appropriate steady states. Controlling a system with multi-steady states poses
more challenges compared to a system with a single stable state. The first
challenge lies in determining which state the system should recover to, and
the second challenge is devising a strategy to facilitate the system recovery
to the respective target stable state. To address the first challenge, we need
to identify the attraction region of the original equilibrium. It is essential
to note that, for nonlinear systems, the attraction region is typically
intricate, and obtaining explicit algebraic expressions for it can be
challenging. In this direction, Gao et al. (2022) derived the attraction
region with explicit algebraic expressions for the single-equilibrium MFD
dynamics. Different from their work, we consider multi-equilibria (including
two-equilibria and four-equilibria) systems in this paper. Given the
involvement of multiple equilibrium points and nontrivial boundary shapes,
their approach can not be applicable here. To this end, our first effort is
the estimation of the attraction region for the multi-equilibria systems (see
Sec. 3 and App. 7). Upon obtaining an estimation of the attraction region with
explicit expressions for the boundaries (including inner and outer
estimations) of the original uncongested equilibrium point, a natural boundary
between the original uncogested equilibrium and the alternative steady state
emerges.

The second challenge, i.e., how to devise a strategy to facilitate the system
recovery to the respective target steady state is also intricate. To tackle
this issue, we make a second effort using a switched controlled system and
design two control schemes based on the inner and outer estimations of the
attraction region (see Sec. 4). It is essential to emphasize that Gao et al.
(2022) also explored this direction. Their switched controlled system is
depicted in Fig. 3(b), where trajectories originating outside the attraction
region are directed towards its boundary before spontaneously converging
towards the single equilibrium. The proposed control, however, addresses a
significantly different and more complex problem. We aim to utilize the
escaped trajectories (from attraction region) and guide them to appropriate
steady points. Therefore, this novel framework faces new technical barriers,
namely how to find a feasible control that can recover system to $\varepsilon-
Neighborhood$ of target point $(\bar{n}_{1},\bar{n}_{2})$ within recovery time
$t_{f}$ under corresponding perturbations, where $\varepsilon$ denotes an
arbitrarily small positive number; and $\varepsilon-Neighborhood$ of a state
$(\bar{n}_{1},\bar{n}_{2})$ is described as the set
$\\{(n_{1},n_{2}):~{}(n_{1}-\bar{n}_{1})^{2}+(n_{2}-\bar{n}_{2})^{2}\leq\varepsilon\\}$.
Note that the target point $(\bar{n}_{1},\bar{n}_{2})$ can be any state in
multiple metastable states, depending on the size of the perturbation; and
system can switch between multiple metastable states in facing with different
perturbations. Moreover, the system can also switch between optimal and
suboptimal states, if the optimal state satisfies the conditions of the
alternative steady states.

###  2.3 Measuring urban traffic resilience

We have depicted the envisioned ecological traffic resilience system and
presented the resilience control framework. However, a quantitative indicator
for traffic ecological resilience is still absent. It is imperative to
establish a quantitative measure for traffic ecological resilience, thereby
enhancing a more intuitive perception of the ecological resilience of traffic
and facilitating a more straightforward evaluation and reference for control
effectiveness.

The quantitative assessment of resilience can be categorized into two main
types: general resilience measure methods and structure-based model-driven
resilience measure (Hosseini, Barker, and Ramirez-Marquez (2016b)). We
prioritize general resilience measure in this paper, given that the second
type is more reliant on specific models. Among general resilience measure, the
“resilience triangle” introduced in Bruneau et al. (2003) has gained notable
attention. Originally designed for assessing a community resilience loss
during earthquakes, it calculates the resilience loss through a definite
integral:

Figure 4: Schematic representation of resilience triangle. Here we consider
hyper-congestion (large $D_{s}(t_{0})$) scenario, and we assume that the
original equilibrium is a non-congested equilibrium (not the optimal state).
CPC is unable to recover the system, and instead leading to a rapid collapse.
Proposed control (RCS-single) recovers the system to a functional level with
deviation $-Ds_{1}$ ($-Ds_{2}$) from the maximum completion flow at
$t_{f}^{1}$ ($t_{f}^{2}$), having a smaller (larger) resilience loss, depicted
as the dark (light) brown shaded region. Note that the resilience loss is the
absolute value of the area enclosed by the resilience triangle. Proposed
control outperforms RCS-single, as proposed control can recover system to a
functional state with less deviation ($Ds_{1}<Ds_{2}$) faster
($t_{f}^{1}<t_{f}^{2}$).

$RL=\int_{t_{0}}^{t_{1}}\left(1-Q(t)\right)\mathrm{d}t,$ (8)

where $Q(t)$ represents community service quality at time $t$, expressed as a
percentage scale from 0 to 100$\%$. $Q(t)$ starts declining at $t_{0}$ and
returns to the normal state at $t_{1}$. The deviation between the quality of
degraded infrastructure and the normal infrastructure quality is quantified
through $1-Q(t)$. A higher RL indicates lower resilience, while a lower RL
suggests higher resilience.

In this direction, Gao et al. (2022) marks one of the initial attempts to
provide a quantitative interpretation of urban traffic resilience. They
defined urban traffic resilience as the integral of the absolute deviation
between vehicle density (vehicle accumulation) and the optimal vehicle density
during the recovery period. However, the defined resilience measure lacks an
evident physical interpretation. Thus, we revise the definition of deviation
$D_{s}(t)$ as the difference between the completion flow $G_{i}(n_{i})(i=1,2)$
and the maximum completion flow $G_{i,max}$. That is:

$D_{s}(t)=\left|G_{1}(n_{1})-G_{1,max}\right|+\left|G_{2}(n_{2})-G_{2,max}\right|.$
(9)

Note that, the unit of $D_{s}(t)$ is $[veh/h]$. Based on the defined
$D_{s}(t)$, we can define the ecological traffic resilience measure as
follows:

$R=\int_{t_{0}}^{t_{f}}\left(0-D_{s}(t)\right)\mathrm{d}t=-\int_{t_{0}}^{t_{f}}D_{s}(t)\,\mathrm{d}t.$
(10)

Here $t_{f}$ is the time when system recovers and stabilizes at functional
steady states. Formula (10) implies that $R$ is negative. Thus, $-R$
corresponds to the area of the resilience triangle, depicted as the dark brown
shaded region in 4. Noteworthily, the unit of resilience measure $R$ is
$[veh]$, and it denotes the loss of completed trips compared to the maximum
potential number of trips that could have been completed during the recovery
period.

We depict a schematic diagram of the resilience triangle under three different
controls, as shown in Fig. 4. Notably, under hyper-congestion, CPC fails to
recover system and instead rapidly collapses system; while RCS-single and the
proposed control can ensure recoverability. The resilience control proposed in
this paper outperforms RCS-single, as it enables the system to recover to an
alternative stable states with higher completion flow within a shorter time.
Thus, the resilience triangle (loss) under our proposed control method (the
dark brown area in Fig. 4) is smaller than that of RCS-single (the light brown
area). Note that the original equilibrium is assumed to be a stable non-
congested state, rather than an optimal state.

##  3 Stability analysis and spontaneous attraction region estimation for
multi-equilibria system

In the previous section, for multi-equilibria system, we defined ecological
resilience of urban traffic with an MFD-based representation. In this section,
we will derive the global phase portrait and the spontaneous attraction region
for the multi-equilibria system under CPC ($u_{1}$ and $u_{2}$). Sec. 3.1
discusses and verifies the local stability of the MFD dynamic; Sec. 3.2
derives and illustrates the global phase portrait; and Sec. 3.3 presents the
estimation of spontaneous attraction region that are used to design the
proposed control in Sec. 4. Note that multiple equilibria systems include both
two-equilibrium systems and four-equilibrium systems. According to Gao et al.
(2022), two equilibria exist if and only if Condition $(\mathbb{K}^{2a})$:
$\frac{{u_{2}}{{d_{2}}}+d_{1}}{1-{u_{1}}{u_{2}}}=G_{1,max}$ $\wedge$
$\frac{{u_{1}}{{d_{1}}}+d_{2}}{1-{u_{1}}{u_{2}}}<G_{2,max}$ or Condition
$(\mathbb{K}^{2b})$:
$\frac{{u_{2}}{{d_{2}}}+d_{1}}{1-{u_{1}}{u_{2}}}<G_{1,max}$ $\wedge$
$\frac{{u_{1}}{{d_{1}}}+d_{2}}{1-{u_{1}}{u_{2}}}=G_{2,max}$ holds; four
equilibria exist if and only if Condition $(\mathbb{K}^{4})$:
$\frac{{u_{2}}{{d_{2}}}+d_{1}}{1-{u_{1}}{u_{2}}}<G_{1,max}$ $\wedge$
$\frac{{u_{1}}{{d_{1}}}+d_{2}}{1-{u_{1}}{u_{2}}}<G_{2,max}$ holds. We in this
section mainly present the theoretical analysis for two-equilibria system. The
theoretical analysis for four-equilibria system can be found in App 7.

###  3.1 Local stability verification

Consider the two-equilibria case with Condition $(\mathbb{K}^{2a})$. According
to Condition $(\mathbb{K}^{2a})$, we have:

$\begin{split}&d_{1}=\left(1-u_{1}u_{2}\right)G_{1,max}-{u_{2}}d_{2},\\\
&d_{2}<G_{2,max}-{u_{1}}G_{1,max}.\end{split}$ (11)

Bringing the first equation of (11) into system (1), the system dynamics
become:

$\displaystyle\frac{\mathrm{d}{n}_{1}(t)}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle-{u_{2}}{{a_{2}}{(n_{2}(t)-\frac{p_{2}}{2})^{2}}}+{{a_{1}}{(n_{1}(t)-\frac{p_{1}}{2})^{2}}}+u_{2}\left(\frac{a_{2}p^{2}_{2}}{4}-u_{1}\frac{a_{1}p^{2}_{1}}{4}-d_{2}\right),$
$\displaystyle\frac{\mathrm{d}{n}_{2}(t)}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle-{u_{1}}{{a_{1}}{(n_{1}(t)-\frac{p_{1}}{2})^{2}}}+{{a_{2}}{(n_{2}(t)-\frac{p_{2}}{2})^{2}}}+\left(\frac{u_{1}a_{1}p^{2}_{1}}{4}-\frac{a_{2}p^{2}_{2}}{4}+d_{2}\right).$
(12)

For systems (3.1), the two equilibria can be denoted as
$P^{2a}_{k}=(\frac{p_{1}}{2},p^{2a}_{k})$ ($k=1,2$), where
$p^{2a}_{k}=\frac{{p_{2}}+(-1)^{k}\times M}{2}$,
$M=\sqrt{{p_{2}}^{2}-\frac{4{d_{2}}+{u_{1}}a_{1}p^{2}_{1}}{{a_{2}}}}$. Upon
substituting Formula (11) into Formula (1), the system (1) reduces by one
variable, $d_{1}$, significantly simplifying our analysis. However, Formula
(3.1) remains cumbersome. To further simplify the system for subsequent
analysis, we shift equilibrium point $P^{2a}_{1}$ of the system (3.1) to
$(0,0)$. That is, letting $x_{1}(t)=n_{1}(t)-\frac{p_{1}}{2}$ and
$x_{2}(t)=n_{2}(t)-p^{2a}_{1}$, we can simplify system (3.1) as:

$\displaystyle\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle
H^{2a}_{1}({x_{1}}(t),{x_{2}}(t)):={u_{2}}{a_{2}}M{x_{2}}(t)+\left({{a_{1}}{x_{1}}^{2}(t)}-{u_{2}}{{a_{2}}{x_{2}}^{2}(t)}\right),$
$\displaystyle\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle
H^{2a}_{2}({x_{2}}(t),{x_{2}}(t)):=-{a_{2}}M{x_{2}}(t)+\left(-{u_{1}}{{a_{1}}{x_{1}}^{2}(t)}+{{a_{2}}{x_{2}}^{2}(t)}\right),$
(13)

The two equilibria for the new system (3.1) are $\hat{P}^{2a}_{1}=(0,0)$ and
$\hat{P}^{2a}_{2}=(0,M)$. Noteworthy, system (3.1) possesses the same
properties as system (3.1) since system (3.1) is obtained by shifting system
(3.1). Due to the simpler formulation of system (3.1), we in the following use
system (3.1) to derive local stability characteristics, global phase
portraits, and attraction region estimations. The properties obtained are
equally valid for system (3.1). Therefore, all our propositions and theorems
are summarized for system (3.1).

From the aforementioned conversion, we derive the following proposition for
its two equilibria.

###### Proposition 3.1

Under Condition $(\mathbb{H})$, the two equilibria $P^{2a}_{1}$ and
$P^{2a}_{2}$ of system (3.1) are both saddle-node points.

Proof: We have that $P^{2a}_{1}$ and $P^{2a}_{2}$ of system (3.1) are
equivalent to $\hat{P}^{2a}_{1}$ and $\hat{P}^{2a}_{2}$ of system (3.1), thus
we consider the stability of $\hat{P}^{2a}_{1}$ and $\hat{P}^{2a}_{2}$. For
$\hat{P}^{2a}_{1}$, we have $\left(\begin{array}[]{ccccc}\frac{\partial
H^{2a}_{1}(x_{1},x_{2})}{\partial x_{1}}&\frac{\partial
H^{2a}_{1}(x_{1},x_{2})}{\partial x_{2}}\\\ \frac{\partial
H^{2a}_{2}(x_{1},x_{2})}{\partial x_{1}}&\frac{\partial
H^{2a}_{2}(x_{1},x_{2})}{\partial x_{2}}\\\
\end{array}\right)\Bigg{|}_{(0,0)}=\left(\begin{array}[]{ccccc}0&{u_{2}}{a_{2}}M\\\
0&-{a_{2}}M\\\ \end{array}\right)$, indicating that:

$Rank\left\\{\left(\begin{array}[]{ccccc}\frac{\partial
H^{2a}_{1}(x_{1},x_{2})}{\partial x_{1}}&\frac{\partial
H^{2a}_{1}(x_{1},x_{2})}{\partial x_{2}}\\\ \frac{\partial
H^{2a}_{2}(x_{1},x_{2})}{\partial x_{1}}&\frac{\partial
H^{2a}_{2}(x_{1},x_{2})}{\partial x_{2}}\\\
\end{array}\right)\right\\}\Bigg{|}_{(0,0)}=1,$

and

$div(H^{2a}_{1},H^{2a}_{2})|(0,0)=\frac{\partial H^{2a}_{1}(0,0)}{\partial
x_{1}}+\frac{\partial H^{2a}_{2}(0,0)}{\partial x_{2}}\neq 0.$

According to Theorem 7.1 in Zhang (2006), we can analyze the nature of
equilibrium points $\hat{P}^{2a}_{1}$ by calculating the bifurcation function
$B(x_{1}(t))$ (defined in Zhang (2006)), yields:

$B(x_{1}(t))=H^{2a}_{1}(x_{1},x_{2}(x_{1}))=H^{2a}_{1}(x_{1},\pm\sqrt{\frac{u_{1}a_{1}x^{2}_{1}}{a_{2}}+\frac{M^{2}}{4}}+\frac{M}{2})=(1-u_{1}u_{2})x^{2}_{1}(t),$

where
$x_{2}(x_{1})=\pm\sqrt{\frac{u_{1}a_{1}x^{2}_{1}}{a_{2}}+\frac{M^{2}}{4}}+\frac{M}{2}$
is obtained by letting $H^{2a}_{2}(x_{1},x_{2})=0$. Since Condition
$(\mathbb{H})$ holds, we have $B(x_{1}(t))>0$ holds for all
$x_{1}(t)\in\left[-\frac{p_{1}}{2}~{}\frac{p_{1}}{2}\right]$ except
$x_{1}(t)=0$. Moreover, the power exponent of function $B(x_{1}(t))$ is 2.
Thus, $(0,0)$ is a saddle-node point (see Zhang (2006) Theorem 7.1(iii)),
indicating that $P^{2a}_{1}$ is a saddle-node point. Similarly, we can prove
that $P^{2a}_{2}$ is a saddle-node point. $\hfill\square$

Proposition 3.1 reveals that, for system (1) under $(\mathbb{K}^{2a})$, the
two equilibria $P^{2a}_{1}$ and $P^{2a}_{2}$ are both saddle-node point.
Similarly, for system (1) under $(\mathbb{K}^{2b})$, we can also prove that
two equilibria $P^{2b}_{1}$ and $P^{2b}_{2}$ are both saddle-node point. These
results will help us to reveal the local stability of the two equilibria,
which will assist in subsequent derivations of attraction region, providing
foundation for understanding the global properties of the system.

For the four-equilibria system, we also try to capture the local stability
through linearization, and the detailed analysis can be found in App. 7.1. In
App. 7.1, we denote the four equilibria as $P^{4}_{m}$ ($m=1,2,3,4$), and we
have that $P^{4}_{1}$ is a locally stable node, $P^{4}_{2}$ and $P^{4}_{3}$
are saddle points, $P^{4}_{4}$ is an unstable node, which can be summarized as
Proposition 7.1 in App. 7.1.

###  3.2 Global phase portrait derivation

To capture more properties of the system, we aim to delve further into its
global nature. Global analysis of dynamical systems is often intricate,
lacking a unified analytical approach. For MFD dynamics, we overcome this gap
by first dividing the phase space into several subregions using demarcation
lines, followed by delineating the trajectory trends within each subregion,
and ultimately analyzing and deriving the global phase portrait leveraging the
trajectory trends in each subregion. For simplicity in analysis, we continue
to consider the more straightforward formulation of Formula (3.1). Note that,
the aforementioned demarcation lines are obtained by solving
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$=0 or
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$=0 respectively. The trajectory
trends in each region are determined by
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$ and
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$. For instance, if
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$ > 0
($\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$> 0), the trajectories move to the
right (upward), while if $\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$< 0
($\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$< 0), the trajectories move to the
left (downward). Furthermore, considering that the demarcation lines
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$=0 or
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$=0 take different shapes under four
distinct parameter conditions, we will separately examine these under four
conditions.

First consider Condition $(\mathbb{H}_{1})$: $0<u_{1},~{}u_{2}\leq 1\wedge
u_{1}u_{2}\neq 1$. Under this condition, both demarcation lines manifest as
hyperbolas, illustrated by the red and blue dotted curves in Fig. 5(a), with
the expressions
$\frac{\left(x_{2}-\frac{M}{2}\right)^{2}}{{\frac{M^{2}}{4}}}-\frac{x^{2}_{1}}{\frac{u_{2}a_{2}M^{2}}{4a_{1}}}=1$
and
$\frac{\left(x_{2}+\frac{M}{2}\right)^{2}}{{\frac{M^{2}}{4}}}-\frac{x^{2}_{1}}{\frac{a_{2}M^{2}}{4u_{1}a_{1}}}=1$.
In addition, we have $\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}=0$ and
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}=\left(1-u_{1}u_{2}\right){a_{2}}{x_{2}(t)}\left({x_{2}(t)-M}\right)>0$
on the red dotted curves except for $(0,0)$ and $(0,M)$; and we have
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}=0$ and
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}=\frac{1-u_{1}u_{2}}{u_{1}}{a_{2}}{x_{2}(t)}\left({x_{2}(t)-M}\right)>0$
on the blue dotted curves except for $(0,0)$ and $(0,M)$. It is worth
mentioning that, since
$\frac{u_{2}a_{2}M^{2}}{4a_{1}}<\frac{a_{2}M^{2}}{4u_{1}a_{1}}$ holds, the
location of these four curves are is uniquely determined, as shown in Fig.
5(a). Then, the plane $\mathbb{R}^{2}$ can be divided into seven unbounded
regions by these four curves. Moreover, Table 1 summarizes the symbols of
$\frac{\mathrm{d}{x}_{i}(t)}{\mathrm{d}t}$ ($i=1,2$) in these seven regions.
On the basis of Table 1, the phase portrait of (3.1) under Condition
$(\mathbb{H}_{1})$ can be obtained, as shown in Fig. 5(a).

(a) $\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$

(b) $\mathbb{K}^{2a}\wedge\mathbb{H}_{2}$

(c) $\mathbb{K}^{2a}\wedge\mathbb{H}_{3}$

(d) $\mathbb{K}^{2a}\wedge\mathbb{H}_{4}$

Figure 5: Theoretical schematic phase portraits of the two equilibria for system (3.1) under Condition $(\mathbb{K}^{2a}\wedge\mathbb{H})$. Red circle: undesired saddle-node equilibrium. Blue circle: desired saddle-node equilibrium. Red dotted curves: $\dot{x}_{1}=0$, blue dotted curves: $\dot{x}_{2}=0$, green arrows: trajectory directions.  Table 1: Symbols of the seven regions for system (3.1) under Condition $(\mathbb{H}_{1})$. Derivative | Symbol in region
---|---
A | B | C | D | E | F | G
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$ | $+$ | $-$ | $+$ | $+$ | $+$ | $-$ | $+$
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$ | $+$ | $+$ | $+$ | $-$ | $+$ | $+$ | $+$

Subsequently, we can verify that regions B and E are both positive invariant
sets by utilizing proof by contradiction. Given that the trajectories starting
from B and E will not escape, we can further obtain that there exist no close
orbits for system (3.1) under Condition $(\mathbb{H}_{1})$.

When Condition $(\mathbb{H}_{2})$: $u_{1}=0\wedge 0<u_{2}\leq 1$,
$(\mathbb{H}_{3})$: $0<u_{1}\leq 1\wedge u_{2}=0$ or $(\mathbb{H}_{4})$:
$u_{1}=0\wedge u_{2}=0$ hold, the two demarcation lines can manifest as one
hyperbola and one straight line (as shown in Fig. 5(b) and 5(c)) or as two
straight lines (as shown in Fig. 5(d)). After obtaining the demarcation lines,
we follow similar steps as above. As such, we verify that no close orbits
exist under these conditions, and the corresponding phase portraits are
illustrated in Fig. 5(b), 5(c) and 5(d). Similarly, we can also obtain the
phase portrait for system (1) under Condition $(\mathbb{K}^{2b})$, and there
also exists no close orbit.

For the four-equilibria system, we also try to capture the global phase
portraits, and the detailed analysis can be found in App. 7.2. Note that the
steps of global phase portrait derivation differ from that used for two-
equilibrium system. Due to the increased number of equilibria, the combination
of demarcation lines becomes intricately complex. On one hand, variations in
demarcation line types (e.g., two hyperbolas, one hyperbola and one straight
line, two straight lines, etc.) lead to different phase portraits. On the
other hand, different relative positions of demarcation lines (e.g., whether
the foci of two hyperbolas lie on the same axis) add another layer of
complexity. The method for two equilibria are insufficient to address the
challenges posed by four equilibria. Furthermore, the increased number of
equilibrium points also hinders us from directly determining the locations of
corresponding separatrices. New methods are required to capture the global
phase portrait. We solve these in App. 7.2 by subdividing Condition
$(\mathbb{K}^{4})$ into five new distinct conditions. That is, Condition
$(\hat{\mathbb{K}}^{4}_{1})$: $M_{1}<0\wedge M_{2}<0$; Condition
$(\hat{\mathbb{K}}^{4}_{2})$: $M_{1}=0\wedge M_{2}<0$; Condition
$(\hat{\mathbb{K}}^{4}_{3})$: $M_{1}<0\wedge M_{2}=0$; Condition
$(\hat{\mathbb{K}}^{4}_{4})$: $M_{1}>0\wedge
M_{2}<\min\left\\{{-\frac{M_{1}}{u_{2}}},{-u_{1}M_{1}}\right\\}$; Condition
$(\hat{\mathbb{K}}^{4}_{5})$:
$M_{1}<\min\left\\{{-u_{2}M_{2}},{-\frac{M_{2}}{u_{1}}}\right\\}\wedge
M_{2}>0$, where
$M_{1}=d_{1}+u_{2}\frac{{a_{2}p_{2}}^{2}}{4}-\frac{{a_{1}p_{1}}^{2}}{4}$; and
$M_{2}=d_{2}+u_{1}\frac{{a_{1}p_{1}}^{2}}{4}-\frac{{a_{2}p_{2}}^{2}}{4}$. The
five conditions mentioned above correspond to different positions of
demarcation lines. In App. 7.2, by discussing the shape of critical
demarcation lines under each condition, we have summarized 20 distinct phase
portraits for four-equilibria system and verifies that no closed orbit exists.

###  3.3 Spontaneous attraction region estimation

In the following, for two-equilibria system (3.1), the attraction region for
the unstable saddle-node point $P^{2a}_{1}$ will be derived. Note that, Gao et
al. (2022) developed a method to find certain special trajectories, e.g.,
$x_{2}(t)=qx_{1}(t)$ for system with single-equilibrium. This form of special
solution cannot apply here for multi-equilibrium systems. In this case, the
theoretical attraction region boundaries are identified by finding certain
separatrices of the equivalent systems (3.1), and the inner and outer
attraction regions can be obtained by finding the positive invariant sets,
under the aforementioned four conditions respectively.

Firstly, we consider Condition $(\mathbb{H}_{1})$. As shown in Fig. 6(a), the
two red dotted demarcation curves can be denoted as
$x_{2}=\hat{r}^{2a}_{i}(x_{1})$ ($i=1,2$), where
$\hat{r}^{2a}_{i}(x_{1})=(-1)^{i}\sqrt{\frac{a_{1}}{u_{2}a_{2}}x_{1}^{2}+\frac{M^{2}}{4}}+\frac{M}{2}$.
Similarly, the two blue dotted demarcation lines can be denoted as
$x_{2}=\hat{b}^{2a}_{i}(x_{1})$ ($i=1,2$), where
$\hat{b}^{2a}_{i}(x_{1})=(-1)^{i}\sqrt{\frac{u_{1}a_{1}}{a_{2}}x_{1}^{2}+\frac{M^{2}}{4}}+\frac{M}{2}$.
Obviously, the vertexes of the demarcation lines are exactly the two
equilibria $\hat{P}^{2a}_{1}=(0,0)$ and $\hat{P}^{2a}_{2}=(0,M)$.

Assuming the lower blue dotted line $x_{2}=\hat{b}^{2a}_{1}(x_{1})$ intersects
$x_{1}=-k_{1}$ with $k_{1}>0$ at $(-k_{1},x^{1}_{D})$ and
$(-k_{1},x^{2}_{D})$, where $x^{1}_{D}<x^{2}_{D}$. The trajectory starting
from $(-k_{1},x^{1}_{D})$ and $(-k_{1},0)$ will go into $\hat{P}^{2a}_{1}$ in
E as $t\rightarrow+\infty$; and the trajectory starting from
$(-k_{1},x^{2}_{D})$ will escape A from the upper boundary of A and stay in B
moving away from $\hat{P}^{2a}_{1}$ as $t\rightarrow+\infty$. Thus, due to the
continuous dependence theorem of solutions on initial conditions (Zhang
(2006)), there exists a $D^{\ast}$ with $0<D^{\ast}<x^{2}_{D}$, such that the
trajectory starting from any $(-k_{1},x_{D})$ with $0<x_{D}\leq D^{\ast}$ will
enter E and then go to $\hat{P}^{2a}_{1}$ in E as $t\rightarrow+\infty$, and
the trajectory starting from any $(-k_{1},x_{D})$ with $D^{\ast}<x_{D}\leq
x^{2}_{D}$ will enter A, then escape A from the upper boundary of A and stay
in B moving away from $\hat{P}^{2a}_{1}$ as $t\rightarrow+\infty$. Thus, the
trajectory starting from $(-k_{1},x_{D^{\ast}})$, donated as
$\hat{\Phi}^{1}(\hat{P}^{2a}_{1})$=
$\hat{\Phi}^{2a,1}_{t}(-k_{1},x_{D^{\ast}})$, is a stable separatrix of
$\hat{P}^{2a}_{1}$, as shown in Fig. 6(a). Similarly, assuming the lower red
dotted line $x_{2}=\hat{r}^{2a}_{1}(x_{1})$ intersect $x_{2}=-k_{2}$ with
$k_{2}>0$ at $(x^{1}_{F},-k_{2})$ and $(x^{2}_{F},-k_{2})$. By similar
analysis, we can obtain that there exists an $F^{\ast}$ with
$0<F^{\ast}<x^{2}_{F}$, such that the trajectory starting from any
$(x_{F},-k_{2})$ with $0<x_{F}\leq F^{\ast}$ will enter E and then go to
$\hat{P}^{2a}_{1}$ in E as $t\rightarrow+\infty$, and the trajectory starting
from any $(x_{F},-k_{2})$ with $F^{\ast}<x_{F}\leq x^{2}_{F}$ will enter G,
then escape G from the upper boundary of G and stay in D moving away from
$\hat{P}^{2a}_{1}$ as $t\rightarrow+\infty$. Thus, the trajectory starting
from $(x_{F^{\ast}},-k_{2})$, donated as
$\hat{\Phi}^{2}(\hat{P}^{2a}_{1})=\hat{\Phi}^{2a,2}_{t}(x_{F^{\ast}},-k_{2})$,
is a stable separatrix of $\hat{P}^{2a}_{1}$, as shown in Fig. 6(a).

Subsequently, combining with phase portrait Fig. 5(a), we can verify that the
region $\hat{S}^{2a}_{in}$ is a positive invariant set utilizing proof by
contradiction, where $\hat{S}^{2a}_{in}=\left\\{(x_{1},x_{2})|x_{i}\leq 0,\
i=1,2\right\\}$, depicted as the light blue zone in Fig. 6(a). Or else, there
will be at least one trajectory originated from $\hat{S}^{2a}_{in}$, and enter
$\hat{V}^{2a}_{1}$ (or $\hat{V}^{2a}_{2}$) through the upper boundary (or
right boundary) of $\hat{S}^{2a}_{in}$, where
$\hat{V}^{2a}_{1}=\left\\{(x_{1},x_{2})|x_{1}<0,~{}0<x_{2}<\hat{b}^{2a}_{2}(x_{1})\right\\}$
and
$\hat{V}^{2a}_{2}=\left\\{(x_{1},x_{2})|x_{1}>0,~{}x_{2}<\hat{r}^{2a}_{1}(x_{1})\right\\}$.
Obviously, it is contradictory with the trajectory direction in phase portrait
5(a). Further, joining Fig. 6(a) with Fig. 5(a), we have that: the trajectory
starting from any point in $D\wedge\hat{S}^{2a}_{in}$ will escape D from the
lower boundary of D and entering E, and the trajectory starting from any point
in $F\wedge\hat{S}^{2a}_{in}$ will escape F from the left boundary of F and
entering E. Since E is a positive invariant set, and we have $\dot{x}_{1}>0$
and $\dot{x}_{2}>0$ in E, the trajectory starting from any point in E will go
to $\hat{P}^{2a}_{1}$ as $t\rightarrow+\infty$. Thus, the trajectory starting
from any point in the positive invariant set $\hat{S}^{2a}_{in}$ will go to
$\hat{P}^{2a}_{1}$ as $t\rightarrow+\infty$, implying
$\hat{S}^{2a}_{in}\subset\mathcal{R}(\hat{P}^{2a}_{1})$, i.e.,
$\hat{S}^{2a}_{in}$ is an inner estimation of attraction region (Zheng et al.
(2018), Wang, She, and Ge (2020b)) for $\hat{P}^{2a}_{1}$.

Similarly, region
$\hat{U}^{2a}=\hat{U}^{2a}_{1}\vee\hat{U}^{2a}_{2}\vee\hat{U}^{2a}_{3}$ (the
light yellow region shown in Fig. 6(a)) is also a positive invariant set,
where

$\displaystyle\hat{U}^{2a}_{1}$ $\displaystyle=$
$\displaystyle\left\\{(x_{1},x_{2})|\left(x_{1}<0\wedge
x_{2}\geq\hat{b}^{2a}_{2}(x_{1})\right)\vee\left(x_{1}>0\wedge
x_{2}\geq\hat{r}^{2a}_{2}(x_{1})\right)\right\\},$
$\displaystyle\hat{U}^{2a}_{2}$ $\displaystyle=$
$\displaystyle\left\\{(x_{1},x_{2})|x_{1}>0\wedge\hat{b}^{2a}_{2}(x_{1})\leq
x_{2}\leq\hat{r}^{2a}_{2}(x_{1})\right\\},$ $\displaystyle\hat{U}^{2a}_{3}$
$\displaystyle=$ $\displaystyle\left\\{(x_{1},x_{2})|x_{1}\geq
0\wedge\hat{r}^{2a}_{1}(x_{1})\leq
x_{2}<\hat{b}^{2a}_{2}(x_{1})\right\\}/(0,0).$

Note that $\hat{U}^{2a}_{1}$ and $\hat{U}^{2a}_{3}$ are also positive
invariant sets. In addition, the trajectory starting from point in
$\hat{U}^{2a}_{2}$ will escape from the upper boundary of $\hat{U}^{2a}_{2}$
and enter $\hat{U}^{2a}_{1}$ moving far from $\hat{P}^{2a}_{1}$, or escape
from the lower boundary of $\hat{U}^{2a}_{2}$ and enter $\hat{U}^{2a}_{3}$
moving far from $\hat{P}^{2a}_{1}$. Thus, we can obtain that
$\mathcal{R}(\hat{P}^{2a}_{1})\subset\hat{S}^{2a}_{out}=\mathbb{R}^{2}\setminus\hat{U}^{2a}$,
i.e., the region $\hat{S}^{2a}_{out}$ is an outer estimation of attraction
region for $\hat{P}^{2a}_{1}$. Further, we can obtain that
$\hat{\Phi}^{1}(\hat{P}^{2a}_{1})\in\hat{V}^{2a}_{1}$ and
$\hat{\Phi}^{2}(\hat{P}^{2a}_{1})\in\hat{V}^{2a}_{2}$, as the yellow lines
shown in Fig. 6(a).

(a) $\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$

(b) $\mathbb{K}^{2a}\wedge\mathbb{H}_{2}$

(c) $\mathbb{K}^{2a}\wedge\mathbb{H}_{3}$

(d) $\mathbb{K}^{2a}\wedge\mathbb{H}_{4}$

Figure 6: Certain separatrices of $\hat{P}^{2a}_{1}$=$(0,0)$ or
$\hat{P}^{2a}_{2}$=$(0,M)$, and the inner and outer estimations of attraction
regions for system (3.1). Red dotted curves $\dot{x}_{1}=0$; blue dotted
curves: $\dot{x}_{2}=0$; yellow curves: separatrices, also serving as
boundaries of attraction region; purple lines: boundary of the inner
estimation of attraction region; bright blue lines: boundary of the outer
estimation of attraction region. The green dashed grid is to distinguish
$\hat{U}^{2a}_{1}$, $\hat{U}^{2a}_{2}$ and $\hat{U}^{2a}_{3}$.

When Condition $(\mathbb{H}_{2})$ hold, we can obtain two stable separatrices
$\hat{\Phi}^{1}(\hat{P}^{2a}_{1})$ and $\hat{\Phi}^{2}(\hat{P}^{2a}_{1})$ of
$\hat{P}^{2a}_{1}$; the inner estimation of attraction region
$\hat{S}^{2a}_{in}$; and outer estimation of attraction region
$\hat{S}^{2a}_{out}$ by similar analysis. Furthermore, we have
$\hat{\Phi}^{1}(\hat{P}^{2a}_{1})\in\hat{V}^{2a}_{1}$ and
$\hat{\Phi}^{2}(\hat{P}^{2a}_{1})\in\hat{V}^{2a}_{2}$, as shown in Figs. 6(b).
When $(\mathbb{H}_{3})$ or $(\mathbb{H}_{4})$ hold, we can first simplify the
system (3.1) by bringing in corresponding conditions. Then, by repeating the
above analysis, we can obtain corresponding results, as shown in Figs. 6(c)
and 6(d).

Consequently, for the system (3.1), the estimation of the attraction region
for the saddle-node points $\hat{P}^{2a}_{1}$ has been obtained. Further,
letting $n_{1}(t)=x_{1}(t)+\frac{p_{1}}{2}$ and
$n_{2}(t)=x_{2}(t)+\frac{p_{2}-M}{2}$, the estimation of attraction region for
$P^{2a}_{1}$ can be obtained. Denoting
$b^{2a}_{i}(n_{1})=\hat{b}^{2a}_{i}(n_{1}-\frac{p_{1}}{2})+\frac{p_{2}-M}{2}$
and
$r^{2a}_{i}(n_{1})=\hat{r}^{2a}_{i}(n_{1}-\frac{p_{1}}{2})+\frac{p_{2}-M}{2}$,
where $i=1,2$, we can obtain the theorem as below.

###### Theorem 3.2

For the unstable saddle-node point $P^{2a}_{1}$, under Condition
$(\mathbb{H})$, we have $S^{2a}_{in}\subseteq
S^{2a}_{A}=\mathcal{R}(P^{2a}_{1})\subseteq S^{2a}_{out}$, where the
attraction region $S^{2a}_{A}$, the inner estimation of attraction region
$S^{2a}_{in}$ and the outer estimation of attraction region $S^{2a}_{out}$ are
defined as follows:

  1. (1)

Under Conditions $(\mathbb{H}_{1})$ and $(\mathbb{H}_{2})$,
$S^{2a}_{A}=\left\\{(n_{1},n_{2})|n_{1}\leq g^{2a}_{2}(n_{2}),~{}n_{2}\leq
g^{2a}_{1}(n_{1})\right\\}$, where $n_{2}=g^{2a}_{1}(n_{1})$ and
$n_{1}=g^{2a}_{2}(n_{2})$ denotes the two stable separatrices
$\Phi^{1}({P}^{2a}_{1})$ and $\Phi^{2}(P^{2a}_{1})$ of $P^{2a}_{1}$,
respectively;

$S^{2a}_{in}=\left\\{(n_{1},n_{2})|n_{1}\leq\frac{p_{1}}{2},~{}n_{2}\leq\frac{p_{2}-M}{2}\right\\};$
$S^{2a}_{out}=S^{2a,1}_{out}\vee S^{2a,2}_{out}$, where

$S^{2a,1}_{out}=\left\\{(n_{1},n_{2})|n_{1}\leq\frac{p_{1}}{2},~{}n_{2}<b^{2a}_{2}(n_{1})\right\\}$
and $S^{2a,2}_{out}=\left\\{(n_{1},n_{2})|\frac{p_{1}}{2}\leq n_{1}\leq
p_{1},~{}n_{2}<r^{2a}_{1}(n_{1})\right\\}.$

  2. (2)

Under Condition $(\mathbb{H}_{3})$,
$S^{2a}_{A}=\left\\{(n_{1},n_{2})|n_{1}\leq\frac{p_{1}}{2},\
n_{2}<g_{3}^{2a}(n_{1})\right\\}$, where $n_{2}=g_{3}^{2a}(n_{1})$ denotes the
stable separatrix $\Phi({P}^{2a}_{2})$ of $P^{2a}_{2}$;
$S^{2a}_{in}=\left\\{(n_{1},n_{2})|n_{1}\leq\frac{p_{1}}{2},~{}n_{2}<\frac{p_{2}+M}{2}\right\\}$;
$S^{2a}_{out}=S^{2a,1}_{out}$;

  3. (3)

Under Condition $(\mathbb{H}_{4})$,
$S^{2a}_{in}=S^{2a}_{A}=S^{2a}_{out}=\left\\{(n_{1},n_{2})|n_{1}\leq\frac{p_{1}}{2},\
n_{2}<\frac{{p_{2}}+M}{2}\right\\}$.

Note that $n_{i}\geq 0$ ($i=1,2$) holds all the time.

The geometries of $S^{2a}_{A}$ in Theorem 3.2 are illustrated as the pale pink
region shown in Fig. 7.

(a) $\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$

(b) $\mathbb{K}^{2a}\wedge\mathbb{H}_{2}$

(c) $\mathbb{K}^{2a}\wedge\mathbb{H}_{3}$

(d) $\mathbb{K}^{2a}\wedge\mathbb{H}_{4}$

Figure 7: Attraction region of the $P^{2a}_{1}$ for system (3.1) under
Condition $\mathbb{K}^{2a}\wedge\mathbb{H}_{i}$ ($i=1,2,\cdots,4$). Red
circle: undesired saddle-node equilibrium. Blue circle: desired saddle-node
equilibrium. The yellow curves: the separatrices, also serving as boundaries
of attraction region. Purple lines and bright blue lines denote the boundary
of inner and outer estimations of attraction regions, respectively.

(a) $\mathbb{K}^{2b}\wedge\mathbb{H}_{1}$

(b) $\mathbb{K}^{2b}\wedge\mathbb{H}_{2}$

(c) $\mathbb{K}^{2b}\wedge\mathbb{H}_{3}$

(d) $\mathbb{K}^{2b}\wedge\mathbb{H}_{4}$

Figure 8: Attraction region of the $P^{2b}_{1}$ for system (1) under Condition
$\mathbb{K}^{2b}\wedge\mathbb{H}_{i}$ ($i=1,2,\cdots,4$). Red circle:
undesired saddle-node equilibrium. Blue circle: desired saddle-node
equilibrium. The yellow curves: the separatrices, also serving as boundaries
of attraction region. Purple lines and bright blue lines denote the boundary
of inner and outer estimations of attraction regions, respectively.

Similarly, for system (1) under Condition $(\mathbb{K}^{2b})$, we can obtain
the inner ($S^{2b}_{in}$) and outer ($S^{2b}_{out}$) estimations of attraction
regions for $P^{2b}_{1}$ under four different conditions, as shown in Fig. 8.

For the four-equilibrium system, we also try to obtain the inner
($S^{4}_{in}$) and outer ($S^{4}_{out}$) estimations of attraction regions,
and the detailed analysis can be found in App. 7.3. Note that for the four-
equilibrium system, the attraction regions become more intricate due to the
increased number of equilibria. In App. 7.3, through a novel classification
method, we obtained 20 distinct inner and outer estimations of attraction
regions, which can be summarized as Theorems 7.2, 7.3, and 7.4 in App. 7.3.
The geometries of the attraction region for the four-equilibrium system are
illustrated as the pale pink region shown in Fig. 9.

(a) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$

(b) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$

(c) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$

(d) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$

(e) $\hat{\mathbb{K}}^{4}_{2}\wedge\mathbb{H}_{1}$

(f) $\hat{\mathbb{K}}^{4}_{2}\wedge\mathbb{H}_{2}$

(g) $\hat{\mathbb{K}}^{4}_{3}\wedge\mathbb{H}_{1}$

(h) $\hat{\mathbb{K}}^{4}_{3}\wedge\mathbb{H}_{3}$

(i) $\hat{\mathbb{K}}^{4}_{4}\wedge\mathbb{H}_{1}$, type 1

(j) $\hat{\mathbb{K}}^{4}_{4}\wedge\mathbb{H}_{2}$, type 1

(k) $\hat{\mathbb{K}}^{4}_{5}\wedge\mathbb{H}_{1}$, type 1

(l) $\hat{\mathbb{K}}^{4}_{5}\wedge\mathbb{H}_{3}$, type 1

(m) $\hat{\mathbb{K}}^{4}_{4}\wedge\mathbb{H}_{1}$, type 2

(n) $\hat{\mathbb{K}}^{4}_{4}\wedge\mathbb{H}_{2}$, type 2

(o) $\hat{\mathbb{K}}^{4}_{5}\wedge\mathbb{H}_{1}$, type 2

(p) $\hat{\mathbb{K}}^{4}_{5}\wedge\mathbb{H}_{3}$, type 2

(q) $\hat{\mathbb{K}}^{4}_{4}\wedge\mathbb{H}_{1}$, type 3

(r) $\hat{\mathbb{K}}^{4}_{4}\wedge\mathbb{H}_{2}$, type 3

(s) $\hat{\mathbb{K}}^{4}_{5}\wedge\mathbb{H}_{1}$, type 3

(t) $\hat{\mathbb{K}}^{4}_{5}\wedge\mathbb{H}_{3}$, type 3

Figure 9: Attraction region of the $P^{4}_{1}$ of system (1) under Condition
$\hat{\mathbb{K}}^{4}_{s}\wedge\mathbb{H}_{i}$ ($s=1,2,\cdots,5$ and
$i=1,2,\cdots,4$). Green circle: desired uncongested stable equilibrium.
Yellow circle: saddle equilibrium. Yellow-brown circle: unstable equilibrium.
Yellow curves with arrows: the separatrices, also serving as boundaries of
attraction region. Purple lines and bright blue lines denote the boundary of
inner and outer estimations of attraction regions, respectively.

##  4 Ecological resilience control design

In the previous two sections, we obtain the inner ($S_{in}$) and outer
($S_{out}$) estimations of attraction regions with explicit algebraic
expressions under each condition for a multi-equilibria MFD-represented
traffic dynamic system (1). Note that for two-equilibria system, $S_{in}$
corresponds to $S^{2a}_{in}$ or $S^{2b}_{in}$, and $S_{out}$ corresponds to
$S^{2a}_{out}$ or $S^{2b}_{out}$; for two-equilibria system, $S_{in}$
corresponds to $S^{4}_{in}$ and $S_{out}$ corresponds to $S^{4}_{out}$. While
In this section, we will design two distinct resilience control schemes,
aiming to facilitate the recoverability of the system to its respective steady
states while minimally altering the original landscape, utilizing the obtained
$S_{in}$ and $S_{out}$. Sec. 4.1 illustrates the resilience control schemes
utilizing $S_{in}$, denoted as RCS-1, and Sec. 4.2 provides the resilience
control schemes utilizing $S_{out}$, denoted as RCS-2.

The trajectories starting from inside attraction region (or inner estimations
of attraction regions) can converge to the original equilibrium, while for
those that starting from outside attraction region (or outer estimations of
attraction regions) will move away from the equilibrium. Our primary task is
to control the trajectories that move away from the equilibrium point and
guide them towards an alternative steady state. For multi-equilibria systems
(1), we set up an alternative steady state $\bm{n}^{1}$=$(ns_{1},ns_{2})$ in
$\mathbb{D}\setminus S_{in}$ or $\mathbb{D}\setminus S_{out}$, respectively.
Then, the control schemes (RCS-1 and RCS-2) shall steer the trajectories to
the $\varepsilon-Neighborhood$ of $\bm{n}^{1}$ within the recovery time
$t_{f}$. Note that the control does not necessarily bring the system to the
equilibrium state $\bm{n}^{0}$, which allows the system to be adaptive (to a
feasible steady state) and recoverable (to less congested state
progressively).

Generally speaking, there are three main differences between the RCS-1 and
RCS-2 schemes. With respect to control “intensity”, RCS-1 is more aggressive
than RCS-2. That is, the control region for RCS-1 is much larger than RCS-2,
as control action $\bm{U}_{1}(\bm{n})$ ($\bm{U}_{2}(\bm{n})$) are imposed in
$\mathbb{D}\setminus S_{in}$ ($\mathbb{D}\setminus S_{out}$) under RCS-1
(RCS-2) and we have $S_{in}\subset S_{out}$. With respect to control
mechanism, RCS-1 has a universal control algorithm applied to all alternative
steady states, while RCS-2 applies three forms of control, depending on the
location of the alternative steady state. Here, we choose three control
methods to prevent the emergence of new equilibrium points on the boundary of
outer estimation of attraction regions. For instance, trajectories from the
boundary of $S_{out}$ tend to escape $S_{out}$, and one control scheme might
lead to intersections between trajectories from $\mathbb{D}\setminus S_{out}$
and trajectories from $S_{out}$. This could potentially generate new
equilibrium points on the boundary of the outer estimation of attraction
regions, complicating theoretical derivations and causing control
inefficiencies. With respect to control effect, the trajectories starting from
$\mathbb{D}\setminus S_{out}$ will move towards $\bm{n}^{1}$, while
trajectories starting from $\mathbb{D}\setminus S_{in}$ will move towards
$\bm{n}^{1}$ (Fig. 10(b)), or $\bm{n}^{0}$ (Fig. 10(c)).

###  4.1 Control formulation of RCS-1

The domain $\mathbb{D}$ can be divided into two subregions: $S_{in}$ and
$\mathbb{D}\setminus S_{in}$, as shown in Fig. 10. If the starting points are
in $S_{in}$, CPC $u_{1}$ and $u_{2}$ that satisfy Conditions
$\mathbb{K}^{2a}$, $\mathbb{K}^{2b}$ or $\mathbb{K}^{4}$ are conducted.
According to Theorems 3.2, 7.2 and 7.4, trajectories from these starting
points will move towards the $\bm{n}^{0}$. Instead, if the initial states are
in $\mathbb{D}\setminus S_{in}$, new control action $\bm{U}^{in}(\bm{n})$ is
activated for recovering trajectories starting from $\mathbb{D}\setminus
S_{in}$ to the $\varepsilon-Neighborhood$ of the target point within the
recovery time $t_{f}$, The target point can be the uncongested equilibrium
$\bm{n}^{0}$ or the alternative steady state $\bm{n}^{1}$ in
$\mathbb{D}\setminus S_{in}$. Specifically, for system (1), if the starting
states are in $\mathbb{D}\setminus S_{in}$, we apply the following control
$\bm{U}^{in}(n_{1},n_{2})=\left(U^{in}_{1}(n_{1},n_{2}),U^{in}_{2}(n_{1},n_{2})\right)$,
defined by:

$\begin{split}U^{in}_{1}(n_{1},n_{2})=-\gamma_{1}\left(n_{1}-ns_{1}\right)-{u_{2}}{G_{2}}{(n_{2}(t))}+{G_{1}}{(n_{1}(t))}-d_{1},\\\
U^{in}_{2}(n_{1},n_{2})=-\gamma_{2}\left(n_{2}-ns_{2}\right)-{u_{1}}{G_{1}}{(n_{1}(t))}+{G_{2}}{(n_{2}(t))}-d_{2}.\end{split}$
(14)

where $\gamma_{i}=\frac{2G_{i,max}}{p_{i}}$. By adding proposed controller
(14) to system (1), the trajectories starting from $\mathbb{D}\setminus
S_{in}$ will move towards $\bm{n}^{1}$ directly (Fig. 10(b)), or intersect
$S_{in}$ (Fig. 10(c)). Once intersecting $S_{in}$, we rest the control, i.e.,
$U^{in}_{1}=U^{in}_{2}=0$ immediately, as such that the trajectories will
follow the initial system (1). According to Theorems 3.2, 7.2 and 7.4,
trajectories from these starting points will move towards the $\bm{n}^{0}$.

(a) CPC

(b) RCS-1

(c) RCS-1

Figure 10: Schematic representation of phase portraits under CPC (a) and RCS-1
((b) and (c)). Notice that in (c), the trajectories under RCS-1 intersect
$S_{in}$. Once intersecting $S_{in}$, we rest the control $U^{in}=0$
immediately, then trajectories spontaneously move towards the $\bm{n}^{0}$.

By doing so, we have developed a switched controlled system (She and Xue
(2014), Skafidas et al. (1999), Gao et al. (2022), She et al. (2020)) as:

$\dot{\bm{n}}=\bm{F}(\bm{n})+\bm{U}_{1}(\bm{n}),$

where $\bm{U}_{1}(\bm{n})\equiv\bm{0}$ when $\bm{n}\in S_{in}$,
$\bm{U}_{1}(\bm{n})=\bm{U}^{in}(n_{1},n_{2})$ when $\bm{n}\in$
$\mathbb{D}\setminus S_{in}$.

###### Remark 4.1

For multi-equilibria systems (1), certain inner estimations of attraction
regions do not contain its partial boundaries $l_{ex}$ (e.g. Figs. 6(d), 8(d),
9(b) etc.). This is because the trajectories starting from the partial
boundaries will deviate from the uncongested equilibrium point. In such cases,
$\bm{U}^{in}(\bm{n})$ are imposed until the trajectories intersect with
$\hat{l}_{ex}$, where $\hat{l}_{ex}$ is obtained by moving the line $l_{ex}$
down (or left) by $\varepsilon$.

###  4.2 Control formulation of RCS-2

RCS-2 is designed based on $S_{out}$. Likewise, the domain $\mathbb{D}$ can be
divided into $S_{out}$ and $\mathbb{D}\setminus S_{out}$, and
$\mathbb{D}\setminus S_{out}$ can be further divided into $R_{1}$, $R_{2}$ and
$R_{3}$, as illustrated in Fig. 11. If the starting states are in $S_{out}$,
CPC $u_{1}$ and $u_{2}$ that satisfy Conditions $\mathbb{K}^{2a}$,
$\mathbb{K}^{2b}$ or $\mathbb{K}^{4}$ will be conducted. According to Theorems
3.2, 7.2 and 7.4, trajectories from these starting points will move towards
$\bm{n}^{0}$ or escape from the upper boundaries of $S_{out}$ and enter
$\mathbb{D}\setminus S_{out}$. Once the states enter in $\mathbb{D}\setminus
S_{out}$, control action $\bm{U}^{out}(\bm{n})$ will be activated for
recovering trajectories to the $\varepsilon-Neighborhood$ of $\bm{n}^{1}$ in
$\mathbb{D}\setminus S_{out}$ within the recovery time $t_{f}$.

(a) CPC

(b) RCS-2, Case 1

(c) RCS-2, Case 2

(d) RCS-2, Case 3

Figure 11: Schematic representation of phase portraits under CPC (a) and RCS-2
((b), (c) and (d)). Note that there are three cases of control action,
depending on the location of $\bm{n}^{1}$. Three cases of control action are
chosen to prevent the emergence of new equilibrium points on the boundary of
$S_{out}$.

Specifically, for a given system (1), there are three cases of control action
$\bm{U}^{out}(\bm{n})$, depending on the location of $\bm{n}^{1}$. First
consider Case 1, where $\bm{n}^{1}$ lies in $R_{1}$, as depicted in Fig.
11(b). In this case, proposed RCS-2 is as follows: in addition to $u_{1}$ and
$u_{2}$, we will apply an additional control
$\bm{U}^{out}_{1}(n_{1},n_{2})=\left(U^{out}_{1,1}(n_{1},n_{2}),U^{out}_{1,2}(n_{1},n_{2})\right)$
in $R_{1}\vee R_{2}$, and
$\bm{U}^{out}_{2}(n_{1},n_{2})=\left(U^{out}_{2,1}(n_{1},n_{2}),U^{out}_{2,2}(n_{1},n_{2})\right)$
in $R_{3}$, where $\bm{U}^{out}_{1}(n_{1},n_{2})$ are defined by:

$\begin{split}U^{out}_{1,1}(n_{1},n_{2})=-\gamma_{1}\left(n_{1}-ns_{1}\right)-{u_{2}}{G_{2}}{(n_{2}(t))}+{G_{1}}{(n_{1}(t))}-d_{1},\\\
U^{out}_{1,2}(n_{1},n_{2})=-\gamma_{2}\left(n_{2}-ns_{2}\right)-{u_{1}}{G_{1}}{(n_{1}(t))}+{G_{2}}{(n_{2}(t))}-d_{2},\end{split}$
(15)

and $\bm{U}^{out}_{2}(n_{1},n_{2})$ are defined by:

$\displaystyle U^{out}_{2,1}(n_{1},n_{2})$ $\displaystyle=$
$\displaystyle-{u_{2}}{G_{2}}{(n_{2}(t))}+{G_{1}}{(n_{1}(t))}-d_{1},$
$\displaystyle U^{out}_{2,2}(n_{1},n_{2})$ $\displaystyle=$
$\displaystyle-\gamma_{2}n_{2}-{u_{1}}{G_{1}}{(n_{1}(t))}+{G_{2}}{(n_{2}(t))}-d_{2}.$
(16)

by adding proposed controller (4.2) to system (1), the trajectories starting
from $R_{3}$ will vertically escape from the upper boundary of $R_{3}$ and
enter $R_{2}$. Moreover, by adding proposed controller (15) to system (1), the
trajectories starting from $R_{1}\vee R_{2}$ will move towards $\bm{n}^{1}$
directly, as shown in Fig. 11(b). Thus, any trajectories entering
$\mathbb{D}\setminus S_{out}$ will move towards $\bm{n}^{1}$ under the
controller.

Then consider Case 2, where $\bm{n}^{1}$ lies in $R_{2}$. In this case,
$\bm{U}^{out}_{1}(n_{1},n_{2})=\left(U^{out}_{1,1}(n_{1},n_{2}),U^{out}_{1,2}(n_{1},n_{2})\right)$
in $\mathbb{D}\setminus S_{out}$ will be applied. Subsequently, any
trajectories going into $\mathbb{D}\setminus S_{out}$ will go to $\bm{n}^{1}$
under our proposed controller, as shown in Fig. 11(c).

Finally consider Case 3, where $\bm{n}^{1}$ lies in $R_{3}$, as shown in Fig.
11(d). similarly, we will apply control action
$\bm{U}^{out}_{1}(n_{1},n_{2})=\left(U^{out}_{1,1}(n_{1},n_{2}),U^{out}_{1,2}(n_{1},n_{2})\right)$
in $R_{3}\vee R_{2}$, and
$\bm{U}^{out}_{3}(n_{1},n_{2})=\left(U^{out}_{3,1}(n_{1},n_{2}),U^{out}_{3,2}(n_{1},n_{2})\right)$
in $R_{1}$, where $\bm{U}^{out}_{3}(n_{1},n_{2})$ are defined by:

$\displaystyle U^{out}_{3,1}(n_{1},n_{2})$ $\displaystyle=$
$\displaystyle-\gamma_{1}n_{1}-{u_{2}}{G_{2}}{(n_{2}(t))}+{G_{1}}{(n_{1}(t))}-d_{1},$
$\displaystyle U^{out}_{3,2}(n_{1},n_{2})$ $\displaystyle=$
$\displaystyle-{u_{1}}{G_{1}}{(n_{1}(t))}+{G_{2}}{(n_{2}(t))}-d_{2}.$ (17)

by adding proposed controller (4.2) to system (1), the trajectories starting
from $R_{1}$ will horizontally escape from the right boundary of $R_{1}$ and
enter $R_{2}$. Moreover, by adding proposed controller (15) to system (1), the
trajectories starting from $R_{2}\vee R_{3}$ will move towards $\bm{n}^{1}$
directly, as shown in Fig. 11(d). Thus, any trajectories entering
$\mathbb{D}\setminus S_{out}$ will move towards $\bm{n}^{1}$ under our
proposed controller.

By doing so, we have developed another switched controlled system (She and Xue
(2014), Skafidas et al. (1999), Gao et al. (2022), She et al. (2020)) as:

$\dot{\bm{n}}=\bm{F}(\bm{n})+\bm{U}_{2}(\bm{n}),$

where $\bm{U}_{2}(\bm{n})\equiv\bm{0}$ when $\bm{n}\in S_{out}$, and
$\bm{U}_{2}(\bm{n})$ features different formulation in different cases when
$\bm{n}\in$ $\mathbb{D}\setminus S_{out}$. Specifically,

  * •

in case 1, $\bm{U}_{2}(\bm{n})=\bm{U}^{out}_{1}(n_{1},n_{2})$ when $\bm{n}\in$
$R_{1}\vee R_{2}$ and $\bm{U}_{2}(\bm{n})=\bm{U}^{out}_{2}(n_{1},n_{2})$ when
$\bm{n}\in$ $R_{3}$;

  * •

in case 2, $\bm{U}_{2}(\bm{n})=\bm{U}^{out}_{1}(n_{1},n_{2})$ when $\bm{n}\in$
$\mathbb{D}\setminus S_{out}$;

  * •

in case 3, $\bm{U}_{2}(\bm{n})=\bm{U}^{out}_{1}(n_{1},n_{2})$ when $\bm{n}\in$
$R_{2}\vee R_{3}$ and $\bm{U}_{2}(\bm{n})=\bm{U}^{out}_{3}(n_{1},n_{2})$ when
$\bm{n}\in$ $R_{1}$.

###### Remark 4.2

$\bm{U}^{in}(n_{1},n_{2})$ and $\bm{U}^{out}_{2}(n_{1},n_{2})$ would have the
same formulation if the alternative steady state $\bm{n}^{1}$ in both
controllers are same. This implies that there may exist initial states, such
that the system may end up with identical recovery trajectories under RCS- 1
and RCS-2, respectively.

Furthermore, for a given system (1), if a stable equilibrium point exists
(i.e., Condition $(\mathbb{K}^{4})$ holds), a local Lyapunov function can be
found (Wang et al. (2021), Willems (1971), Blanchini (1995)) for the new
controlled system with RCS-2. Specifically, denoting
$\bm{n}=(n_{1},n_{2})^{\mathrm{T}}$, we can prove that the function
$V=\left(\bm{n}-\bm{n}^{0}\right)^{\mathrm{T}}\bm{B}\left(\bm{n}-\bm{n}^{0}\right)\left(\bm{n}-\bm{n}^{1}\right)^{\mathrm{T}}\left(\bm{n}-\bm{n}^{1}\right)$
is a local Lyapunov function for system (1) with RCS-2, where
$\bm{B}=\left(\bm{P}^{-1}\right)^{\mathrm{T}}\bm{P}^{-1}$ and

$\bm{P}=\left(\begin{array}[]{cc}p_{11}&p_{12}\\\
p_{21}&p_{22}\end{array}\right)=\left(\begin{array}[]{cc}\frac{-a_{1}\,s_{1}+a_{2}\,s_{2}+Q}{2\,a_{1}\,s_{1}\,u_{1}}&\frac{-a_{1}\,s_{1}+a_{2}\,s_{2}-Q}{2\,a_{1}\,s_{1}\,u_{1}}\\\
1&1\end{array}\right)$

with
$Q=\sqrt{{a_{1}}^{2}\,{s_{1}}^{2}+{a_{2}}^{2}\,{s_{2}}^{2}-2\,a_{1}\,a_{2}\,s_{1}\,s_{2}+4\,a_{1}\,a_{2}\,s_{1}\,s_{2}\,u_{1}\,u_{2}}$;
$s_{1}=\frac{1}{2}\sqrt{{p_{1}}^{2}-\frac{4{d_{1}}+4{u_{2}}{d_{2}}}{{a_{1}}(1-{u_{1}}{u_{2}})}}$;
$s_{2}=\frac{1}{2}\sqrt{{p_{2}}^{2}-\frac{4{d_{2}}+4{u_{1}}{d_{1}}}{{a_{2}}(1-{u_{1}}{u_{2}})}}$.
Obviously, $V>0$ holds for all $\bm{n}\in\mathbb{D}$ except for
$\bm{n}=\bm{n}^{0}$ and $\bm{n}=\bm{n}^{1}$. Moreover, $V=0$ holds if and only
if $\bm{n}=\bm{n}^{0}$ and $\bm{n}=\bm{n}^{1}$. In addition, in the
neighborhood of $\bm{n}^{0}$, we have $\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}<0$
holds, and detailed proofs are as follows:

Letting $\bm{n}-\bm{n}^{0}=\bm{P}\bm{y}$, where
$\bm{y}=(y_{1},y_{2})^{\mathrm{T}}$, we can transform the system (1) into:

$\begin{split}\frac{\mathrm{d}{y}_{1}(t)}{\mathrm{d}t}=\lambda_{1}{y_{1}(t)}+H_{1}({y}_{1}(t),{y}_{2}(t)),\\\
\frac{\mathrm{d}{y}_{2}(t)}{\mathrm{d}t}=\lambda_{2}{y_{2}(t)}+H_{2}({y}_{1}(t),{y}_{2}(t)),\end{split}$
(18)

where $\lambda_{1}=-a_{1}\,s_{1}-a_{2}\,s_{2}+Q<0$;
$\lambda_{2}=-a_{1}\,s_{1}-a_{2}\,s_{2}-Q<0$; $H_{1}({y}_{1}(t),{y}_{2}(t))$
and $H_{2}({y}_{1}(t),{y}_{2}(t))$ represent second-order remained terms.

By the above coordinate transformation, denoting
$\bm{\delta}=\bm{n}^{1}-\bm{n}^{0}=(\delta_{1},\delta_{2})^{\mathrm{T}}$, we
can obtain that:

$\displaystyle{V}(t)$ $\displaystyle=$
$\displaystyle{\left(\bm{n}-\bm{n}^{0}\right)^{\mathrm{T}}\bm{B}\left(\bm{n}-\bm{n}^{0}\right)}\left(\bm{n}-\bm{n}^{1}\right)^{\mathrm{T}}\left(\bm{n}-\bm{n}^{1}\right),$
(19) $\displaystyle=$
$\displaystyle{(\bm{P}\bm{y})^{\mathrm{T}}\bm{B}(\bm{P}\bm{y})}\left(\bm{P}\bm{y}-\bm{\delta}\right)^{\mathrm{T}}\left(\bm{P}\bm{y}-\bm{\delta}\right),$
$\displaystyle=$
$\displaystyle{\bm{y}^{\mathrm{T}}\bm{y}}\left(\bm{P}\bm{y}-\bm{\delta}\right)^{\mathrm{T}}\left(\bm{P}\bm{y}-\bm{\delta}\right),$
$\displaystyle=$ $\displaystyle V_{1}\cdot V_{2}$

where $V_{1}={y_{1}}^{2}+{y_{2}}^{2}$ and
$V_{2}=(p_{11}y_{1}+p_{12}y_{2}-\delta_{1})^{2}+(p_{21}y_{1}+p_{22}y_{2}-\delta_{2})^{2}$.
Then we have:

$\displaystyle\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}$ $\displaystyle=$
$\displaystyle\frac{\partial{V}}{\partial{y}_{1}}\cdot\frac{\mathrm{d}{y}_{1}(t)}{\mathrm{d}{t}}+\frac{\partial{V}}{\partial{y}_{2}}\cdot\frac{\mathrm{d}{y}_{2}(t)}{\mathrm{d}{t}},$
(20) $\displaystyle=$ $\displaystyle
2(\delta_{1}^{2}+\delta_{2}^{2})(\lambda_{1}y_{1}(t)^{2}+\lambda_{2}y_{2}(t)^{2})+H_{v}^{1}(y_{1}(t),y_{2}(t)).$

where $H_{v}^{1}(y_{1}(t),y_{2}(t))$ represents the high-order remained terms,
and the lowest power of it is cubic. Obviously,
$\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}<0$ holds for arbitrary $(y_{1},y_{2})$
in the neighborhood of $(0,0)$ except $(0,0)$, as we have
$\lambda_{2}<\lambda_{1}<0$. That is, $\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}<0$
holds for arbitrary $(n_{1},n_{2})$ in the neighborhood of $\bm{n}^{0}$,
implying function ${V}$ serves as a local Lyapunov function for the system in
the neighbor of $\bm{n}^{0}$.

Similarly, we can prove that $\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}<0$ holds in
the neighborhood of $\bm{n}^{0}$. Specifically, letting $z_{1}=n_{1}-ns_{1}$
and $z_{2}=n_{2}-ns_{2}$, we have

$\displaystyle\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}=-2(\delta_{1}\,\delta_{2})B(\delta_{1}\,\delta_{2})^{\mathrm{T}}(z_{1}(t)^{2}+z_{2}(t)^{2})+H_{v}^{2}(z_{1}(t),z_{2}(t)).$
(21)

where $H_{v}^{2}(z_{1}(t),z_{2}(t))$ represents the high-order remained terms,
and the lowest power of it is cubic. Obviously,
$\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}<0$ holds for arbitrary $(z_{1},z_{2})$
in the neighborhood of $(0,0)$ except $(0,0)$, as $B$ is a symmetric positive
definite matrix. In other words, $\frac{\mathrm{d}{V}(t)}{\mathrm{d}t}<0$
holds for arbitrary $(n_{1},n_{2})$ in the neighborhood of $\bm{n}^{1}$,
implying function ${V}$ is also local Lyapunov function for system in the
neighbor of $\bm{n}^{1}$.

##  5 Case studies

In this section, the control performance of the proposed resilient control by
comparative case studies are demonstrated. In particular, we benchmark with
the classical perimeter control (CPC) and the traffic-resilience oriented
control proposed in Gao et al. (2022).

###  5.1 Numerical verification

A two-region system is set up for case studies, where the network and traffic
settings follow the Downtown San Francisco road network from Aboudolas and
Geroliminis’ work (Aboudolas and Geroliminis (2013)). We adopt the function
form of the MFD models from Gao et al. (2022) for comparison purposes, i.e.,
$G_{1}(n_{1})={-\frac{28}{289}}{n_{1}}(n_{1}-850)$ for region 1 and
$G_{2}(n_{2})={-\frac{32}{961}}{n_{2}}(n_{2}-1550)$ for region 2, with maximum
capacity $G_{1,max}=7\times 10^{4}$, $G_{2,max}=8\times 10^{4}$ $[veh/h]$ and
jam accumulations $p_{1}=1700$, $p_{2}=3100$ $[veh]$.

Based on this MFD data, we consider traffic flow dynamics following Eq. (1)
and (2) as defined in Sec. 2.1. We first validate the theoretical results of
inner and outer estimation of attraction regions for both two and four
equilibria derived in Sec. 3 (refer to Theorems 3.2, 7.2 and 7.4, Figs. 7, 8
and 9), accumulating to a total of 28 distinct phase portraits (8 for two-
equilibria system and 20 for four-equilibria system). To achieve this, we
consider 28 different scenarios, labeled as Scenario 1 to Scenario 28. Each
scenario has distinct allowed flow rates $u_{i}$ and net fixed demand $d_{i}$
($i=1,2$), corresponding to the 28 conditions in Sec. 3. For example, Scenario
1 corresponds to Condition $\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$. We showcase
two scenarios (scenario 1 and scenario 9) in Fig. 12. Fig. 12(a) depicts the
phase portrait for Scenario 1: $u_{1}=0.4,u_{2}=0.5128,d_{1}=3\times
10^{4},d_{2}=5\times 10^{4}$ $[veh/h]$, satisfying Condition
$\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$, corresponding to the two-equilibrium
case. It verifies Theorem 1 (or Fig. 7(a)). Fig. 12(b) presents the phase
portrait for Scenario 9: $u_{1}=0.3,u_{2}=0.4,d_{1}=3\times
10^{4},d_{1}=5\times 10^{4}$ $[veh/h]$, satisfying Condition
$\mathbb{K}^{4}_{1}\wedge\mathbb{H}_{1}$, corresponding to the four-
equilibrium case. It verifies Theorem 2 (or Fig. 9(a)). Parameters and
numerical results for other scenarios can be found in App. 8: Fig. 22, Fig. 23
and Table 4. The results are also consistent with our theoretical findings.

Next, numerical phase portraits of the controlled systems under RCS-1 and
RCS-2 are demonstrated, as depicted in Fig. 13(b), 13(c), and 13(e), 13(f).
Here we use Scenario 9 as a showcase, which corresponds to the complex four-
equilibria case. For the simpler two-equilibria case, the control scheme can
also be applicable. Note that Case $i(i=1,2,3)$ denotes the case where the
alternative steady state $\bm{n}^{1}\in R_{i}$. The illustration of $R_{i}$
can be found in Fig. 11. Here, we take Case 1 ($\bm{n}^{1}\in R_{1}$) and Case
3 ($\bm{n}^{1}\in R_{3}$) as examples. Fig. 13(a) represents the results under
CPC control, serving as a benchmark control. Let us focus on the part
$\mathbb{D}\setminus S_{in}$ in Fig. 13(a), (b), and (e). In Fig. 13(a), a
small portion of the trajectories in $\mathbb{D}\setminus S_{in}$ enters
$S_{in}$, while a larger portion of the trajectories moves towards the upper
or right boundaries of the phase portrait, approaching grid-lock, i.e., high
traffic density but low outflow. In contrast, Fig. 13(b) and 13(e) demonstrate
that RCS-1 is capable of guiding trajectories starting from
$\mathbb{D}\setminus S_{in}$ towards either the non-congested equilibrium
point $\bm{n}^{0}$ or the alternative steady state $\bm{n}^{1}$ in
$\mathbb{D}\setminus S_{in}$, avoiding grid-lock. Similarly, in Fig. 13(a),
trajectories in $\mathbb{D}\setminus S_{out}$ all move towards the upper and
right boundaries of the phase portrait, while Fig. 13(c) and 13(f) confirm
that RCS-2 directs trajectories in $\mathbb{D}\setminus S_{out}$ towards the
alternative steady state $\bm{n}^{1}$ in $\mathbb{D}\setminus S_{out}$. These
results demonstrate the recoverability under RCS-1 and RCS-2. Furthermore, the
similar trajectories in $R1\vee R2$ (upper region in $\mathbb{D}\setminus
S_{out}$) in Fig. 13(b) and 13(c) confirm that there exist some initial states
such that the recovery trajectories starting from them under RCS-1 and RCS-2
are identical, as illustrated in Remark 2. This is because, in $R1\vee R2$,
RCS-1 and RCS-2 have the same control formulation.

(a) $\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$ (scenario 1)

(b) $\mathbb{K}^{4}_{1}\wedge\mathbb{H}_{1}$ (scenario 9)

Figure 12:  Numerical simulation phase portraits for system (1) under
Conditions (a) $\mathbb{K}^{2a}\wedge\mathbb{H}_{1}$ (scenario 1) and (b)
$\mathbb{K}^{4}_{1}\wedge\mathbb{H}_{2}$ (scenario 9). Red circle: undesired
saddle-node equilibrium. Blue circle: desired uncongested saddle-node
equilibrium. Green circle: desired uncongested stable equilibrium. Yellow
circle: saddle equilibrium. Yellow-brown circle: unstable equilibrium. Green
line with arrows: the trajectories. Purple line and light blue line represent
the inner and outer estimations of attraction regions, respectively. The
yellow lines represent the numerical attraction region boundaries, solid
yellow lines indicating the boundary line belongs to attraction region, while
dashed yellow lines indicating the boundary line does not belong to attraction
region. The agreement of this figure with Figs. 7(a) and Fig. 9(a) validate
the theoretical analysis. Figure 13: Numerical simulation phase portraits of
system (1) with (a) CPC, (b) RCS-1 (case 1), (c) RCS-2 (case 1), (d) RCS-
single, (e) RCS-1 (case 3) and (f) RCS-2 (case 3) under condition
$\mathbb{K}^{4}_{1}\wedge\mathbb{H}_{1}$ (scenario 9). Green circle:
uncongested equilibrium $\bm{n}^{0}$ and alternative steady state
$\bm{n}^{1}$. Pentagrams denote three situations (i.e., three initial recovery
states), black lines departing from the pentagrams are the corresponding
recovery trajectories. CPC cannot recover the system; RCS-single recovers to
the original equilibrium $\bm{n}^{0}=(481,926)$, while RCS-1 and RCS-2 can
recover to alternative stable states $\bm{n}^{1}$. Here we choose
$\bm{n}^{1}=(850,2274)$ for Case 1 and $\bm{n}^{1}=(1269,1550)$ for Case 3,
both satisfying the selectable range of alternative stable states Eq. 7. Note
that Case $i(i=1,2,3)$ denotes the case where the alternative steady state
$\bm{n}^{1}\in R_{i}$. The illustration of $R_{i}$ can be found in Fig. 11.

Moreover, for the case where the system has a stable equilibrium (we here use
scenario 9 as a showcase), we can find the function $V(t)$ defined in Sec.
4.2, as shown in Fig. 14(a), is a local Lyapunov function for the controlled
system with RCS-2. Moreover, the projection of $V(t)$ in the plane
$n_{1}-n_{2}$ is shown in the Fig. 14(b). Note that the local Lyapunov is only
applicable to four-equilibria system, and is not suitable for two-equilibria
system, as the there exist no stable equilibria in two-equilibria system.

(a) V function

(b) Projection of the V function

Figure 14: (a) V function with respect to $n_{1}$ and $n_{2}$. (b) Projection
of the V function on the $n_{1}-n_{2}$ plane. The black line with arrows
represents several trajectories of system with RCS-2. Green circle: the
equilibrium points of system with RCS-2 (case 2). Note that scenario 9 is
considered. Here, $\bm{n}^{0}=(481,926)$ and $\bm{n}^{1}=(1219,2174)$.

###  5.2 Comparison on proposed two control schemes with CPC and RCS-single

In this section, we compare our proposed control schemes (RCS-1 and RCS-2)
with two classical perimeter control strategies (RCS-single in Gao et al.
(2022) and CPC in Haddad and Geroliminis (2012)) to show the resilience of our
approaches. We maintain the same MFD setup as in the previous section, and use
scenario 9 as an example. Fig. 13(d) depicts the numerical simulation of phase
portrait for the control system with RCS-single, where trajectories within the
whole region recover to the original equilibrium $\bm{n}^{0}$. The key
distinction between RCS-single and proposed schemes lies in the trajectories
within the part: our proposed control schemes can recover to the alternative
steady state $\bm{n}^{1}$, while RCS-single can only recover to the original
equilibrium $\bm{n}^{0}$.

Then, we calculate the resilience measure (defined in Eq. (10)) under
different control schemes. Resilience measures are situation-dependent, where
situations referring to states (vehicle accumulations) after disturbances,
also known as initial recovery states. For instance, we assume state
perturbations suddenly occur at moment $t=0.5$ $[h]$. Situation $i$ denotes
the disturbed state at $t=0.5$ $[h]$ as $A_{i}$ ($i=1,2,3$). Here, we consider
three situations, i.e., $A_{1}=(500,2800)$, $A_{2}=(1500,2800)$ and
$A_{3}=(1428,1395)$, marked as three pentagrams in Fig. 13. Following these
trajectories (the black lines in Fig. 13) departing from the pentagrams, we
plot the vehicle accumulation $n_{1}(t)$, completion flow rate
$G_{1}(n_{1}(t))$, and the numerical diagrams of deviation from maximum
completion flow ($G_{1,max}$ and $G_{2,max}$) evolving with time under each
control scheme (CPC, RCS-1, RCS-2, RCS-single) in Fig. 15. Here we use Case 1
(corresponding to Fig. 13(b) and 13(c)) as an example to illustrate
performance the proposed control schemes, and the results for Case 2 and Case
3 can be found in App. 8, Fig. 24 and Fig. 25.

Figure 15: Vehicle accumulations (dotted line) and completion flow rates
(solid line) of Region 1 (((a), (b) and (c))) and Region 2 ((d), (e) and (f)).
The numerical diagrams of deviation from maximum complete flow $-D_{s}$ ((g),
(h) and (i)) evolving with time $t$ under CPC, RCS-1, RCS-2 and RCS-single.
Scenario 9 is considered. Note that for (a)-(f), we employ dual y-axes. The
left y-axis represents vehicle accumulation (density), while the right y-axis
represents completion flow. $A_{i}$ ($i=1,2,3$) are the same as the three
pentagrams in Fig. 13. The results in (a)(d)(g) ((b)(e)(h) and (c)(f)(i)) are
obtained by following the recovery trajectories (black lines in Fig. 13)
departing from the pentagram $A_{1}$ ($A_{2}$ and $A_{3}$). Reaffirming that
RCS-single, RCS-1, and RCS-2 can ensure recoverability, while CPC can not. The
resilience triangles in (g), (h) and (i) reveal that RCS-1 and RCS-2
outperform RCS-single, as RCS-1 and RCS-2 have smaller resilience triangles
(loss) in most scenarios.

Clearly, as depicted in Fig. 15(g), (h), and (i), traffic system collapses
rapidly under CPC (blue solid line) (within 0.006 $[h]$ for $A_{1}$, 0.01
$[h]$ for $A_{2}$ and 0.007 $[h]$ for $A_{3}$); under RCS-single (red solid
line), system initially rises slowly approaching maximum completion flow, then
rebounds and stabilizes at a deviation of 26140 $[veh/h]$ (within 0.70 $[h]$
for $A_{1}$, 1.61 $[h]$ time for $A_{2}$ and 1.04 $[h]$ for $A_{3}$); under
the proposed RCS-1 (yellow solid line), the system swiftly recovers to a
deviation of 17441 $[veh/h]$ (within 0.16 $[h]$ for $A_{1}$ and 0.16 $[h]$ for
$A_{2}$, coinciding with the green solid line), or it first rises rapidly
approaching the maximum completion flow, then rebounds and stabilizes at a
deviation of 26140 $[veh/h]$ (within 0.21 $[h]$ for $A_{3}$); under the
proposed RCS-2 (green solid line), the system swiftly recovers to a deviation
of 17441 $[veh/h]$ (within 0.16 $[h]$ for $A_{1}$ and 0.16 $[h]$ for $A_{2}$),
or it resists the continued deterioration of the system, guiding a swift
recovery, and stabilizes at a deviation of 17441 $[veh/h]$ (within 0.19 $[h]$
for $A_{3}$). Here, we take $A_{1}$ as an example to elucidate the reasons for
the aforementioned phenomena. Under CPC, the vehicle accumulation $n_{2}(t)$
(Fig. 15(d), blue dotted line) develops from 2800 to a jam vehicle
accumulation of 3100 $[veh]$ within 0.006 $[h]$, leading to the completion
flow $G_{2}(n_{2}(t))$ (blue solid line) tending towards 0. This causes Region
2 to enter a gridlock state from which it cannot recover. Under RCS-single,
the vehicle accumulation $n_{1}(t)$ (Fig. 15(a), red dotted line) initially
decreases slowly from 500 $[veh]$, then slowly increases before decreasing
again, finally stabilizing at the original equilibrium of 481 $[veh]$.
Correspondingly, the completion flow decreases slowly, then increases slowly,
before decreasing again and stabilizing at 56818 $[veh/h]$. Simultaneously,
the vehicle accumulation $n_{2}(t)$ (Fig. 15(d), red dotted line) decreases
slowly at first, then rapidly decreases and stabilizes at 926 $[veh]$.
Consequently, the completion flow $G_{2}(n_{2}(t))$ (Fig. 15(d), red solid
line) first rises slowly to maximum completion flow, then rapidly decreases
and stabilizes at 67045 $[veh/h]$. Therefore, the recoverability of RCS-single
is guaranteed, but with a longer duration. Under the proposed RCS-1 and RCS-2,
due to high travel demand, the vehicle accumulation $n_{1}(t)$ (Fig. 15(a),
green dotted line) rapidly increases from 500 $[veh]$ and stabilizes at 850,
resulting in the completion flow (Fig. 15(a), green solid line) rapidly rising
and stabilizing at the maximum completion flow of $7\times 10^{4}$ $[veh/h]$.
Simultaneously, the vehicle accumulation $n_{2}(t)$ (Fig. 15(b), green dotted
line) rapidly decreases from 2800 to 2274 $[veh]$, leading the completion flow
$G_{2}(n_{2}(t))$ quickly stabilized at 62559 $[veh/h]$. Consequently, RCS-1
and RCS-2 can ensure rapid recoverability, outperforming RCS-single. It is
noteworthy that in this scenario, RCS-1 and RCS-2 exhibit identical recovery
trajectories. We also provide the cases of RCS-1 and RCS-2 employing different
recovery trajectories, as shown in Fig. 15(c), (f), and (i).

For a more accurate comparison, based on Fig. 15 and Eq. (10), we integrate
the curve of $-D_{s}$ versus $t$ to obtain the value of resilience measure,
$R$ under RCS-1, RCS-2, and RCS-single. The results are outlined in Table 2.
Noteworthy, the absolute value of $R$ corresponds to the area of the
resilience triangle during the recovery period, and it signifies the loss of
actual completed trips compared to the maximum number of trips that could have
been completed during the recovery period. A larger resilience measure $R$
corresponds to a smaller area of the resilience triangle (less loss),
indicating greater resilience. As presented in Table 2, for all three
Situations ($A_{1}$, $A_{2}$ and $A_{3}$), resilience measures for RCS-1 and
RCS-2 are significantly larger than RCS-single, demonstrating exceptional
recovery capabilities.

Table 2: Resilience measure $[veh]$.

resilience control | $A_{1}$ | $A_{2}$ | $A_{3}$
---|---|---|---
Case 1 | Case 2 | Case 3 | Case 1 | Case 2 | Case 3 | Case 1 | Case 2 | Case 3
RCS-1 | -4031 | -4384 | -3110 | -4284 | -5646 | -3681 | -4250 | -3302 | -1710
RCS-2 | -4031 | -4384 | -5800 | -4284 | -6558 | -6615 | -3776 | -3302 | -1744
RCS-single | -24568 | -56501 | -20584 | -24568 | -56501 | -20584 | -24568 | -56501 | -20584

  * •

Case $i(i=1,2,3)$ denotes the case where the alternative steady state
$\bm{n}^{1}\in R_{i}$, as shown in Fig. 11. $A_{i}$ denotes three situations,
corresponding to three pentagrams in Fig. 13 and Fig. 15.

##  6 Conclusion

In this paper, we define the ecological resilience of urban traffic, and
further propose a resilience control methodology from the perspective of
ecological resilience. Specifically, the ecological resilience of urban
traffics is defined by the ability for a traffic system to resist uncertain
perturbations by shifting to alternative states. The resilience control
methodology comprises three aspects: portraying the recoverable scopes,
designing alternative steady states, controlling system to shift to
alternative steady states for adapting large disturbances. Among them, the
recoverable scopes are portrayed by inner and outer estimations of attraction
regions; the alternative steady states are set close to the optimal state and
outside the recoverable scopes of the original equilibrium; the controller
ensures the local stability of the alternative steady states, without changing
the trajectories inside the recoverable scopes of the original equilibrium as
much as possible. We implemented the proposed control framework in a two-
region traffic network described by parabolic MFD dynamic, and designed
proposed control schemes (RCS-1 and RCS-2). Comparisons with the classical
urban traffic resilience control schemes (CPC and RCS-single, etc.) show that,
proposed resilience control schemes (RCS-1 and RCS-2) has a better
adaptability and generates a greater resilience measure. Noteworthy, RCS-1
differs with RCS-2 in control intensity, control formulation and control
effect. Note that, for multi-region (more than two-region) MFD dynamics, our
resilience control idea can still be applied. The difficulty in implementing
the resilience control idea is that for N-dimensional system, the equilibrium
point will be more complex, and the attraction region maybe more difficult to
obtain. Moreover, the inner estimation of attraction region we can obtain may
shrink and the outer estimation of attraction region may expand, potentially
affecting the effectiveness of the regulation scheme.

There is still work to be done in the future. Given this paper considers
mainly constant MFD dynamics, and it is worth to extend the methodology of
this paper to more general time-varying MFD dynamics in the future. Moreover,
a feasible resilience control is provided in this paper, and future work will
be extended to find an optimal resilience control scheme.

This work was supported by the National Key R$\&$D Program of China
(No.2022YFA1005103), the National Natural Science Foundation of China
(12371452, 72225012, 72288101, 71822101, 12101030) and the Fundamental
Research Funds for the Central Universities.

## References

  * Aboudolas and Geroliminis (2013) Aboudolas K, Geroliminis N, 2013 _Perimeter and boundary flow control in multi-reservoir heterogeneous networks_. _Transportation Research Part B: Methodological_ 55:265–281.
  * Aghamohammadi and Laval (2020) Aghamohammadi R, Laval JA, 2020 _Dynamic traffic assignment using the macroscopic fundamental diagram: A review of vehicular and pedestrian flow models_. _Transportation Research Part B: Methodological_ 137:99–118.
  * Ampountolas, Zheng, and Geroliminis (2017) Ampountolas K, Zheng N, Geroliminis N, 2017 _Macroscopic modelling and robust control of bi-modal multi-region urban road networks_. _Transportation Research Part B: Methodological_ 104:616–637.
  * Blanchini (1995) Blanchini F, 1995 _Nonquadratic lyapunov functions for robust control_. _Automatica_ 31(3):451–461.
  * Bruneau et al. (2003) Bruneau M, Chang SE, Eguchi RT, Lee GC, O’Rourke TD, Reinhorn AM, Shinozuka M, Tierney K, Wallace WA, Von Winterfeldt D, 2003 _A framework to quantitatively assess and enhance the seismic resilience of communities_. _Earthquake spectra_ 19(4):733–752.
  * Daganzo, Gayah, and Gonzales (2011) Daganzo CF, Gayah VV, Gonzales EJ, 2011 _Macroscopic relations of urban traffic variables: Bifurcations, multivaluedness and instability_. _Transportation Research Part B: Methodological_ 45(1):278–288.
  * Dakos and Kéfi (2022) Dakos V, Kéfi S, 2022 _Ecological resilience: What to measure and how_. _Environmental Research Letters_ 17(4):043003.
  * Dantsuji, Fukuda, and Zheng (2021) Dantsuji T, Fukuda D, Zheng N, 2021 _Simulation-based joint optimization framework for congestion mitigation in multimodal urban network: a macroscopic approach_. _Transportation_ 48:673–697.
  * Folke (2006) Folke C, 2006 _Resilience: The emergence of a perspective for social ecological systems analyses_. _Global Environ Change_ 16(3):253–267.
  * Folke et al. (2010) Folke C, Carpenter SR, Walker B, Scheffer M, Chapin T, Rockström J, 2010 _Resilience thinking: integrating resilience, adaptability and transformability_. _Ecology and society_ 15(4).
  * Gao et al. (2022) Gao S, Li D, Zheng N, Hu R, She Z, 2022 _Resilient perimeter control for hyper-congested two-region networks with mfd dynamics_. _Transportation Research Part B: Methodological_ 156:50–75.
  * Geroliminis and Daganzo (2008) Geroliminis N, Daganzo CF, 2008 _Existence of urban-scale macroscopic fundamental diagrams: Some experimental findings_. _Transportation Research Part B: Methodological_ 42(9):759–770.
  * Gunderson (2000) Gunderson LH, 2000 _Ecological resilience–in theory and application_. _Annual review of ecology and systematics_ 425–439.
  * Haddad (2015) Haddad J, 2015 _Robust constrained control of uncertain macroscopic fundamental diagram networks_. _Transportation Research Part C_ 59(OCT.):323–339.
  * Haddad and Geroliminis (2012) Haddad J, Geroliminis N, 2012 _On the stability of traffic perimeter control in two-region urban cities_. _Transportation Research Part B: Methodological_ 46(9):1159–1176.
  * Haddad and Mirkin (2020) Haddad J, Mirkin B, 2020 _Resilient perimeter control of macroscopic fundamental diagram networks under cyberattacks_. _Transportation Research Part B: Methodological_ 132:44–59.
  * Haddad and Zheng (2020) Haddad J, Zheng Z, 2020 _Adaptive perimeter control for multi-region accumulation-based models with state delays_. _Transportation Research Part B: Methodological_ 137:133–153.
  * Holling (1973) Holling CS, 1973 _Resilience and stability of ecological systems_. _Annual review of ecology and systematics_ 1–23.
  * Holling (1996) Holling CS, 1996 _Engineering resilience versus ecological resilience_. _Engineering within ecological constraints_ 31(1996):32.
  * Hosseini, Barker, and Ramirez-Marquez (2016a) Hosseini S, Barker K, Ramirez-Marquez JE, 2016a _A review of definitions and measures of system resilience_. _Reliability Engineering & System Safety_ 145:47–61.
  * Hosseini, Barker, and Ramirez-Marquez (2016b) Hosseini S, Barker K, Ramirez-Marquez JE, 2016b _A review of definitions and measures of system resilience_. _Reliability Engineering & System Safety_ 145:47–61.
  * Huang et al. (2020a) Huang HJ, Xia T, Tian Q, Liu TL, Wang C, Li D, 2020a _Transportation issues in developing china’s urban agglomerations_. _Transport Policy_ 85:A1–A22.
  * Huang et al. (2020b) Huang Y, Xiong J, Sumalee A, Zheng N, Lam W, He Z, Zhong R, 2020b _A dynamic user equilibrium model for multi-region macroscopic fundamental diagram systems with time-varying delays_. _Transportation Research Part B: Methodological_ 131:1–25.
  * Johari et al. (2021) Johari M, Keyvan-Ekbatani M, Leclercq L, Ngoduy D, Mahmassani HS, 2021 _Macroscopic network-level traffic models: Bridging fifty years of development toward the next era_. _Transportation Research Part C: Emerging Technologies_ 131:103334.
  * Khwais and Haddad (2022) Khwais M, Haddad J, 2022 _Optimal presignal control for two-mode traffic at isolated signalized intersections_. _Transportation science_ 57(2):376–398.
  * Kouvelas, Saeedmanesh, and Geroliminis (2023) Kouvelas A, Saeedmanesh M, Geroliminis N, 2023 _A linear-parameter-varying formulation for model predictive perimeter control in multi-region mfd urban networks_. _Transportation Science_ 57(6):1496–1515.
  * Laval, Leclercq, and Chiabaut (2017) Laval JA, Leclercq L, Chiabaut N, 2017 _Minimal parameter formulations of the dynamic user equilibrium using macroscopic urban models: Freeway vs city streets revisited_. _Transportation research procedia_ 23:517–530.
  * Li, Yildirimoglu, and Ramezani (2021) Li Y, Yildirimoglu M, Ramezani M, 2021 _Robust perimeter control with cordon queues and heterogeneous transfer flows_. _Transportation Research Part C: Emerging Technologies_ 126:103043.
  * Liao (2012) Liao KH, 2012 _A theory on urban resilience to floods-a basis for alternative planning practices_. _Ecology and society_ 17(4).
  * Mariotte, Leclercq, and Laval (2017) Mariotte G, Leclercq L, Laval JA, 2017 _Macroscopic urban dynamics: Analytical and numerical comparisons of existing models_. _Transportation Research Part B: Methodological_ 101:245–267.
  * Martin-Breen and Anderies (2011) Martin-Breen P, Anderies JM, 2011 _Resilience: A literature review_ .
  * Mohajerpoor, Cai, and Ramezani (2023) Mohajerpoor R, Cai C, Ramezani M, 2023 _Optimal traffic signal control of isolated oversaturated intersections using predicted demand_. _IEEE Transactions on Intelligent Transportation Systems_ 24(1):815–826.
  * O’Neill et al. (1986) O’Neill RV, Deangelis DL, Waide JB, Allen TF, Allen GE, 1986 _A hierarchical concept of ecosystems_ (Princeton University Press).
  * Pimm (1984) Pimm SL, 1984 _The complexity and stability of ecosystems_. _Nature_ 307(5949):321–326.
  * Ramezani, Haddad, and Geroliminis (2015) Ramezani M, Haddad J, Geroliminis N, 2015 _Dynamics of heterogeneity in urban networks: aggregated traffic modeling and hierarchical control_. _Transportation Research Part B: Methodological_ 74:1–19.
  * Ratschan and She (2010) Ratschan S, She Z, 2010 _Providing a basin of attraction to a target region of polynomial systems by computation of lyapunov-like functions_. _SIAM Journal on Control and Optimization_ 48(7):4377–4394.
  * Rölfer, Celliers, and Abson (2022) Rölfer L, Celliers L, Abson D, 2022 _Resilience and coastal governance: knowledge and navigation between stability and transformation_. _Ecology and Society_ 27(2).
  * Scheffer et al. (1993) Scheffer M, Hosper SH, Meijer ML, Moss B, Jeppesen E, 1993 _Alternative equilibria in shallow lakes_. _Trends in ecology & evolution_ 8(8):275–279.
  * Scoones et al. (2020) Scoones I, Stirling A, Abrol D, Atela J, Charli-Joseph L, Eakin H, Ely A, Olsson P, Pereira L, Priya R, et al., 2020 _Transformations to sustainability: combining structural, systemic and enabling approaches_. _Current Opinion in Environmental Sustainability_ 42:65–75.
  * She and Xue (2014) She Z, Xue B, 2014 _Discovering multiple lyapunov functions for switched hybrid systems_. _SIAM Journal on Control and Optimization_ 52(5):3312–3340.
  * She et al. (2020) She Z, Zhang A, Lu J, Hu R, Sam Ge S, 2020 _Design proportional-integral-derivative/proportional-derivative controls for second-order time-varying switched nonlinear systems_. _International Journal of Robust and Nonlinear Control_ 30(5):1979–2000.
  * Skafidas et al. (1999) Skafidas E, Evans RJ, Savkin AV, Petersen IR, 1999 _Stability results for switched controller systems_. _Automatica_ 35(4):553–564.
  * Su et al. (2020) Su Z, Chow AH, Zheng N, Huang Y, Liang E, Zhong R, 2020 _Neuro-dynamic programming for optimal control of macroscopic fundamental diagram systems_. _Transportation Research Part C: Emerging Technologies_ 116:102628.
  * Tilman and Downing (1994) Tilman D, Downing JA, 1994 _Biodiversity and stability in grasslands_. _Nature_ 367(6461):363–365.
  * Tsitsokas, Kouvelas, and Geroliminis (2023) Tsitsokas D, Kouvelas A, Geroliminis N, 2023 _Two-layer adaptive signal control framework for large-scale dynamically-congested networks: Combining efficient max pressure with perimeter control_. _Transportation Research Part C: Emerging Technologies_ 152:104128.
  * Walker et al. (2004) Walker B, Holling CS, Carpenter SR, Kinzig A, 2004 _Resilience, adaptability and transformability in social–ecological systems_. _Ecology and society_ 9(2).
  * Walker and Salt (2012) Walker B, Salt D, 2012 _Resilience thinking: sustaining ecosystems and people in a changing world_ (Island press).
  * Wang, She, and Ge (2020a) Wang S, She Z, Ge SS, 2020a _Estimating minimal domains of attraction for uncertain nonlinear systems_. _IEEE Transactions on Systems, Man, and Cybernetics: Systems_ .
  * Wang, She, and Ge (2020b) Wang S, She Z, Ge SS, 2020b _Inner-estimating domains of attraction for nonpolynomial systems with polynomial differential inclusions_. _IEEE transactions on cybernetics_ .
  * Wang et al. (2021) Wang S, Wu W, Lu J, She Z, 2021 _Inner-approximating domains of attraction for discrete-time switched systems via multi-step multiple lyapunov-like functions_. _Nonlinear Analysis: Hybrid Systems_ 40:100993.
  * Willems (1971) Willems JC, 1971 _The generation of lyapunov functions for input-output stable systems_. _SIAM Journal on Control_ 9(1):105–134.
  * Yang, Menendez, and Zheng (2019) Yang K, Menendez M, Zheng N, 2019 _Heterogeneity aware urban traffic control in a connected vehicle environment: A joint framework for congestion pricing and perimeter control_. _Transportation Research Part C: Emerging Technologies_ 105:439–455.
  * Yang, Zheng, and Menendez (2018) Yang K, Zheng N, Menendez M, 2018 _Multi-scale perimeter control approach in a connected-vehicle environment_. _Transportation Research Part C: Emerging Technologies_ 94:32–49.
  * Zeng et al. (2020) Zeng G, Gao J, Shekhtman L, Guo S, Lv W, Wu J, Liu H, Levy O, Li D, Gao Z, et al., 2020 _Multiple metastable network states in urban traffic_. _Proceedings of the National Academy of Sciences_ 117(30):17528–17534.
  * Zeng et al. (2019) Zeng G, Li D, Guo S, Gao L, Gao Z, Stanley HE, Havlin S, 2019 _Switch between critical percolation modes in city traffic dynamics_. _Proceedings of the National Academy of Sciences_ 116(1):23–28.
  * Zeng et al. (2022) Zeng X, Yu Y, Yang S, Lv Y, Sarker MNI, 2022 _Urban resilience for urban sustainability: Concepts, dimensions, and perspectives_. _Sustainability_ 14\.
  * Zhang (2006) Zhang Z, 2006 _Qualitative theory of differential equations_ , volume 101 (American Mathematical Soc.).
  * Zheng et al. (2018) Zheng X, She Z, Liang Q, Li M, 2018 _Inner approximations of domains of attraction for a class of switched systems by computing lyapunov-like functions_. _International Journal of Robust and Nonlinear Control_ 28(6):2191–2208.
  * Zhong et al. (2020a) Zhong R, Cai H, Xu D, Chen C, Sumalee A, Pan T, 2020a _Dynamic feedback control of day-to-day traffic disequilibrium process_. _Transportation Research Part C: Emerging Technologies_ 114:297–321.
  * Zhong et al. (2018a) Zhong R, Chen C, Huang Y, Sumalee A, Lam W, Xu D, 2018a _Robust perimeter control for two urban regions with macroscopic fundamental diagrams: A control-lyapunov function approach_. _Transportation Research Part B: Methodological_ 117:687–707.
  * Zhong et al. (2018b) Zhong R, Huang Y, Chen C, Lam W, Xu D, Sumalee A, 2018b _Boundary conditions and behavior of the macroscopic fundamental diagram based network traffic dynamics: A control systems perspective_. _Transportation Research Part B: Methodological_ 111:327–355.
  * Zhong et al. (2020b) Zhong R, Xie X, Luo J, Pan T, Lam W, Sumalee A, 2020b _Modeling double time-scale travel time processes with application to assessing the resilience of transportation systems_. _Transportation research part B: methodological_ 132:228–248.

##  7 Global phase portrait and attraction region for four-equilibria system

we will consider the four-equilibria cases, derive its qualitative
characteristics under CPC ($u_{1}$ and $u_{2}$), and further provide an inner
and outer estimations of attraction regions. Specifically, local stability
verification of the four-equilibria system will be shown in App. 7.1; global
phase portrait derivation is performed in App. 7.2; and spontaneous attraction
region estimation of four-equilibria system will be included in App. 7.3.

###  7.1 Local stability verification

According to Condition $(\mathbb{K}^{4})$, we have:

$\begin{split}{u_{2}}d_{2}+d_{1}-\left(1-u_{1}u_{2}\right)G_{1,max}<0,\\\
{u_{1}}d_{1}+d_{2}-\left(1-u_{1}u_{2}\right)G_{2,max}<0.\end{split}$ (22)

For system (1) satisfying (22), the four equilibria are denoted as
$P^{4}_{m}=\left(p^{4}_{m,1},p^{4}_{m,2}\right)$, where $m=1,2,3,4$,
$p^{4}_{m,1}=\frac{{p_{1}}+(-1)^{m}\sqrt{{p_{1}}^{2}-\frac{4{d_{1}}+4{d_{2}}{u_{2}}}{{a_{1}}\left(1-{u_{1}}{u_{2}}\right)}}}{2}$;
$p^{4}_{m,2}=\frac{{p_{2}}-\sqrt{{p_{2}}^{2}-\frac{4{d_{2}}+4{d_{1}}{u_{1}}}{{a_{2}}\left(1-{u_{1}}{u_{2}}\right)}}}{2}$
when $m=1,2$ and
$p^{4}_{m,2}=\frac{{p_{2}}+\sqrt{{p_{2}}^{2}-\frac{4{d_{2}}+4{d_{1}}{u_{1}}}{{a_{2}}\left(1-{u_{1}}{u_{2}}\right)}}}{2}$
when $m=3,4$.

We first verify the local stability of the four equilibria for system (1)
satisfying (22). Under Condition ($\mathbb{K}^{4}\wedge\mathbb{H}$), we have
the following proposition for its four equilibria.

###### Proposition 7.1

Under Condition ($\mathbb{K}^{4}\wedge\mathbb{H}$), for system (1),
$P^{4}_{1}$ is a locally stable node, $P^{4}_{2}$ and $P^{4}_{3}$ are saddle
points, $P^{4}_{4}$ is an unstable node.

Proof: For system (1), we can get the derivative operator $A|_{P^{4}_{m}}$
($m=1,2,3,4$) at the equilibrium $P^{4}_{m}$ as follows:

$A|_{P^{4}_{m}}=Df|_{P^{4}_{m}}=\left(\begin{array}[]{ccccc}\frac{\partial
F_{1}(n_{1},n_{2})}{\partial n_{1}}&\frac{\partial
F_{1}(n_{1},n_{2})}{\partial n_{2}}\\\ \frac{\partial
F_{2}(n_{1},n_{2})}{\partial n_{1}}&\frac{\partial
F_{2}(n_{1},n_{2})}{\partial n_{2}}\\\
\end{array}\right)\Bigg{|}_{P^{4}_{m}}=\left(\begin{array}[]{ccccc}2a_{1}(p^{4}_{m,1}-\frac{p_{1}}{2})&-2u_{2}a_{2}(p^{4}_{m,2}-\frac{p_{2}}{2})\\\
-2u_{1}a_{1}(p^{4}_{m,1}-\frac{p_{1}}{2})&2a_{2}(p^{4}_{m,2}-\frac{p_{2}}{2})\\\
\end{array}\right).$

Thus, the two eigenvalues $\lambda_{m,j}$ ($j=1,2$) of $A|_{P^{4}_{m}}$
satisfy:

$\displaystyle\lambda^{2}_{m,j}-2\left[a_{1}(p^{4}_{m,1}-\frac{p_{1}}{2})+a_{2}(p^{4}_{m,2}-\frac{p_{2}}{2})\right]\lambda_{m,j}+4(1-u_{1}u_{2})a_{1}a_{2}(p^{4}_{m,1}-\frac{p_{1}}{2})(p^{4}_{m,2}-\frac{p_{2}}{2})=0.$
(23)

Note that the discriminant $\Delta_{m}$ for Eq. (23) is:

$\displaystyle\Delta_{m}=4\left[a_{1}(p^{4}_{m,1}-\frac{p_{1}}{2})+a_{2}(p^{4}_{m,2}-\frac{p_{2}}{2})\right]^{2}-16(1-u_{1}u_{2})a_{1}a_{2}(p^{4}_{m,1}-\frac{p_{1}}{2})(p^{4}_{m,2}-\frac{p_{2}}{2}).$

  1. 1.

For $P^{4}_{1}$=($p^{4}_{1,1}$,$p^{4}_{1,2}$), we can get $\Delta_{1}>0$, then
Eq. (23) has two real roots $\lambda_{1,1}$ and $\lambda_{1,2}$. Moreover, we
have:
$\lambda_{1,1}+\lambda_{1,2}=2\left[a_{1}(p^{4}_{1,1}-\frac{p_{1}}{2})+a_{2}(p^{4}_{1,2}-\frac{p_{2}}{2})\right]<0$,
$\lambda_{1,1}\lambda_{1,2}=4(1-u_{1}u_{2})a_{1}a_{2}(p^{4}_{1,1}-\frac{p_{1}}{2})(p^{4}_{1,2}-\frac{p_{2}}{2})>0$,
which implies $\lambda_{1,1}<0,\lambda_{1,2}<0$. Thus, we can conclude the
equilibrium point $P^{4}_{1}$ is a locally stable node.

  2. 2.

For $P^{4}_{2}$=($p^{4}_{2,1}$,$p^{4}_{2,2}$), we can get $\Delta_{2}>0$, then
Eq. (23) has two real roots $\lambda_{2,1}$ and $\lambda_{2,2}$, where
$\lambda_{2,1}<\lambda_{2,2}$. Moreover, we have:
$\lambda_{2,1}\lambda_{2,2}=4(1-u_{1}u_{2})a_{1}a_{2}(P_{21}-\frac{p_{1}}{2})(p^{4}_{2,2}-\frac{p_{2}}{2})<0$,
which implies $\lambda_{2,1}<0<\lambda_{2,2}$. Thus, we have that the
equilibrium point $P^{4}_{2}$ is a saddle point.

  3. 3.

For $P^{4}_{3}$=($p^{4}_{3,1}$,$p^{4}_{3,2}$), we also have $\Delta_{3}>0$,
then Eq. (23) has two real roots $\lambda_{3,1}$ and $\lambda_{3,2}$, where
$\lambda_{3,1}<\lambda_{3,2}$. Moreover, we have:
$\lambda_{3,1}\lambda_{3,2}<0$, which indicates
$\lambda_{3,1}<0<\lambda_{3,2}$. Thus, $P^{4}_{3}$ is a saddle point.

  4. 4.

For $P^{4}_{4}$=($p^{4}_{4,1}$,$p^{4}_{4,2}$), we have $\Delta_{4}>0$,
Moreover, the two real roots $\lambda_{41}$ and $\lambda_{4,2}$ of Eq. (23)
satisfy $\lambda_{4,1}+\lambda_{4,2}>0$, $\lambda_{4,1}\lambda_{4,2}>0$, which
implies $\lambda_{4,1}>0,\lambda_{4,2}>0$. Thus, $P^{4}_{4}$ is an unstable
node.

$\hfill\square$

###  7.2 Global phase portrait derivation

Subsequently, for the four-equilibria system (1) under Condition
($\mathbb{K}^{4}\wedge\mathbb{H}$), the global phase portrait will be derived.
First letting $x_{1}(t)=n_{1}(t)-p^{4}_{1,1}$ and
$x_{2}(t)=n_{2}(t)-p^{4}_{1,2}$, we simplify system (1) as:

$\begin{split}&\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}=-{u_{2}}{a_{2}}(x_{2}-s_{2})^{2}+a_{1}(x_{1}-s_{1})^{2}+M_{1},\\\
&\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}=-{u_{1}}{a_{1}}(x_{1}-s_{1})^{2}+a_{2}(x_{2}-s_{2})^{2}+M_{2},\end{split}$
(24)

where
$s_{1}=\frac{1}{2}\sqrt{{p_{1}}^{2}-\frac{4{d_{1}}+4{u_{2}}{d_{2}}}{{a_{1}}(1-{u_{1}}{u_{2}})}}$,
$s_{2}=\frac{1}{2}\sqrt{{p_{2}}^{2}-\frac{4{d_{2}}+4{u_{1}}{d_{1}}}{{a_{2}}(1-{u_{1}}{u_{2}})}}$,
$M_{1}=d_{1}+u_{2}\frac{{a_{2}p_{2}}^{2}}{4}-\frac{{a_{1}p_{1}}^{2}}{4}$, and
$M_{2}=d_{2}+u_{1}\frac{{a_{1}p_{1}}^{2}}{4}-\frac{{a_{2}p_{2}}^{2}}{4}$.
Obviously, the corresponding four equilibria of system (24) are
$\hat{P}^{4}_{1}=(0,0)$, $\hat{P}^{4}_{2}=(2s_{1},0)$,
$\hat{P}^{4}_{3}=(0,2s_{2})$ and $\hat{P}^{4}_{4}=(2s_{1},2s_{2})$. Note that
the four equilibria of the original system (1) are
${P}^{4}_{1}=(\frac{p_{1}-2s_{1}}{2},\frac{p_{2}-2s_{2}}{2})$,
${P}^{4}_{2}=(\frac{p_{1}+2s_{1}}{2},\frac{p_{2}-2s_{2}}{2})$,
${P}^{4}_{3}=(\frac{p_{1}-2s_{1}}{2},\frac{p_{2}+2s_{2}}{2})$, and
${P}^{4}_{4}=(\frac{p_{1}+2s_{1}}{2},\frac{p_{2}+2s_{2}}{2})$. Moreover, since
$M_{1}+u_{2}M_{2}=d_{1}+u_{2}d_{2}-\left(1-u_{1}u_{2}\right)G_{1,max}$ and
$M_{2}+u_{1}M_{1}=d_{2}+u_{1}d_{1}-\left(1-u_{1}u_{2}\right)G_{2,max}$, we
have that Condition $\mathbb{K}^{4}$ (or inequality (22)) holds if and only if
Condition ($\hat{\mathbb{K}}^{4}$):

$\begin{split}M_{1}+u_{2}M_{2}<0,\\\ M_{2}+u_{1}M_{1}<0,\end{split}$ (25)

holds. Further, Condition ($\hat{\mathbb{K}}^{4}$) can be divided into five
sub-conditions:

  1. (1)

Condition $(\hat{\mathbb{K}}^{4}_{1})$: $M_{1}<0\wedge M_{2}<0$;

  2. (2)

Condition $(\hat{\mathbb{K}}^{4}_{2})$: $M_{1}=0\wedge M_{2}<0$

  3. (3)

Condition $(\hat{\mathbb{K}}^{4}_{3})$: $M_{1}<0\wedge M_{2}=0$

  4. (4)

Condition $(\hat{\mathbb{K}}^{4}_{4})$: $M_{1}>0\wedge
M_{2}<\min\left\\{{-\frac{M_{1}}{u_{2}}},{-u_{1}M_{1}}\right\\}$;

  5. (5)

Condition $(\hat{\mathbb{K}}^{4}_{5})$:
$M_{1}<\min\left\\{{-u_{2}M_{2}},{-\frac{M_{2}}{u_{1}}}\right\\}\wedge
M_{2}>0$

Thus in the following, we will discuss the global qualitative characteristics
and attraction region for equivalence system (24) under Condition
($\hat{\mathbb{K}}^{4}_{i}\wedge\mathbb{H}$) respectively, where
$i=1,2,\cdots,5$.

To begin with, we consider the phase portrait under Condition
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}$), which can be partitioned into
four sub-Conditions: ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$) and
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$).

First consider Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$).
Following the steps in Sec. 3.2, we have that: the demarcation lines are
$\frac{(x_{1}-s_{1})^{2}}{-\frac{M_{1}}{a_{1}}}-\frac{(x_{2}-s_{2})^{2}}{-\frac{M_{1}}{u_{2}a_{2}}}=1$
and
$\frac{(x_{2}-s_{2})^{2}}{-\frac{M_{2}}{a_{2}}}-\frac{(x_{1}-s_{1})^{2}}{-\frac{M_{2}}{u_{1}a_{1}}}=1$,
portrayed as the red and blue dotted curve respectively in Fig. 16(a). In
addition, we have $\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}=0$ and
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}=\left(1-u_{1}u_{2}\right){a_{2}}x_{2}(t)\left({x_{2}(t)-2s_{2}}\right)$
on these two red curves, indicating
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}>0$ when $x_{2}(t)<0$ or
$x_{2}(t)>2s_{2}$ holds and $\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}<0$ when
$0<x_{2}(t)<2s_{2}$ holds. Similarly, we have
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}=0$ and
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}=\left(1-u_{1}u_{2}\right){a_{1}}x_{1}(t)\left({x_{1}(t)-2s_{1}}\right)$
on these two blue curves, implying
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}>0$ when $x_{1}(t)<0$ or
$x_{1}(t)>2s_{1}$ holds and $\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}<0$ when
$0<x_{1}(t)<2s_{1}$ holds. Obviously, we can obtain nine unbounded regions in
$\mathbb{R}^{2}$ divided by these four curves. Moreover, Table 3 summarizes
the symbols of $\frac{\mathrm{d}{x}_{i}(t)}{\mathrm{d}t}$ ($i=1,2$) in these
nine regions. On the basis of Table 3, the phase portrait of (24) under
Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$) can be obtained, as
shown in Fig. 16(a). Subsequently, based on the phase portrait 5(a), by
utilizing proof by contradiction, we can verify that regions B, G and F are
positive invariant sets. Since the trajectories starting from B, G and F will
not escape, there exist no close orbits for system (24) under Condition
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$).

(a) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$

(b) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$

(c) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$

(d) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$

Figure 16: Theoretical schematic phase portraits of system (24) under Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}$). Green circle: desired uncongested stable equilibrium. Yellow circle: saddle equilibrium. Yellow-brown circle: unstable equilibrium. Red dotted curves: $\dot{x}_{1}=0$, blue dotted curves: $\dot{x}_{2}=0$, green arrows: trajectory directions.  Table 3: Symbols of the nine regions for system (24) under Condition $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$. Derivative | Symbol in the region
---|---
A | B | C | D | E | F | G | H | I
$\frac{\mathrm{d}{x}_{1}(t)}{\mathrm{d}t}$ | $+$ | $-$ | $+$ | $+$ | $-$ | $+$ | $+$ | $-$ | $+$
$\frac{\mathrm{d}{x}_{2}(t)}{\mathrm{d}t}$ | $+$ | $+$ | $+$ | $-$ | $-$ | $-$ | $+$ | $+$ | $+$

Then, for Conditions ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$) and
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$), we can also simplify the
system (24) and then perform the above analysis. In a similar way, we can
obtain four corresponding curves under each Condition, which separate
$\mathbb{R}^{2}$ to nine unbounded regions. Then the corresponding phase
portraits of (24) under these three Conditions can be acquired, as shown in
Figs. 16(b), 16(c) and 16(d). Further, we can also obtain the positive
invariant sets under each Condition, which verify that there is no close orbit
for system (24) under Conditions
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$) and
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$). Thus, there exists no close
orbit, for the system (1) under Condition $\hat{\mathbb{K}}^{4}_{1}\wedge H$.

###  7.3 Spontaneous attraction region estimation

Based on the nonexistence of close orbits and Fig. 16, we will next analyze
the attraction region of the locally stable node $P^{4}_{1}$. For this
purpose, we wish to find certain corresponding separatrices of system (24)
under Conditions ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$) and
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$). However, due to the
increased number of equilibria, it seems very difficult to find certain
separatrices directly. Thus, we will first find inner and outer estimations of
attraction regions, and then find certain corresponding separatrices under the
above four conditions.

To start with, we consider Condition
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$). The two red dotted
demarcation lines can be denoted as $x_{1}=\hat{r}^{4,1}_{i}(x_{2})$
($i=1,2$), where
$\hat{r}^{4,1}_{i}(x_{2})=(-1)^{i}\sqrt{\frac{u_{2}a_{2}}{a_{1}}\left(x_{2}-s_{2}\right)^{2}-\frac{M_{1}}{a_{1}}}+s_{1}$.
Similarly, the two blue dotted demarcation lines can be denoted as
$x_{2}=\hat{b}^{4}_{i}(x_{1})$ ($i=1,2$), where
$\hat{b}^{4}_{i}(x_{1})=(-1)^{i}\sqrt{\frac{u_{1}a_{1}}{a_{2}}\left(x_{1}-s_{1}\right)^{2}-\frac{M_{2}}{a_{2}}}+s_{2}$.
Obviously, the vertex of the line $x_{1}=\hat{r}^{4,1}_{i}(x_{2})$ is
$Ve_{i}=((-1)^{i}\sqrt{-\frac{M_{1}}{a_{1}}}+s_{1},s_{2})$ and the vertex of
the line $x_{2}=\hat{b}^{4}_{i}(x_{1})$ ($i=1,2$) is
$Ve_{i+2}=(s_{1},(-1)^{i+2}\sqrt{-\frac{M_{2}}{a_{2}}}+s_{2})$. Subsequently,
in light of Fig. 16(a), we can similarly verify that the region $\hat{U}^{4}$
(the light yellow region shown in Fig. 17(a)) is a positive invariant set,
where $\hat{U}^{4}=\left\\{(x_{1},x_{2})|x_{1}\geq\hat{l}^{2}_{out}(x_{2})\vee
x_{2}\geq\hat{l}^{1}_{out}(x_{1})\right\\}$ with

$\hat{l}^{i}_{out}(x_{i})=\begin{cases}\sqrt{\frac{u_{i}a_{i}}{a_{j}}\left(x_{i}-s_{i}\right)^{2}-\frac{M_{j}}{a_{j}}}+s_{j}&\text{$x_{i}<0$}\\\
2s_{j}&\text{$0<x_{i}<2s_{i}$}\end{cases}$

Moreover, by similar analysis to $\hat{U}^{2a}$ in two-equilibria case, we can
arrive that
$\mathcal{R}(\hat{P}^{4}_{1})\subset\mathbb{R}^{2}\setminus\hat{U}^{4}$, i.e.,
$\hat{S}^{4}_{out}=\mathbb{R}^{2}\setminus\hat{U}^{4}=\left\\{(x_{1},x_{2})|x_{1}<\hat{l}^{2}_{out}(x_{2}),~{}x_{2}<\hat{l}^{1}_{out}(x_{1})\right\\}$
is an outer estimation of attraction region for $\hat{P}^{4}_{1}$.

In addition, $\hat{S}^{4}_{in}$ (the light blue region shown in Fig. 17(a)) is
also a positive invariant set, where

$\hat{S}^{4}_{in}=\left\\{(x_{1},x_{2})|x_{1}\leq\hat{l}^{2}_{in}(x_{2}),~{}x_{2}\leq\hat{l}^{1}_{in}(x_{1})\right\\}$

with

$\hat{l}^{i}_{in}(x_{i})=\begin{cases}2s_{j}&\text{$x_{i}<0$}\\\
\sqrt{\frac{u_{i}a_{i}}{a_{j}}\left(x_{i}-s_{i}\right)^{2}-\frac{M_{j}}{a_{j}}}+s_{j}&\text{$0<x_{i}<s_{i}$}\\\
\sqrt{-\frac{M_{j}}{a_{j}}}+s_{j}&\text{$s_{i}\leq n_{2}<2s_{i}$}\end{cases}$

Further, we can prove that $\hat{S}^{4}_{in}$ is an inner estimation of
attraction region for $\hat{P}^{4}_{1}$ by following two steps. Firstly, the
trajectory starting from any point in $\hat{S}^{4}_{in}$ will entering the
bounded region $\hat{S}^{4}_{in}\wedge\mathbb{R}^{2}_{1}$ as
$t\rightarrow+\infty$, where

$\mathbb{R}^{2}_{1}=\left\\{(x_{1},x_{2})|x_{1}\geq-\frac{p_{1}-2s_{1}}{2},~{}x_{2}\geq-\frac{p_{2}-2s_{2}}{2}\right\\},$

as we have $\dot{x}_{1}>0$ on the left boundary of
$\hat{S}^{4}_{in}\wedge\mathbb{R}^{2}_{1}$ and $\dot{x}_{2}>0$ on the lower
boundary of $\hat{S}^{4}_{in}\wedge\mathbb{R}^{2}_{1}$. Secondly, the
trajectory starting from any point in bounded region
$\hat{S}^{4}_{in}\wedge\mathbb{R}^{2}_{1}$ will go to $\hat{P}^{4}_{1}$ as
$t\rightarrow+\infty$, since there is no closed orbit, and there exists only
one stable node $\hat{P}^{4}_{1}$ in the bounded region (Zhang (2006)).

(a) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$

(b) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$

(c) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$

(d) $\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$

Figure 17: Qualitative characteristics of system (24) ((a)-(d)) and system (1)
((e)-(h)) under Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}$). Green
circle: desired uncongested stable equilibrium. Yellow circle: saddle
equilibrium. Yellow-brown circle: unstable equilibrium. Red dotted curves:
$\dot{x}_{1}=0$, blue dotted curves: $\dot{x}_{2}=0$, yellow curves: the
separatrices, also serving as boundaries of attraction region. Light pink
region: attraction region of $P^{4}_{1}$.

As for Conditions ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$),
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{3}$) and
($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{4}$), through similar analysis, we
can obtain the inner and outer estimations of attraction regions for
$\hat{P}^{4}_{1}$, as shown in the Figs. 17(b), 17(c) and 17(d). Further,
letting $n_{1}(t)=x_{1}(t)+\frac{p_{1}-2s_{1}}{2}$ and
$n_{2}(t)=x_{2}(t)+\frac{p_{2}-2s_{2}}{2}$, the inner and outer estimations of
attraction regions for $P^{4}_{1}$ can be obtained. Denoting
$l^{i}_{in}(n_{i})=\hat{l}^{i}_{in}(n_{i}-\frac{p_{i}-2s_{i}}{2})+\frac{p_{j}-2s_{j}}{2}$
and
$l^{i}_{out}(n_{i})=\hat{l}^{i}_{out}(n_{i}-\frac{p_{i}-2s_{i}}{2})+\frac{p_{j}-2s_{j}}{2}$,
where $i,j=1,2$ and $i\neq j$, we can arrive at the following theorem.

###### Theorem 7.2

Under Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}$), for the local
stable point $P^{4}_{1}$, we have
$S^{4}_{in}\subset\mathcal{R}(P^{4}_{1})\subset S^{4}_{out}$, where the inner
estimation of attraction region $S^{4}_{in}$ and the outer estimation of
attraction region $S^{4}_{out}$ are defined as follows:

  1. (1)

Under Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{1}$),

$S^{4}_{in}=\left\\{(n_{1},n_{2})|n_{1}\leq l^{2}_{in}(n_{2}),~{}n_{2}\leq
l^{1}_{in}(n_{1})\right\\};S^{4}_{out}=\left\\{(n_{1},n_{2})|n_{1}\leq
l^{2}_{out}(n_{2}),~{}n_{2}\leq l^{1}_{out}(n_{1})\right\\}.$

  2. (2)

Under Condition ($\hat{\mathbb{K}}^{4}_{1}\wedge\mathbb{H}_{2}$),