Analyzing the Impact of Vandalism, Hoarding, and Strikes on Fuel Distribution in Nigeria

1. Introduction

Petroleum products, the lifeblood of Nigeria’s economy, are critical for powering industries, transportation, and households across the nation. However, the persistent scarcity of these refined products has posed significant challenges, hampering economic growth and adversely affecting the general populace. A deeper understanding of the distribution system is essential for formulating effective strategies to eradicate scarcity and ensure a consistent and reliable supply of petroleum products.

Crude oil exploration in Nigeria dates back to 1908, when the Shell Darcy Development Petroleum Company commenced operations in Okitipupa and Araromi, in present-day Ondo State [1,2,3,4]. These efforts continued until 1938 [4]. In 1956, Shell-BP discovered crude oil in commercial quantities in Oloibiri, Rivers State [2]. This landmark discovery positioned Nigeria among the ranks of oil-producing nations, with its first oil field going operational in 1958, producing 5,100 barrels per day [5,6,7]. Since then, oil production has grown significantly [8].

Nigeria boasts four refineries with a combined installed production capacity of 445,000 barrels per day [9,10,11], initially adequate to meet domestic consumption and allow for exports. Recently, six additional refineries with a total refining capacity of 41,000 barrels per day were established between 2020 and 2022 [12]. Despite this infrastructure and Nigeria’s proven oil reserves [1,3], the country has yet to fully harness its potential to meet domestic demand for refined products.

Paradoxically, while Nigeria is Africa’s largest crude oil producer and one of the continent’s largest exporters [11,12], it remains plagued by recurrent shortages of refined petroleum products. The country imports the majority of its refined products (e.g., kerosene, PMS, aviation kerosene, and diesel) to meet domestic needs [13]. Overreliance on imported refined products, coupled with an aging and inefficient distribution infrastructure, has contributed to chronic scarcity. This scarcity disrupts economic activities, increases operational costs, and fuels public dissatisfaction.

Efforts to address the problem, which spans from military administrations to the current democratic dispensation, have not yielded lasting solutions [14]. Nigeria remains one of the few oil producing countries that imports 95% of its domestic petroleum products [5,12,15]. The smuggling of refined products, particularly PMS, to neighboring countries where pump prices are significantly higher, exacerbates the situation. Fuel-loaded trucks intended for Nigerian service stations are often diverted to these border towns by smugglers to make more money [16].

To address these systemic issues, the Nigerian government recently deregulated the petroleum sector by removing fuel (PMS) subsidies and enacting the Petroleum Industry Act (PIA) [17]. These measures aim to minimize, if not eradicate fuel scarcity and attract private sector investment [2,10,13,18,19]. However, challenges such as pipeline vandalism, road transportation inefficiencies, hoarding, smuggling, and Nigeria’s growing population continue to strain the downstream sector.

The Dangote Refinery, a multibillion-dollar project nearing operational commencement, is widely anticipated to transform Nigeria’s fuel supply landscape. Many citizens view the refinery as a potential game-changer for addressing the nation’s domestic fuel demand. Despite significant government investment, concerns persist regarding the state of supporting infrastructure. Poor road conditions around the refinery pose serious risks to safety and logistics. Moreover, the Petroleum Industry Act (PIA) lacks provisions for road infrastructure improvements or alternative transportation routes, leaving critical gaps in planning and preparedness.

This study aims to address these challenges by examining the historical and current state of Nigeria’s downstream sector, with a focus on issues such as pipeline vandalism, industrial actions by oil and gas workers, and the hoarding and diversion of petroleum products. It provides insights from a novel mathematical model designed to guide decision-makers in managing the downstream sector more effectively. Additionally, the study seeks to demystify the operations of the NNPC’s downstream activities for public understanding while emphasizing the need for sustainable and effective policies to benefit all stakeholders.

The remainder of the manuscript is organized as follows: a glossary of terms and abbreviations clarifies technical terminology in Table 1. Section 2 describes the model and its formulation. Section 3 explores the mathematical properties of the model and presents analytical results. Numerical simulations and their interpretations are detailed in Section 4 . Finally, Section 5 summarizes the key findings, draws conclusions, and outlines the study’s limitations.

2. Scientific Approach and Model Formulation

This section outlines the scientific methodology employed, reviews related studies, and presents the formulation of the proposed model.

2.1. Scientific Approach, Literature Review, and Justification

Mathematical models provide a systematic approach to representing real-world systems, enabling critical analysis and informed decision-making [20,21,22,23]. According to Rutherford Aris: "a set of mathematical equations that provide an adequate description of a physical system."

Mathematical modeling has been extensively applied in epidemiology to study population dynamics and disease transmission [24,25]. Its applications have expanded into social sciences, including economic and energy systems [26,27,28,29,30], with studies emphasizing global stability analysis and compartmental frameworks for complex systems [31,32].

Modeling has also been instrumental in optimizing energy systems. For example, a bilevel mixed-integer nonlinear programming model integrates supply chain design with competitive market dynamics to optimize fuel pricing [33]. Gasification models for solid fuels incorporate temperature, pressure, and fuel composition to enhance syngas production [34,35]. Additionally, research on fuel mixture ignition highlights the impact of droplet size distribution on combustion efficiency [36], while graph theory and differential equations improve operational planning for underground gas storage [37].

Despite these advancements, studies on Nigeria’s downstream petroleum sector primarily focus on operational challenges rather than intrinsic efficiency and effectiveness, which are crucial for ensuring value for money. Key issues include pipeline vandalism, governance inefficiencies, and supply disruptions [14,38,39,40,41,42]. Weak regulatory frameworks hinder coordination among stakeholders [9], while transportation and storage losses further exacerbate shortages [9,43,44].

2.2. Justification for the Compartmental Model

A compartmental model is well-suited for analyzing Nigeria’s petroleum sector due to its complexity, as highlighted by the Organization of Economic Cooperation (OEC) [15]. These models divide systems into manageable sub-populations governed by ordinary differential equations (ODEs), allowing for quantitative analysis and prediction [45,46]. Widely used in epidemiology, biology, and social sciences, they offer a structured way to analyze system dynamics, integrate diverse data, and derive analytical solutions [47,48,49].

Unlike Agent-Based Models (ABMs), which simulate individual interactions, or statistical models, which focus on correlations, compartmental models emphasize population-level dynamics and mechanistic relationships [50,51]. While they involve simplifying assumptions and challenges in parameter estimation, they remain effective for studying nonlinear systems [46,48].

Fuel scarcity in Nigeria affects nearly all citizens, with inflationary consequences sparing only a few elites and top officials. Given the sector’s complexity and the limited availability of detailed data, a compartmental approach provides a practical and efficient framework for studying fuel distribution challenges. By capturing key factors influencing supply and disruptions, this model offers a structured approach for evaluating policy interventions and improving system resilience.

2.3. A Guide to Model Formulation

The operational framework of the Petroleum Products Marketing Company (PPMC), as outlined on the Nigerian National Petroleum Corporation (NNPC) webpage [52], serves as the primary guide for this model formulation, supported by additional resources such as [2,10,17,38,53,54]. The petroleum distribution network in Nigeria involves pipelines connecting 21 white product depots and tank farms nationwide [43,55]. However, frequent acts of vandalism have led the Nigerian government and NNPC management to restructure the old PPMC into the Nigerian Pipelines and Storage Company (NPSC), PPMC, and NNPC Shipping [3,52]. This model is based on that restructuring, aiming to eliminate product shortages at the consumer level.

2.4. Model Formulation

The proposed model captures the flow of petroleum products across key nodes in Nigeria’s downstream sector. The system consists of six compartments:

• Product Refinement/Importation $\left(Q\right(t\left)\right)$: Initial supply of petroleum products
• Storage Facilities ($A\left(t\right)$): Large-scale holding facilities managed by refineries or independent marketers (124 storage depots in 21 locations nationwide)[55,56].
• Transporter/NUPENG/PENGASSAN/PPMC: Distributors responsible for delivering products to retail stations.
• Retail Stations $(E_1\left(t\right),E_2\left(t\right))$: Points of sale, classified into major marketers (NNPC Retail/Mega Stations) and independent marketers.
• Consumers ($M\left(t\right)$): The final stage where petroleum products are utilized.

The interactions among these compartments are modeled using ordinary differential equations (ODEs), incorporating key operational, policy, and environmental factors such as vandalism rates ($V\left(t\right)$), natural loss rates ($\mu$), and stakeholders contributions ($x,y$).

The compartments, represented as ($QAZ{E}_1{E}_{2}M$) evolve over time t. At time t, the NNPC downstream infrastructure for refined petroleum products, represented by $Q\left(t\right)$ holds a quantity of petroleum product, replenished daily at a rate $\mathrm{\Lambda }$ (representing the amount refined or imported by PPMC). The proportion of products transported via pipelines per day is denoted by $\alpha$, while vandalism, modeled as a constant function ($V\left(t\right)$), results in a loss of $\alpha V\left(t\right)$ ([39,40,41,42,57,58,59,60,61]).

Losses occur at various stages; Depot losses, $\mu$, Retail station loses ${\mu }_1,$ and consumer-level loses (${\mu }_2$). The depot losses are documented in NNPC’s 2021 Monthly FAAC reports and illustrated in Figure 1.

The quantity of fuel dispensed to retail marketers depends on the refined products from NNPC local refineries and PPMC [9]. The percentage distributed from NNPC retail depots via PENGASSAN/NUPENG to filling stations is denoted by c, resulting in a net distribution rate of $cA\left(t\right)Z\left(t\right)$.

To facilitate the model formulation, the following fundamental assumptions are considered:

2.4.1. Basic Assumption of the Model

The model assumptions are as follows:

• Multiple factors contribute to petroleum product scarcity in Nigeria, including NUPENG strikes (Z(t)), fuel diversion, hoarding ([2,62]), or persistent vandalization of NNPC pipelines.
• Vandalization is treated as a constant function ($V\left(t\right)=k$) and satisfied one of the properties of probability $0<V\left(t\right)\le 1$ ([39,40,41,57,58,59]), ([61,63]) except in case where $V\left(t\right)$ is zero.
• Literature highlights that vandalization increases fuel distribution costs ([14,43,64]), resulting in scarcity for consumers. An alternative route for transportation is incorporated in the model.
• Fuel scarcity leads to price increases, either directly at the pump or through parallel markets (black markets), prompting hoarding and diversion.
• Ideally, retail marketers are meant to be the sole suppliers to consumers unless tanker drivers siphon products, so a siphon parameter ($\eta$) is included in the model.
• All parameters are non-negative; any negative results in the model will be interpreted as impractical.
• Fuel scarcity can be resolved when consumer demand is lower than supply from major and independent marketers. Therefore, policy and good implementation are embedded in the model simulation

According to Ambituuni et al., fuel distribution in Nigeria’s downstream sector is shared among NNPC Retail Ltd (23.43%), six major marketers (25.47%), and independent marketers (51%). Thus, the combined share of major marketers and NNPC Retail is x, while independent marketers contribute y. The distribution formula allocates 48.90% to NNPC retail stations and major marketers, with independent marketers receiving 51%. A siphon factor $\eta$ and natural losses at retail stations ${\mu }_1$ account for the remaining 0.1%. Despite differing infrastructures [56], major marketers and NNPC Retail are treated as a single entity due to their similar pricing regulations and operations. The rate of product distribution to consumers is expressed as:

$${\displaystyle \frac{dZ}{dt}}=cAZ-(x+y+\eta +\mu )Z$$

(1)

In cases of severe vandalism, an alternative railway transport route is introduced, modeled as $b_{1}Q\left(t\right)$ to offset losses. Since rail transport is independent of road infrastructure (Ambituuni et al.), the adjusted fuel distribution dynamics during crises are:

$${\displaystyle \frac{dZ}{dt}}=b_{1}Q+cAZ-(b_2+x+y+\eta +\mu )Z$$

(2)

The product volume sold to consumers by major and independent marketers, denoted by $k_1$ and $k_2$, respectively, is modeled as:

$$\begin{array}{cc}& {\displaystyle \frac{d{E}_1}{dt}}=xZ-(k_1+\mu )E_1\hfill \end{array}$$

(3)

$$\begin{array}{c}\hfill {\displaystyle \frac{d{E}_2}{dt}}=yZ-(k_2+\mu )E_2\end{array}$$

(4)

2.4.2. General Form of the Model

The model formulation, along with the assumptions, is visually summarized in the diagram (Figure 2) below, followed by the corresponding mathematical expressions:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dQ}{dt}}=\mathrm{\Lambda }-(a(1-V)+b_1+\mu )Q\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dA}{dt}}=a(1-V)Q+b_{2}Z-(\mu +cZ)A\hfill \end{array}$$

(6)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dZ}{dt}}=b_{1}Q+cAZ-(b_2+\mu +x+y+\eta )Z\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{E}_1}{dt}}=xZ-({\mu }_1+k_1)E_1\hfill \end{array}$$

(8)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{d{E}_2}{dt}}=yZ-({\mu }_1+k_2)E_2\hfill \end{array}$$

(9)

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dM}{dt}}=-{\mu }_{2}M+k_1{E}_1+k_2{E}_2+\eta Z\hfill \end{array}$$

(10)

Additionally, the total population (i.e. the amount (quantity) of refined products) at time t is represented by the equation: $N\left(t\right)=Q\left(t\right)+A\left(t\right)+Z\left(t\right)+E_1\left(t\right)+E_2\left(t\right)+M\left(t\right)$

This formulation captures the flow of products, losses, and interactions between various subsystems, providing a framework for exploring the effects of policy changes and operational interventions. The description of the model parameters are presented in Table 2 Further analysis, including stability and sensitivity, follows in subsequent sections.

3. Properties and Analysis of the Mathematical Model

3.1. Model Properties

Since the model describes fuel distribution, where fuel quantities cannot be negative, all model parameters and state variables are assumed to be non-negative for any time t. This assumption is formalized in the following theorem:

Theorem 1. 

Let:

• ${\mathrm{\Omega }}_1\in {\mathbb{R}}_{+}^3$ represents a set of storage facilities owned by NNPC/MOMAN for distributing products.
• ${\mathrm{\Omega }}_2\in {\mathbb{R}}_{+}^2$ represents the storage at retail (refilling) stations.
• $M\left(t\right)$ represents the quantity of fuel stored by consumers, and
• Ω is the combined quantity of refined products in the country.

Then

$$\mathrm{\Omega }=\left\{Q\left(t\right),A\left(t\right),Z\left(t\right),E_1\left(t\right),E_2\left(t\right),M\left(t\right)\in {\mathbb{R}}_{+}^6:{\mathrm{\Omega }}_1+{\mathrm{\Omega }}_2+M\left(t\right)=N\left(t\right)\mathrm{and}0<N\left(t\right)\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\right\}$$

(11)

bounded and it is non-negative for all $t\in [0,t_f]$ and at $t=0,N\left(0\right)=N_0$. Hence, Ω is positively invariant and attracting for the proposed model.

Proof. 

The total quantity of refined product in circulation is denoted by $N\left(t\right)$, summing (5) to (10)

$$\begin{array}{cc}\hfill {\displaystyle \frac{dN}{dt}}& ={\displaystyle \frac{dQ}{dt}}+{\displaystyle \frac{dA}{dt}}+{\displaystyle \frac{dZ}{dt}}+{\displaystyle \frac{d{E}_1}{dt}}+{\displaystyle \frac{d{E}_2}{dt}}+{\displaystyle \frac{dM}{dt}}\hfill \\ \hfill {\displaystyle \frac{dN}{dt}}& =\mathrm{\Lambda }-\mu N\left(t\right)\hfill \end{array}$$

The integral solution of this differential equation is

$$N\left(t\right)=N_0{e}^{-\mu t}+{\displaystyle \frac{\mathrm{\Lambda }}{\mu }}\left(1-e^{-\mu t}\right)$$

Therefore,

$$0<\underset{t\to \infty }{lim}supN\left(t\right)\le {\displaystyle \frac{\mathrm{\Lambda }}{\mu }}$$

This shows that $N\left(t\right)$ is ultimately bounded for all t.

Next, to prove the non-negativity of the state variables, consider

$${\displaystyle \frac{dQ}{dt}}=\mathrm{\Lambda }-(a(1-V)+b_1+\mu )Q$$

$${\displaystyle \frac{dQ}{dt}}+(a(1-V)+b_1+\mu )Q=\mathrm{\Lambda }$$

Using integrating factors

$$I.F.=e^{(a(1-V)+b_1+\mu )t}$$

the solution for $Q\left(t\right)$ becomes:

$$Q\left(t\right)=Q\left(0\right)e^{(a(1-V)+b_1+\mu )t}$$

(12)

Since $Q\left(t\right)\ge 0\forall t\ge 0$, a similar approach applied to the other state variables shows that they are all non-negative, completing the proof. □

3.2. Existence and Uniqueness of Solutions

To establish that the model is mathematically well-posed, we verify the Lipschitz condition. For clarity, we summarize this in the following theorem and provide its interpretation. Detailed mathematical proofs are available in Appendix A.1.

Theorem 2. 

The system (5)–(10) has a unique solution, provided the functions are continuously differentiable in $\mathrm{\Omega }\in {\mathbb{R}}^6$, as defined in Equation (11), and the system satisfies the Lipschitz condition.

Interpretation: Continuous differentiability and the Lipschitz condition ensure that the model is mathematically well-posed, meaning the model reliably describes Nigeria’s fuel crisis dynamics. This guarantees its applicability in managing NNPCL’s downstream operations and supports decision-making for policymakers.

The next section will explore the model’s stability under various scenarios.

4. Model Analysis: Equilibrium Points, Threshold Parameter, and Stability Analysis

This section examines key scenarios affecting petroleum product distribution in oil-producing countries, particularly Nigeria, highlighting factors that may disrupt supply and potential mitigation strategies.

• High Pipeline Vandalism:
  • Significant vandalism severely disrupts product distribution,
  • The system relies solely on pipelines, with no secured alternative distribution routes.
  • Despite the stability of other system components, high vandalism creates a major bottleneck.
• No Pipeline Vandalism: Assumes no vandalism, ensuring uninterrupted pipeline distribution. However, industrial actions may still disrupt supply despite functional infrastructure.
• Mitigation Measures and Alternative Routes:
  • Policies and interventions are implemented to curb vandalism.
  • Alternative distribution routes (e.g., railways for inter-state transportation and long distance transport) reduce dependence on pipelines.
  • A diversified distribution strategy enhances system resilience against disruptions.

These scenarios underscore the need for proactive measures to secure fuel distribution, addressing vandalism, industrial actions, and security threats to oil and gas workers. Evaluating these factors helps strengthen the supply chain and ensure a more reliable distribution system. The following subsection presents the mathematical analysis of these scenarios.

4.1. POINTS OF EQUILIBRIUM

Theorem 3. 

The system of Equations (5)–(10) exhibits five equilibrium points, denoted as $\mathrm{\Gamma }=E_o^i,E_o^{ii},E_o^{iii},E_o^{iv},E_o^v$, subject to the following conditions:

• $0\le V\left(t\right)\le 1$ (vandalism level),
• $0\le Z\left(t\right)<\infty$ (industrial action level),

implies vandalism reaches its maximum level when $V=1$, and NUPENG embarks on industrial action (i.e., $Z=0$).

The proof of this theorem is in the Appendix A.2. However, we present the summary of the prove below for continual readability.

To establish Theorem 3, the right-hand sides (RHS) of Equations (5)–(10) are set to zero.

The equilibrium points, denoted as $\mathrm{\Gamma }=\left\{E_0^i,E_0^{ii},E_0^{iii},E_0^{iv},E_0^v\right\}$ are obtained under various scenarios outlined at the beginning of this Section 4 and the elements of the set $\mathrm{\Gamma }$ are respectively defined as:

• Critical Collapse - Severe fuel scarcity with widespread economic and social disruption.
• Partial Collapse - Significant shortages are causing moderate disruptions to daily activities.
• Partial Functionality - Limited fuel availability with manageable impacts.
• Moderate Functionality - Consistent but suboptimal fuel supply meeting most consumer needs.
• Optimal Distribution - Smooth and efficient fuel supply, fully satisfying demand.

Interpretation:

• Critical Collapse ($E_0^i$): When vandalism reaches its peak and industrial actions (e.g., strikes by key transporters or workers such as NUPENG) occur, the entire fuel distribution system collapses. In this state, fuel is confined to the White Tank storage with no distribution across the supply chain. This scenario mirrors recent events where industrial actions triggered nationwide fuel crises [65].
• Partial Collapse ($E_0^{ii}$): In this scenario, vandalism remains at maximum levels $(V\approx 1)$, and no alternative transportation routes exist $(b_1=0,b_2=0)$, but transporter activities continue $(Z\ne 0)$. This state may give the illusion of operational functionality while the system struggles to maintain supply. Further details are provided in Equation (A10).
• Partial Functionality ($E_0^{iii}$): This event occurs when vandalism is significantly reduced, but the absence of transporters $(Z=0)$ disrupts distribution. Although the system is less strained, the lack of transporter involvement continues to impair functionality.
• Moderate Functionality ($E_0^{iv}$): In this case, vandalism is minimal $(V\ll 1)$, and while transporters are active $(Z\ne 0)$, no alternative transportation routes are available $(b_1=0,b_2=0)$. However, this scenario does not account for the anticipated population growth, which will inevitably lead to increased fuel demand.
• Optimal Functionality ($E_0^v$): This scenario represents a state of crisis preparedness and ensures optimal operation at all times. It is characterized by minimal vandalism, active transporters, and alternative routes that provide resilience against disruptions.

More details are provided in Appendix (Appendix A.2)

4.2. THRESHOLD FOR PRODUCT DISTRIBUTION

Nigeria’s growing population [66] without proportional infrastructure growth raises concerns. To address fuel scarcity, the Nigerian National Petroleum Corporation (NNPC) must proactively analyze demand and, where necessary, improve the supply chain. Collaboration with other stakeholders to resolve bottlenecks in the downstream sector is essential for hitch-free fuel supply. The daily consumption (threshold) to prevent PMS scarcity is derived using the Next Generation Matrix (NGM) approach in [47], the parameters that constituted this threshold (including $\mathrm{\Lambda }$) are estimated from macroeconomic trends ([67]) and NNPC data ([39,40,41,42]).

While Figure 3 illustrates Nigeria’s population growth and annual percentage change, reports from Premium Times and other Nigerian news outlets indicate a decline in petrol (PMS) consumption following the removal of the fuel subsidy by the current Nigerian government. This trend is depicted in Figure 4. However, the actual daily fuel consumption before the subsidy removal on May 29, 2023, remains uncertain.

This subsection analytically derives a threshold parameter for daily fuel distribution to meet consumption needs and numerically estimates it using available data from [39,40,41,42]. The necessary conditions to eliminate fuel scarcity are summarized in the following theorem:

Theorem 4. 

Fuel scarcity can be eliminated if threshold parameters $R_0$ and $R_0^a$ satisfy:

• $R_0>1$, $R_0^a>1$, in the event of a pipeline crisis.
• $R_0\approx R_0^a\ge 1$ under smooth transportation conditions.

The interested reader will find the proof of this Theorem 4 at Appendix A.3, but the interpretation is:

Interpretation:

Equation (A16) represents the analytic threshold derived from the proposed model, capturing two scenarios: transportation without an alternative route and with an alternative route.

Equation (A16) defines the analytic threshold derived from the proposed model, capturing both transportation scenarios with and without alternative routes.

This threshold is critical for addressing corruption in Nigeria’s downstream oil and gas sector. Despite NNPC Limited and its partners being responsible for fuel distribution, they have not provided a clear or reliable estimate of the country’s actual daily PMS consumption. Investigative journalism has played a crucial role in exposing discrepancies in reported figures. Therefore, deriving this threshold through a scientific approach is essential.

In Section 5, we will use data from FAAC reports, NMDPRA, and other available sources earlier mentioned to estimate the threshold value numerically.

4.3. STABILITY ANALYSIS OF THE MODEL

This section establishes the criteria for evaluating the viability of the proposed model in decision-making. Stability analysis, a fundamental tool in mathematical biology [47], can be performed locally around equilibrium points or globally across the entire system. Here, we examine both local (detail in Appendix A.4) and global stability (Appendix A.5) of the equilibrium points identified in Theorem 3.

4.3.1. Methodology for Stability Analysis

To determine the stability at each equilibrium point, we employ the following approach:

• Linearization: Linearize the system Equations (5) to (10) around the equilibrium points by computing the Jacobian matrix,
• Eigenvalue Computation: Evaluate the eigenvalues of the Jacobian matrix,
• Stability Characterization: Classify the equilibrium’s stability based on the eigenvalues, using the established criteria in [68], summarized in Table 3,
• Conditions for Stability: Derive conditions on system parameters that ensure the eigenvalues meet the stability criteria.

Table 3. Characterization of Stability Analysis. [68]

Eigenvalues of J	Type of Critical Point	Remark
${\lambda }_1\ge {\lambda }_2>0$	Source	Unstable
${\lambda }_1\le {\lambda }_2<0$	Sink	Asymptotically Stable
${\lambda }_1>{\lambda }_2>0$	Saddle Point	Unstable
${\lambda }_1,{\lambda }_2=r\pm i\mu ,r>0$	Outward Spiral	Unstable
${\lambda }_1,{\lambda }_2=r\pm i\mu ,r<0$	Inward Spiral	Asymptotically Stable
${\lambda }_1$,${\lambda }_2=\pm i\mu$	Centre	Indeterminate

Table 3. Characterization of Stability Analysis. [68]

Eigenvalues of J	Type of Critical Point	Remark
${\lambda }_1\ge {\lambda }_2>0$	Source	Unstable
${\lambda }_1\le {\lambda }_2<0$	Sink	Asymptotically Stable
${\lambda }_1>{\lambda }_2>0$	Saddle Point	Unstable
${\lambda }_1,{\lambda }_2=r\pm i\mu ,r>0$	Outward Spiral	Unstable
${\lambda }_1,{\lambda }_2=r\pm i\mu ,r<0$	Inward Spiral	Asymptotically Stable
${\lambda }_1$,${\lambda }_2=\pm i\mu$	Centre	Indeterminate

4.4. Stability Result Summary

The stability results are summarized in the below theorem and the proof in Appendix A.4.1:

Theorem 5. 

The system (5)–(10) is locally asymptotically stable (LAS) near Scenario II if the following conditions hold:

• $\alpha V<\alpha +\mu +b_1$,
• $R_0>1$ and $R_0^a>1$, ensuring product demand does not exceed supply, regardless of the distribution method.

Additionally,

• $E_0^{iii}$ is globally asymptotically stable $R_0<1$
• $E_0^{iv}$ is globally asymptotically stable whenever the arithmetic mean exceeds the geometric mean of Lyapnouv function for $R_0>1$

4.4.1. Interpretation

The results suggest that fuel scarcity can be eliminated if government policies ensure a balanced operational framework among stakeholders. Aligning policy interventions with the model’s stability conditions will enhance fuel supply consistency and resilience against disruptions.

4.4.2. Terminology Consistency

To maintain coherence with established terminology in mathematical epidemiology, the following conventions are adopted:

• Scenario IIa is referred to as the Products Distribution-Free Equilibrium (PDE), representing conditions where no products are distributed.
• Scenario IIb is referred to as the Products Distribution Persistent Equilibrium (PDP), representing conditions where product distribution continues over time.
• The parameter $R_e$ is introduced as a generalized threshold encompassing both $R_0$ and $R_0^a$, serving as a comprehensive indicator of system stability.
• Analytical results use conditions $R_0>1$ or $R_0<1$ to establish equilibrium stability, following classical approaches.
• Numerical simulations compute actual $R_e$ values to provide a realistic representation of fuel distribution dynamics.

4.4.3. Appendix Reference

The mathematical proof of Theorem 5 is provided in Appendix A.4.1 and Global stability analysis in Appendix A.5, which offer a detailed exploration of the stability conditions.

5. Numerical Simulation and Discussion of Results

This section outlines the baseline values used for numerical simulations and discusses the results. The baseline parameters are derived from empirical data sourced from the NNPC and other credible organizations to ensure accuracy and reliability.

5.1. Parametrization and Base Value

The initial condition for the state variables used for simulation are:

$Q\left(0\right)=150Mltr,A\left(0\right)=50-100Mltr,Z\left(0\right)=45-80Mltr,E_1\left(0\right)=20-26Mltr,E_2\left(0\right)=23-36Mltr,M\left(0\right)=45-56Mltr$

for time $t=0$. These estimates are based on data from the 2021 NNPC-FAAC reports [39]–61], NNPC-ASB reports [69], and the NMDPRA [70].

The baseline parameter values are computed as the mean of data collected over at least five consecutive months from the NNPC-FAAC reports. These values in Table 4 provide the foundation for the numerical simulations, enabling an accurate representation of the real-world dynamics in the downstream petroleum sector.

5.1.1. Assumption of Parameters Value

• Values for $\alpha ,k_1,\mathrm{and}{k}_2$ are assumed due to insufficient direct empirical data.
• Loss rates at $E_1\left(t\right)$ and $E_2\left(t\right)$ are denoted as (${\mu }_1$) and is treated uniformly in their respective subsystems for simplicity.

The value of $V\left(t\right)$ used for simulation is the mean (average) of 5- years (2016-2020) vandalization on pipeline

5.1.2. Simulation Setup

The model Equations (5)–(10) were solved numerically using a standard ODE solver. Results were generated for different scenarios, reflecting variations in vandalism activity, government policy, and transporter dynamics. These simulations help visualize the interplay between stakeholders in the downstream sector and assess how system interventions affect fuel distribution efficiency and scarcity resolution.

The next subsection presents the results of the simulations and their implications for addressing the challenges in Nigeria’s downstream sector.

5.1.3. Brief About the Data

The data presented in Figure 5 is from the NNPC FAAC reports highlight the erratic operations of Nigeria’s refineries over the past decade. In 2016, the highest crude oil input to Nigerian refineries was recorded at 20%, but this figure steadily declined to 0% by 2020. This downward trend underscores the diminishing local refining capacity, culminating in a complete halt in refinery operations as of March 2024.

Figure 6 illustrates the Petroleum Products Distribution (PPD) through the NNPC retail network, further confirming that the majority of petroleum products are imported rather than locally refined. The high cost of importation, as noted in [2], remains a significant contributor to the persistent fuel scarcity in Nigeria.

Given that Premium Motor Spirit (PMS) accounts for the highest demand among petroleum products, its distribution data from the NNPC FAAC reports is used as the basis for this study’s simulations. This focus ensures the analysis captures the most critical aspect of Nigeria’s downstream challenges.

5.2. Scenario Simulation and Results Discussion

This section presents numerical simulations for the scenarios described in Section 4 to validate the analytical results. Using baseline values from Table 4, we simulate the minimum threshold (A16) required for petroleum product distribution and consumption. The analysis identifies 42 million liters per day as the minimum supply needed to prevent scarcity before subsidy removal, approximately 3 million liters less than the post-subsidy demand reported in [71].

5.2.1. Effect of Vandalism and Alternative Routes on ($R_e$)

The stability of the minimum supply threshold $R_e$ (Equation (A16))

The minimum supply required to eliminate fuel scarcity, denoted as $R_e$ (Equation (A16)), (Equation (A16)) was analytically proven in Section 5. Simulations assess how pipeline vandalism and alternative transportation routes affect this threshold.

Figure 7 presents the results. Panel (a) shows a contour plot, while panel (b) provides a 3D perspective. The analysis considers $b_1$ values between 0 and 1 and scenarios of high vandalism ($V\approx 1$). The results indicate that even under severe vandalism (blue regions), a 10% provision of effective alternative routes can prevent fuel scarcity. When alternative routes are 100% effective (yellow regions), $R_e$ consistently exceeds demand, ensuring fuel supply stability and minimizing economic disruptions.

5.2.2. Effect of Hoarding and Diversion of Product on the $R_e$

Simulations also examine the impact of hoarding and product diversion under varying levels of vandalism. Figure 8 shows that when vandalism is low (10%), a 10% contribution from NNPC retail hubs and major marketers is sufficient to maintain stable fuel distribution and increase PMS reserves.

Under extreme vandalism (unity), stability is maintained as long as product diversion is minimized. Meeting the 42 million liters per day threshold effectively prevents fuel crises, underscoring the need for robust distribution strategies and efficient supply chain management.

5.2.3. General Model Simulation

The simulation of the entire model provided valuable insight, as shown in Figure 9, Figure 10, Figure 11, Figure 12 and Figure 13.

When vandalism is extremely high (90%) and there are no alternative routes, operators at $Q\left(t\right)$ stop pipeline transport, causing fuel accumulation in the Q compartment (dashed blue plot, Figure 9). This bottleneck disrupts downstream compartments A, Z, $E_1$, $E_2$, and M (Figure 9, Figure 10 and Figure 11), resulting in a critical fuel scarcity crisis.

In contrast, moderate operations ("smooth operation") or alternative transportation routes ($b_1\ne 0$) significantly reduce disruptions. Here, the Q compartment becomes less critical, highlighting the urgent need for alternative transport infrastructure, such as railways, to complement Nigeria’s aging pipeline network.

Introducing alternative routes improves distribution efficiency, reducing reliance on road transportation and minimizing risks such as petroleum truck accidents. These findings align with existing literature [72].

Hoarding by Retail Marketers

Hoarding, primarily by privately owned independent marketers ($E_2$ class), poses a greater threat to consumers than extreme vandalism.

• Figure 11 shows that hoarded fuel accumulates in the $E_2$ compartment (orange plot) instead of reaching consumers (M, Figure 11).
• Without stricter regulations, hoarding disrupts distribution and fuels smuggling to neighboring countries where PMS prices are higher will on the daily basis increase.
• This highlights the urgent need for regulatory measures targeting independent marketers to ensure equitable fuel distribution.

Although the Department of Petroleum Resources (DPR) oversees licensing, safety, and quality control in Nigeria’s oil sector [73], stronger enforcement is needed to curb hoarding and improve supply chain transparency.

 Transportation Efficiency and Industrial Actions 

The model estimates 99% transportation efficiency under low to moderate vandalism (0-90%), even with minimal reliance on alternative routes. However, industrial actions (strikes) by transporters (Z class) introduce severe supply disruptions, as shown in Figure 12 and Figure 13.

• During strikes, the red plots in Figure 12 and Figure 13 indicate severe fuel shortages, demonstrating the system’s vulnerability to industrial actions.
• These findings emphasize the need for policymakers to address stakeholder grievances in the Z class to maintain supply chain stability.

Ensuring a cooperative and functional Z compartment is crucial for sustaining fuel distribution and preventing recurrent crises.

6. Summary, Conclusions, and Limitation

This section synthesizes the findings and implications of the mathematical model developed to address challenges in Nigeria’s downstream petroleum sector. The study emphasizes critical aspects of fuel distribution, including factors such as vandalism, industrial actions, and hoarding, which exacerbate scarcity. It concludes with actionable recommendations for transportation optimization, stakeholder policies, and the potential adaptability of the model to broader contexts. Limitations are acknowledged to inform future research avenues, offering insights into enhancing the robustness of the model and its applications.

6.1. Summary from mathematical Analysis

This study developed a mathematical model to analyze the downstream petroleum sector in Nigeria, specifically focusing on the distribution of petroleum products for domestic consumption. The model incorporated critical factors such as vandalism, industrial actions (strikes), diversion, and hoarding of products by the NNPC Retail hub, independent marketers, and private marketers, all of which contribute to fuel scarcity.

Key Findings:

• Model Validation: The feasibility region, positivity of solutions, and existence and uniqueness of the model were established, ensuring alignment with real-world scenarios and validating the model’s applicability.
• Equilibrium Points: Five equilibrium points were identified, representing conditions under which fuel scarcity could be mitigated or resolved.
• Threshold Analysis: A critical supply PMS threshold of 42 million liters per day was derived, indicating the minimum volume required to meet consumer demand and eliminate queues at service stations.
• Stability Analysis: Using the Lyapunov function and Routh-Hurwitz criteria, the study demonstrated that effective fuel distribution could achieve stability, provided government policies create a conducive environment for stakeholders such as NUPENG, MOMAN, and tanker drivers.
• Simulation Results: Numerical simulations, conducted using the classical Runge-Kutta method in Maple 18, highlighted that:
  • Hoarding by private retail stations exacerbates fuel scarcity and drives up prices of fuel and other commodities.
  • Industrial actions by marketers and transporters (Z) have a more significant impact on fuel distribution crises than pipeline vandalism, underscoring the need for robust government intervention.

6.2. Insights, Challenges, and Future Directions in Fuel Distribution Modeling

This study employs a deterministic, autonomous, and continuous ODE model to analyze fuel distribution in Nigeria’s downstream petroleum sector, focusing on vandalism, industrial action, hoarding, and fuel diversion. Future research could explore hybrid models combining deterministic and stochastic approaches for a more comprehensive analysis.

The model’s insights are adaptable to other countries facing similar fuel distribution challenges, provided robust data is available. However, data limitations constrained this study. Future work should incorporate sensitivity analysis, optimal control, and cost-effectiveness evaluations.

This study considers both pre- and post-subsidy removal scenarios, though simulations were primarily based on the subsidy era to establish demand and supply thresholds for PMS. Future research will extend this work to capture post-subsidy fuel distribution trends more comprehensively.

6.3. Conclusion, Recommendation, and Limitation

This study concludes that alternative transportation methods (e.g., rail), strong policies, and political will are essential for addressing challenges in Nigeria’s downstream sector. Key interventions include:

• Advanced security systems (e.g., smart surveillance) to protect pipelines.
• Infrastructure improvements to mitigate vandalism, corruption, and poor road conditions.
• Effective regulatory oversight of private marketers (Z) to ensure supply stability.

For oil-producing developing nations, reliance on a single transport route poses risks during crises. Implementing alternative routes, as proposed in this study, enhances resilience. Conversely, non-oil-producing nations face significant risks from road-based tanker transportation. Rail infrastructure offers a safer alternative, and governments can apply this model to inform energy security policies.

With the Dangote Refinery now operational, the Nigerian government should explore formal partnerships with neighboring countries to establish legitimate fuel supply channels. ECOWAS could also facilitate regional fuel refining and transportation, reducing smuggling and associated security risks.

Finally, this research highlights critical gaps in fuel consumption data, as NNPCL has yet to transparently disclose actual daily PMS usage. Future studies should focus on data acquisition and validation to enhance policy formulation and sectoral efficiency.

Appendix A.1. Existence and Uniqueness of Solutions

To ensure that the model is mathematically well-posed, we check the Lipschitz condition. The system of Equations (5)–(10) can be written as:

$$\begin{array}{cc}\hfill & F_1=\mathrm{\Lambda }-(a(1-V)+b_1+\mu )Q\hfill \end{array}$$

(A1)

$$\begin{array}{cc}\hfill & F_2=a(1-V)Q+b_{2}Z-(\mu +cZ)A\hfill \end{array}$$

(A2)

$$\begin{array}{cc}\hfill & F_3=b_{1}Q+cAZ-(b_2+\mu +x+y+\eta )Z\hfill \end{array}$$

(A3)

$$\begin{array}{cc}\hfill & F_4=xZ-({\mu }_1+k_1)E_1\hfill \end{array}$$

(A4)

$$\begin{array}{cc}\hfill & F_5=yZ-({\mu }_1+k_2)E_2\hfill \end{array}$$

(A5)

$$\begin{array}{cc}\hfill & F_6=-{\mu }_{2}M+k_1{E}_1+k_2{E}_2+\eta Z\hfill \end{array}$$

(A6)

For Lipschitz continuity, we compute the partial derivatives of (A1), with respect to each state variable $X=\{Q,A,Z,E_1,E_2,M\}$ i.e.

$$\left|{\displaystyle \frac{\partial F_1}{\partial Q}}\right|=|\alpha (1-V)+b_1+\mu |,\left|{\displaystyle \frac{\partial F_1}{\partial A}}\right|=\left|{\displaystyle \frac{\partial F_1}{\partial Z}}\right|=\left|{\displaystyle \frac{\partial F_1}{\partial E_1}}\right|=\left|{\displaystyle \frac{\partial F_1}{{\partial }_2}}\right|=\left|{\displaystyle \frac{\partial F_1}{\partial M}}\right|=0$$

$$\left|{\displaystyle \frac{\partial F_3}{\partial Q}}\right|=b_1,\left|{\displaystyle \frac{\partial F_3}{\partial A}}\right|=cZ,\left|{\displaystyle \frac{\partial F_3}{\partial Z}}\right|=|cA-(b_2+\mu +x+y+\eta )|,\left|{\displaystyle \frac{\partial F_3}{\partial E_1}}\right|=\left|{\displaystyle \frac{\partial F_3}{{\partial }_2}}\right|=\left|{\displaystyle \frac{\partial F_3}{\partial M}}\right|=0$$

Similar steps are applied to the remaining state variables, ensuring that all partial derivatives are bounded. This shows that the model is bounded and satisfies the Lipschitz condition.

We investigate each of the model state variables solutions at any two points in the region $\mathrm{\Omega }\in {\mathbb{R}}^6$ as follows

$$\begin{array}{cc}\hfill |F_1(t,X^2)-F_1(t,X^1)|& =|(\mathrm{\Lambda }-(a(1-V)+b_1+\mu )Q^2)-(\mathrm{\Lambda }-(a(1-V)+b_1+\mu )Q^1)|\hfill \\ \hfill & =|-(a(1-V)+b_1+\mu )(Q^2-Q^1)|\hfill \\ \hfill & \le |(a(1-V)+b_1+\mu )\left|\right|(Q^2-Q^1)|\hfill \end{array}$$

(A7)

Similarly for $F_2$,

$$\begin{array}{cc}\hfill |F_2(t,X^2)-F_2(t,X^1)|& =|(a(1-V)Q+b_{2}Z-(\mu +cZ)A^2)-(a(1-V)Q+b_{2}Z-(\mu +cZ)A^1)|\hfill \\ \hfill & =|-(\mu +cZ)(A^2-A^1)|\hfill \\ \hfill & \le \left|(\mu +cZ)\right||(A^2-A^1)|\hfill \end{array}$$

(A8)

Appendix A.2. Equilibrium Points

To establish Theorem 3, the right-hand sides (RHS) of Equations (5)–(10) are set to zero. On the simplification, we have:

Proof. 

Scenario I: Severe Vandalism

• Critical Collapse: Vandalism is at its maximum ($V\approx 1$), with no alternative transportation route ($b_1=0$) and $Z=0$. Simplifying system (1)–(10) when equated to zero gives the equilibrium point:

  $$E_0^i=\left\{\left(Q,A,Z,E_1,E_2,M\right):\left({\displaystyle \frac{\mathrm{\Lambda }}{\mu }},0,0,0,0,0\right)\right\}.$$

  (A9)
• Partial Collapse: Vandalism remains at maximum ($V\approx 1$), $b_1=0$, but $Z\left(t\right)\ne 0$. Similarly, simplifying system (1)–(10) when equated to zero gives yields:

  $$E_0^{ii}=\left\{{\displaystyle \frac{\mathrm{\Lambda }}{\mu }},{\displaystyle \frac{\mu +x+y+\eta }{c}},{\displaystyle \frac{-\mu }{c}},{\displaystyle \frac{-x\mu }{c({\mu }_1+k_1)}},{\displaystyle \frac{-y\mu }{c({\mu }_1+k_2)}},-{\displaystyle \frac{1}{{\mu }_2}}\left({\displaystyle \frac{\mu +x+y+\eta }{c({\mu }_1+k_1)}}+{\displaystyle \frac{y\mu }{c({\mu }_1+k_2)}}+{\displaystyle \frac{\mu \eta }{c}}\right)\right\}.$$

  (A10)

Scenario II: Reduced Vandalism

• Partial functionality: When vandalism is significantly reduced ($V\ll 1$) and $Z=0$, the system stabilizes at:

  $$E_0^{iii}=\left\{{\displaystyle \frac{\mathrm{\Lambda }}{a(1-V)+\mu }},{\displaystyle \frac{a(1-V)\mathrm{\Lambda }}{\mu \left(a\right(1-V)+\mu )}},0,0,0,0\right\}.$$

  (A11)
• Moderate functionality: Vandalism is reduced ($V\ll 1$), $b_1=0$, and $Z\ne 0$. The equilibrium is:

  $$E_0^{iv}=\left\{{\displaystyle \frac{\mathrm{\Lambda }}{a(1-V)+\mu }},A^{ii},Z^{\ast },{\displaystyle \frac{x{Z}^{\ast }}{{\mu }_1+k_1}},{\displaystyle \frac{y{Z}^{\ast }}{{\mu }_1+k_2}},M^{\ast }\right\},$$

  (A12)

  where:

  $$Z^{\ast }={\displaystyle \frac{\alpha (1-V)(\mathrm{\Lambda }{R}_0-\mu )-{\mu }^2}{c\alpha (1-V)+\mu }},M^{\ast }={\displaystyle \frac{Z^{\ast }}{{\mu }_2}}\left({\displaystyle \frac{x{k}_1}{{\mu }_1+k_1}}+{\displaystyle \frac{y{k}_2}{{\mu }_1+k_2}}+{\displaystyle \frac{\eta }{1}}\right).$$

Scenario III (Optimal Distribution): Increased Fuel Demand and Systemic Challenges

For this scenario, $b_1\ne 0,b_2\ne 0$, $V\ll 1$, and $Z\ne 0$. These are considered while solving (1)–(10) and the equilibrium is:

$$\begin{array}{cc}\hfill E_0^v=& \{{\displaystyle \frac{\mathrm{\Lambda }}{a(1-V)+b_1+\mu }},{\displaystyle \frac{a(1-V)\mathrm{\Lambda }+b_2{Z}^{\ast \ast }}{(\mu +c{Z}^{\ast \ast })(a(1-V)+b_1+\mu )}},Z^{\ast \ast },\hfill \\ \hfill & {\displaystyle \frac{x{Z}^{\ast \ast }}{{\mu }_1+k_1}},{\displaystyle \frac{y{Z}^{\ast \ast }}{{\mu }_2+k_2}},M^{iv}\},\hfill \end{array}$$

(A13)

where:

$$M^{iv}={\displaystyle \frac{Z^{\ast \ast }\left(k_{1}x({\mu }_2+k_2)+({\mu }_1+k_1)\left(k_{2}y+\eta ({\mu }_1+k_2)\right)\right)}{{\mu }_2({\mu }_1+k_1)({\mu }_1+k_2)}}.$$

□

Appendix A.3. Proof of the Threshold Parameter

Proof. 

Two scenarios are considered: product distribution under high vandalism (with and without alternative routes).

$$\begin{array}{cc}\hfill & F=\left(\begin{array}{ccc}cA& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right)\mathrm{and}\hfill \\ & V=\left(\begin{array}{ccc}{b}_2+\mu +x+y+\eta & 0& 0\\ -x& {\mu }_1+k_1& 0\\ -y& 0& {\mu }_1+k_2\end{array}\right)\hfill \end{array}$$

(A14)

Matrix F represents the volume of petroleum transported via pipelines to storage facilities, while matrix V accounts for transitions leading to product delivery to consumers.

If alternative routes are unavailable $b_1=b_2=0$, the system relies solely on pipelines. The inverse of V is given by:

$$V^{-1}=\left(\begin{array}{ccc}{\left(\mu +x+y+\eta \right)}^{-1}& 0& 0\\ {\displaystyle \frac{x}{\left(\mu +x+y+\eta \right)\left({\mu }_1+k_1\right)}}& {\left({\mu }_1+k_1\right)}^{-1}& 0\\ {\displaystyle \frac{y}{\left(b_2+\mu +x+y+\eta \right)\left({\mu }_1+k_2\right)}}& 0& {\left({\mu }_1+k_2\right)}^{-1}\end{array}\right)$$

The threshold parameter, which accounts for both smooth distribution and disruptions, is determined by the spectral radius of $F{V}^{-1}$. The results for the two scenarios are:

$$R_e^a={\displaystyle \frac{cA+b_{1}Q}{(b_2+\mu +x+y+\eta )}}\mathrm{and}{R}_e={\displaystyle \frac{cA}{(\mu +x+y+\eta )}}$$

(A15)

At the Petroleum Distribution-Free Equilibrium (PDE), substituting the value of A in Equation (A15) yields:

$$R_e={\displaystyle \frac{c\mathrm{\Lambda }\alpha (1-V)}{\mu k\left(\alpha (1-V)+\mu \right)}}$$

(A16)

□

Appendix A.4. Stability Analysis

Appendix A.4.1. Application of the Methodology

The stability theorem (5) is proved using the methodology outlined in Section 3. The study focuses on improving fuel distribution to mitigate scarcity, leveraging concepts from mathematical epidemiology. Accordingly, stability is analyzed around the equilibrium points corresponding to Scenarios II and III, encompassing partial functionality, moderate functionality, and optimal distribution states.

Two specific equilibrium points are analyzed:

• Products Distribution-Free Equilibrium (PDE): Equation (A12), representing a scenario where products are distributed without an alternative route.
• Products Distribution Persistent Equilibrium (PDP): Equation (A13), representing a scenario where products are distributed with an alternative route.

Proof of Local Stability

To demonstrate local asymptotic stability, the following steps are undertaken: The Jacobian matrix J of the system is defined as:

$$J=\left(\begin{array}{cccccc}-(\alpha (1-V)+b_1+\mu )& 0& 0& 0& 0& 0\\ \alpha (1-V)& -(\mu +cZ)& b_2-cA& 0& 0& 0\\ b_1& cZ& cA-(b_2+\mu +x+y+\eta )& 0& 0& 0\\ 0& 0& x& -({\mu }_1+k_1)& 0& 0\\ 0& 0& y& 0& -({\mu }_1+k_2)& 0\\ 0& 0& \eta & k_1& k_2& -{\mu }_2\end{array}\right)$$

(A17)

• Step 1: Linearization: The Jacobian matrix of the system is obtained and evaluated near the partial functionality equilibrium and for a specific case; in the absence of alternative routes, the Jacobian reduces to:

  $$J\left(E_0^{\left(iii\right)}\right)=\left[\begin{array}{cccccc}-\alpha (1-V)-\mu & 0& 0& 0& 0& 0\\ \alpha (1-V)& -\mu & -cA& 0& 0& 0\\ 0& 0& cA-\eta -\mu -x-y& 0& 0& 0\\ 0& 0& x& -{\mu }_1-k_1& 0& 0\\ 0& 0& y& 0& -{\mu }_1-k_2& 0\\ 0& 0& \eta & k_1& k_2& -\mu \end{array}\right]$$

  (A18)
• Step 2: Eigenvalue Computation: The eigenvalues are determined from the characteristic equation $\left|J\right(E_0^{\left(iii\right)})-\lambda I|=0$.
• Step 3: Characterization of Eigenvalues: The eigenvalues of the Jacobian are:

  $${\lambda }_1=-\mu ,{\lambda }_2=-{\mu }_1-k_1,{\lambda }_3=-{\mu }_1-k_2,{\lambda }_4=-\mu ,{\lambda }_5=-\alpha (1-V)-\mu ,{\lambda }_6=cA-(\mu +x+y+\eta ).$$
• Step 4: Stability Criteria: By the Routh-Hurwitz stability criterion, all eigenvalues must be negative for local asymptotic stability.

  -

  ${\lambda }_1,{\lambda }_2,{\lambda }_3,{\lambda }_4<0$ hold true.

  -

  For ${\lambda }_5<0$, the vandalism rate must satisfy $1-V>0$.

  -

  For ${\lambda }_6<0$, the consumer demand must exceed the supply at retailer filling stations, implying $cA-(\mu +x+y+\eta )>0$.

  Consequently, stability depends on controlling vandalism and ensuring the supply chain meets consumer demand. The model is locally asymptotically stable under these conditions, affirming Theorem 5.

Reduced System Analysis

To simplify the stability analysis, fuel distribution is narrowed to major petroleum dealers and suppliers in the community. This enables qualitative analysis of the subsystem:

$$\begin{array}{cc}\hfill & {\displaystyle \frac{dQ}{dt}}=\mathrm{\Lambda }-(a(1-V)+\mu )Q,\hfill \\ \hfill & {\displaystyle \frac{dA}{dt}}=a(1-V)Q-(\mu +cZ)A,\hfill \\ \hfill & {\displaystyle \frac{dZ}{dt}}=cAZ-kZ,\hfill \end{array}$$

(A19)

where $k=\mu +x+y+\eta$.

The equilibrium points for this subsystem are:

• PDE ($P_0$): $\left({\displaystyle \frac{\mathrm{\Lambda }}{g_1+\mu }},{\displaystyle \frac{\mathrm{\Lambda }{g}_1}{\mu (g_1+\mu )}},0\right)$,
• PDP ($P_1$): $\left({\displaystyle \frac{\mathrm{\Lambda }}{g_1+\mu }},{\displaystyle \frac{k}{c}},{\displaystyle \frac{\mu (R_0-1)}{c}}\right)$,

where $g_1=a(1-V)$ and $R_0={\displaystyle \frac{\mathrm{\Lambda }c{g}_1}{\mu k(g_1+\mu )}}$.

Global Stability Analysis

$$J(Q,A,Z)=\left(\begin{array}{ccc}-(g_1+\mu )& 0& 0\\ g_1& -(\mu +cZ)& -cA\\ 0& cZ& cA-k\end{array}\right)$$

(A20)

• At $P_0$: The Jacobian matrix evaluated at $P_0$ yields the characteristic equation:

  $$(\lambda +g_1+\mu )(\lambda +\mu )\left(\lambda -{\displaystyle \frac{c\mathrm{\Lambda }{g}_1}{\mu (g_1+\mu )}}+k\right)=0.$$

  (A21)

  Eigenvalues are:

  $${\lambda }_1=-(g_1+\mu ),{\lambda }_2=-\mu ,{\lambda }_3=k(R_0-1).$$

  Stability holds if $R_0<1$, establishing $P_0$ as globally asymptotically stable under this condition.
  Theorem A1. 

  The stability at $P_0$ is Globally Asymptotically stable if $R_0<1$ and the system transits to instability $R_0>1$.
  Proof:
  Consider the Lyapnouv function $F=Z$ with derivative of

  $$\begin{array}{cc}\hfill \dot{F}& =\dot{Z}\hfill \end{array}$$

  (A22)

  $$\begin{array}{cc}\hfill & =cAZ-kZ=Z(cA-k)\hfill \end{array}$$

  (A22)

  from PDE (partial functionality), $A\le {\displaystyle \frac{\mathrm{\Lambda }{g}_1}{\mu (g_1+\mu )}}$, (A23) becomes

  $$\dot{F}\le Z\left({\displaystyle \frac{c\mathrm{\Lambda }{g}_1}{\mu (g_1+\mu )}}-k\right)$$

  (A24)

  simplification gives

  $$\begin{array}{cc}\hfill \dot{F}& \le kZ\left({\displaystyle \frac{c\mathrm{\Lambda }{g}_1}{k\mu (g_1+\mu )}}-1\right)\hfill \\ \hfill & =kz(R_0-1)\hfill \end{array}$$

  (A25)

  Therefore $\dot{F}\le 0$ if $R_0\le 1$ and $\dot{F}=0$ if $z=0$. Thus, F is a Lyapnouv function and by LaSellei invariance principle, every solution of (A19) approaches $P_0$ as $t\to \infty$.
• At $P_1$: Using the Goh-Volterra nonlinear Lyapunov function, global stability at $P_1$ is confirmed if $R_0>1$.

Appendix A.5. Global Stability at Moderate Functionality

Proof. 

To establish the system’s global stability, it is a suffix to use the reduced system (A19) at $P_1$ by formulating the Goh-Volterra nonlinear Lyapnouv function

$$\begin{array}{cc}\hfill L& =\left(Q-Q^{\ast \ast }-Q^{\ast \ast }ln{\displaystyle \frac{Q}{Q^{\ast \ast }}}+\right)+\left(A-A^{\ast \ast }-A^{\ast \ast }ln{\displaystyle \frac{A}{A^{\ast \ast }}}\right)\hfill \\ & +\left(Z-Z^{\ast \ast }-Z^{\ast \ast }ln{\displaystyle \frac{Z}{Z^{\ast \ast }}}\right)\hfill \end{array}$$

(A26)

Differentiate the Lyapnouv function (A26) gives:

$$\begin{array}{cc}\hfill \dot{L}& =\left(1-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)\dot{Q}+\left(1-{\displaystyle \frac{A^{\ast \ast }}{A}}\right)\dot{A}+\left(1-{\displaystyle \frac{Z^{\ast \ast }}{Z}}\right)\dot{Z}\hfill \end{array}$$

(A27)

substituting (A19) in (A27) to get

$$\begin{array}{cc}\hfill \dot{L}& =\left(1-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)(\mathrm{\Lambda }-(g_1+\mu )Q)\hfill \\ & +\left(1-{\displaystyle \frac{A^{\ast \ast }}{A}}\right)(g_{1}Q-(cZ+\mu )A)+\left(1-{\displaystyle \frac{Z^{\ast \ast }}{Z}}\right)(cZA-kZ)\hfill \end{array}$$

(A28)

substitute $\mathrm{\Lambda }=(g_1+\mu )Q^{\ast \ast }$ into (A28) to obtain

$$\begin{array}{cc}\hfill \dot{L}& =\left(1-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)\left((g_1+\mu )Q^{\ast \ast }-(g_1+\mu )Q\right)\hfill \\ & +\left(1-{\displaystyle \frac{A^{\ast \ast }}{A}}\right)(g_{1}Q-(cZ+\mu )A)+\left(1-{\displaystyle \frac{Z^{\ast \ast }}{Z}}\right)(cZA-kZ)\hfill \end{array}$$

(A29)

Expanding and simplification of (A29) gives

$$\begin{array}{cc}\hfill \dot{L}& =(g_1+\mu )Q^{\ast \ast }\left(2-{\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)+g_{1}Q-\mu A-g_{1}Q{\displaystyle \frac{A^{\ast \ast }}{A}}\hfill \\ & +(cZ+\mu )A^{\ast \ast }-kZ-A{Z}^{\ast \ast }+k{Z}^{\ast \ast }\hfill \end{array}$$

(A30)

Let $k{Z}^{\ast \ast }=c{A}^{\ast \ast }{Z}^{\ast \ast }$, we have

$$\begin{array}{cc}\hfill \dot{L}& =(g_1+\mu )Q^{\ast \ast }\left(2-{\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)+g_1{Q}^{\ast \ast }\left({\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q{A}^{\ast \ast }}{Q^{\ast \ast }A}}\right)\hfill \\ & +(cZ+\mu )A^{\ast \ast }-cZ{A}^{\ast \ast }-A{Z}^{\ast \ast }+c{A}^{\ast \ast }{Z}^{\ast \ast }-\mu A\hfill \end{array}$$

(A31)

factorizing the common terms and further simplification gives

$$\begin{array}{cc}\hfill \dot{L}& =(g_1+\mu )Q^{\ast \ast }\left(2-{\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)+g_1{Q}^{\ast \ast }\left({\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q{A}^{\ast \ast }}{Q^{\ast \ast }A}}\right)\hfill \\ & +(\mu +c{Z}^{\ast \ast })A^{\ast \ast }\left(1-{\displaystyle \frac{A}{A^{\ast \ast }}}\right)\hfill \end{array}$$

(A32)

we observed from PDE (partial functionality) that $g_1{Q}^{\ast \ast }=(\mu +c{Z}^{\ast \ast })A^{\ast \ast }$, thus, we substitute in (A32) to get

$$\begin{array}{cc}\hfill \dot{L}& =\mu Q^{\ast \ast }\left(2-{\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\right)+g_1{Q}^{\ast \ast }\left(3-{\displaystyle \frac{Q^{\ast \ast }}{Q}}-{\displaystyle \frac{A}{A^{\ast \ast }}}-{\displaystyle \frac{Q{A}^{\ast \ast }}{Q^{\ast \ast }A}}\right)\hfill \end{array}$$

(A33)

Since the arithmetic mean exceeds the geometric mean, it follows that

$$2-{\displaystyle \frac{Q}{Q^{\ast \ast }}}-{\displaystyle \frac{Q^{\ast \ast }}{Q}}\le 0$$

and

$$3-{\displaystyle \frac{Q^{\ast \ast }}{Q}}-{\displaystyle \frac{A}{A^{\ast \ast }}}-{\displaystyle \frac{Q{A}^{\ast \ast }}{Q^{\ast \ast }A}}\le 0$$

In this view, the relationship between the arithmetic mean of being greater than or equal to the geometric mean supports our conclusion about $\dot{L}\le 0$ for $R_0>1$ and $\dot{L}=0$ for $R_0=1$.

This completed the last part of the Theorem 5 □