A Hybrid Fuzzy Mathematical Programming Approach for Manufacturing Inventory Models with Partial Trade Credit Policy and Reliability

Abstract

This study introduces an inventory model for manufacturing that prioritizes product quality and cost efficiency. Utilizing fuzzy logic and mathematical programming, the model integrates fuzzy numbers to describe uncertainties associated with manufacturing costs and quality control parameters. The model extends beyond conventional inventory systems by incorporating a dynamic mechanism to halt production, employing fuzzy decision variables to optimize the economic order quantity and minimize total costs. Key innovations include the application of approaches related to graded mean integration for defuzzification and the use of Kuhn–Tucker conditions to ensure optimal solutions under complex constraints. These approaches facilitate the precise management of production rates, inventory levels, and cost factors, which are essential in achieving a balance between supply and demand. A computational analysis validates the model’s effectiveness, demonstrating cost reductions while maintaining optimal inventory levels. This underscores the potential of integrating fuzzy arithmetic with traditional optimization techniques to enhance decision making in inventory management. The model’s adaptability and accuracy indicate its broad applicability across various sectors facing similar challenges, offering a valuable tool for operational managers and decision makers to improve efficiency and reduce waste in production cycles.

1. Introduction

Trade credit policies have emerged as pivotal strategies for manufacturers aiming to stimulate retailer demand and boost product sales, especially in current highly competitive markets. These strategies afford retailers a specified period to settle payments for acquired goods, typically referred to as a single-level credit policy or permissible payment delay [1]. However, manufacturing processes are often imperfect, resulting in the generation of defective items alongside flawless products. Identified during quality inspections, these defective items are generally repaired, thereby recovering potential revenue [2,3,4,5]. A wealth of research has considered inventory models with a focus on contexts related to the reliability of the production process and maintenance policies [6,7,8,9,10,11,12,13,14]. Techniques from biological and mathematical modeling have been successfully applied to these contexts, providing a diverse range of solutions [15,16,17,18,19,20,21,22].

Research on lot-size models with uncertain demand, accounting for skewness and kurtosis and employing stochastic programming, has been applied to complex scenarios such as hospital pharmacies using COVID-19 data [23]. Additionally, inventory management for new products with triangularly distributed demand and lead time has also been studied [24]. Further research has explored the optimization of inventory costs and supply chain success under various complex conditions [25,26,27]. Recent studies have also addressed nonlinear quantity discounts and price-dependent demand, adding valuable insights for the management of inventories [28,29]. In certain circumstances, retailers extend credit facilities to their customers, resulting in a two-level trade credit policy. In such cases, fuzzy set theory provides a robust framework for managing imprecise or vague data by assigning degrees of membership to each object [30,31,32,33,34]. It distinguishes itself from traditional methods through its ability to handle uncertainty in a more nuanced manner. Furthermore, the graded mean integration (GMI) method plays a key role in reducing the uncertainty associated with ambiguous data, making it especially helpful for inventory models and optimization processes in uncertain environments [35,36,37,38,39,40,41]. Recent advancements have demonstrated its utility in diverse applications, including the management of pandemic models [42]. Moreover, contributions have been made in methodologies for the consolidation of the effects of inventory management with serially dependent random demand, offering robust solutions for the optimization of inventory strategies in environments with fluctuating demand [43]. A case study on drug supply further demonstrated the applicability of stochastic programming to effectively manage time-dependent demand [44].

In the current fiercely competitive manufacturing landscape, ensuring product reliability is paramount. To enhance accuracy in handling uncertainties associated with costs and quality controls, fuzzy set theory is an essential tool, offering an improvement over single-valued numbers (traditional or crisp set theory) that may provide exact but less realistic values. Fuzzy methodologies such as those applied in reliability improvement processes [45] and product design systems [46] have proven effective in managing uncertainties in complex industrial systems. For defuzzification, the GMI method can be implemented, allowing for more precise interpretation of fuzzy data. Moreover, hybrid multi-dimensional fuzzy approaches have shown potential to improve reliability and safety in various sectors [47], further proving the applicability of fuzzy logic in enhancing decision-making processes. A notable gap in the literature lies in the development of a manufacturing inventory model that effectively manages deteriorating goods while incorporating price-dependent demands and a production rate that is sensitive to reliability under a partial trade credit policy. The present research aims to address this gap by exploring the complex relationship between reliability and production total cost (TC), providing valuable insights for strategic decision making in inventory management. To achieve this, we employ logistic regression (LR), a powerful supervised machine learning (ML) technique specifically used here to predict the optimal cost conditions of the system. This technique not only supports decision making but also demonstrates the advantage of integrating fuzzy arithmetic with ML to improve the precision and adaptability of inventory management.

Our main objective is to develop a manufacturing inventory model that manages deteriorating products and incorporates price-dependent demands, as well as a production rate sensitive to reliability under a partial trade credit policy. Specifically, we implement an economic order quantity (EOQ) model, a fundamental tool in inventory management [48], and illustrate how fuzzy arithmetic is utilized within the model to handle the inherent complexities of real-world scenarios. Additionally, we apply nonlinear mathematical programming using the Lagrangian method to solve the EOQ model and evaluate its efficiency compared to other techniques. By leveraging the Weka software (version 3.8.6), we analyze the TCs by classifying them into profitable and non-profitable segments based on cycle length, optimizing inventory costs across varying market conditions.

The remainder of this article is organized as follows. Section 2 describes the methodology used to develop the fuzzy inventory model, including background on fuzzy set theory and its application in the model. In Section 3, we detail the mathematical formulation of the TC and optimization methods employed in our methodology. In Section 4, the findings derived from our research are discussed. Section 5 provides the conclusions drawn from the analysis and suggests ideas for future research.

2. Fuzzy Inventory Model

This section presents the methodology for developing a manufacturing inventory model aimed at optimizing costs and managing inventory under uncertain conditions, utilizing fuzzy numbers and nonlinear programming techniques.

2.1. Fuzzy Inventory Model Components

We recall that the primary objective of this study is to develop a manufacturing inventory model that manages deteriorating products and incorporates price-dependent demands, as well as a production rate sensitive to reliability, under a partial trade credit policy. We apply fuzzy set theory to address uncertainties in costs and quality control, using the GMI method for defuzzification to provide more accurate and reliable interpretations of the data. Before presenting the formulation of the proposed model, we state the notations utilized in the model by means of

Table 1. Notations used in the inventory model.

Symbol	Description	Unit
$\alpha$	Portion of the purchase price covered by the buyer, where $0\le \alpha \le 1$	proportion
d	Holding cost per unit per time	$/unit
e	Requirement scale component	constant
f	Regular production cost	$/unit
G	Manufacturing rate	units/time
H	Production capacity or target production level	units
h	Efficiency factor affecting production halt time	constant
$J_c$	Interest charged by the producer	$/unit year
$J_e$	Manufacturer interest earned	$/unit
$J\left(t\right)$	Inventory rate function at time t, where $0\le t\le T$	units/time
L	Cycle period	years
M	Sale price	$/unit
o	Supplier trade credit period	years
⌀	Exponentially growing deterioration rate	constant
$\vartheta$	Adjustment factor for deterioration or efficiency in credit terms	constant
p	Retailer credit term period for customers	years
$R_0$	Resupply cost per order	$/unit
$\rho$	Supply-limiting factor	constant
${\sigma }_e$	Reliability rate for production of quality goods	percentage
$t_1$	Manufacturer halt period	years
$\mathrm{TC}\left(L\right)$	Total cost of the system as a function of the cycle period L	$/year

.

Next, we outline the following primary components of our optimization framework, that is, the decision variable, objective function, and constraints that govern the model:

• Decision variable (L), which corresponds to the production cycle period;
• Objective function, which must minimize the TC of the system ($\mathrm{TC}\left(L\right)$), including the ordering cost, carrying cost, deterioration cost, interest paid, and interest earned, which is formulated as

  $$min\left\{\mathrm{TC}\left(L\right)\right\}={\displaystyle \frac{R_0}{T}}+{\displaystyle \frac{(d+⌀f+f{J}_c)({\sigma }_{e}G{t}_1-LT)}{⌀T}}+{\displaystyle \frac{f{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-H)}{{⌀}^{2}T}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M{J}_e}{2T}};$$

  (1)
• The constraint is stated as

  $$t_1={\displaystyle \frac{hL}{{\sigma }_{e}G}}+{\displaystyle \frac{1}{2}}\left({\displaystyle \frac{h⌀}{{\sigma }_{e}G}}-{\displaystyle \frac{h^2⌀}{{\sigma }_e^{2}G}}\right)L^2,$$

  (2)

  where $0\le t\le T$ represents the time constraints ensuring that the production and inventory processes are bounded within the cycle period; ${\sigma }_{e}G{t}_1\ge LT$ ensures that the production meets or exceeds the demand over the cycle; $L>0$ is the cycle period, which must be positive; $t_1\ge 0$ is the production halt time, which must be non-negative; and $0\le \alpha \le 1$ ensures the portion of the purchase price covered by the buyer is between 0 and 1.

The constraint presented in (2) defines the production halt time ($t_1$), which is essential for balancing the production rate and demand. From a practical standpoint, $t_1$ represents the time during which the manufacturer halts production to avoid overproduction or accumulation of inventory beyond capacity. In real-world scenarios, manufacturing systems must balance the rates of production and demand fulfillment. For example, in industries where products are perishable or have a limited shelf life, such as pharmaceuticals or food elaboration, producing too much can lead to excessive costs associated with spoilage or storage.

The first term of the constraint presented in (2), $hL/{\sigma }_{e}G$ say, shows that the halt time is proportional to the production cycle length (L) and inversely proportional to the production rate (G). As production reliability (${\sigma }_e$) improves, the system becomes more efficient, reducing halt times. Conversely, lower reliability requires longer halt periods to minimize defects.

The second term of the constraint stated in (2), $(h⌀/{\sigma }_{e}G-h^2⌀/{\sigma }_e^{2}G{L}^2)/2$ namely, introduces nonlinearity into the system, reflecting the increasing marginal cost or effort required to adjust production in response to growing inefficiencies over time. This is particularly relevant in scenarios where the rate of deterioration (⌀) strongly impacts production and storage decisions. Industries that handle fragile goods, such as electronics or high-tech equipment, can leverage this to optimize production cycles and minimize wastage due to quality deterioration over time.

In practical terms, the constraint presented in (2) ensures that the production process is dynamically adjusted according to demand fluctuations, ensuring that the cycle period (L) is aligned with the demand rate and quality control mechanisms. It also ensures that the production halt time ($t_1$) is non-negative, guaranteeing that the system remains feasible, even under varying production conditions.

The values of the constants in the model (see Table 1), such as $\vartheta$, are generally estimated from empirical data, historical studies, or sensitivity analyses tailored to specific industry contexts. For instance, the constant $\vartheta$ reflects the sensitivity of the deterioration rate to credit terms and can be adjusted according to the characteristics of the products (such as perishable goods or electronics) and the observed trade credit policies. In practice, these values may vary depending on the specific operational conditions and can be refined through simulations to match real-world scenarios more closely.

The inventory model relies on various mathematical and fuzzy logic concepts. For fuzzy logic, we use trapezoidal (TFNs), pentagonal (PFNs), and hexagonal (HFNs) fuzzy numbers to represent uncertain variables. These numbers define the boundaries of the fuzzy set, capturing the range of possible values for uncertain variables.

To provide a clearer context, we present a comparison between the proposed model and previous works in

Table 2. Comparison of proposed the fuzzy inventory model with previous works (part 1).

Criteria	Lee and Rosenblatt (1987) [1]	Urban (1992) [2]	Panda et al. (2019) [4]
Objective and methodology	Joint control of production cycles and maintenance by inspection (deterministic approach)	Optimal lot size, price mark-up, and advertising expenditure (separable programming)	Two-warehouse inventory model with price- and stock-dependent demand (solved using Lingo software, version 10.0)
Demand type and deterioration	Constant; no deterioration	Price-dependent deterministic; no deterioration	Price- and stock-dependent; deterioration considered with partial backlogging
Trade credit policy	Not considered	Not considered	Partial trade credit incorporated
Fuzzy logic and reliability	No fuzzy logic or reliability considered	No fuzzy logic or reliability considered	No fuzzy logic or reliability considered
Optimization technique	Approximation to cost function	Separable programming	Nonlinear mathematical modeling (Lingo software)

,

Table 3. Comparison of proposed fuzzy inventory model with previous works (part 2).

Criteria	Sarkar et al. (2014) [7]	Singh and Sharma (2017) [9]	Das et al. (2021) [49]
Objective and methodology	Economic manufacturing quantity with imperfect production and reliability (Euler–Lagrange method)	Decaying items with stochastic demand and inflation effects (stochastic programming)	Production inventory with partial trade credit and reliability (Lagrangian method)
Demand type and deterioration	Time-dependent demand; deterioration considered	Stochastic demand; deterioration considered	Not explicitly specified; reliability integrated into inventory management
Trade credit policy	Not explicitly specified	Not explicitly specified	Partial trade credit incorporated
Fuzzy logic and reliability	Reliability considered; no fuzzy logic	Reliability considered; no fuzzy logic	Reliability and fuzzy logic both considered
Optimization technique	Euler–Lagrange method	Stochastic programming	Lagrangian method (fuzzy set theory)

and

Table 4. Comparison of proposed fuzzy inventory model (part 3).

Criteria	Proposed Model
Objective and methodology	Develop fuzzy inventory model for deteriorating products, price-dependent demand, and partial trade credit with reliability (fuzzy set theory, GMI, LR, and ML)
Demand type and deterioration	Price-dependent with fuzzy uncertainty; deterioration adjusted by reliability factor ($\vartheta$)
Trade credit policy	Partial trade credit policy incorporated
Fuzzy logic and reliability	Fuzzy logic and reliability both considered
Optimization technique	Nonlinear mathematical programming, Lagrangian method, and fuzzy logic integration

. These tables highlight the key contributions and differences, particularly in terms of handling uncertainty, trade credit policies, and inventory management under deteriorating conditions.

2.2. Fuzzy Numbers

Before defining specific fuzzy numbers, it is important to understand the concept of fuzzy sets and membership functions, as well as to establish the difference between fuzzy sets and fuzzy numbers.

Let $\mathcal{X}$ be the universe of discourse, which is the set of all possible elements under consideration. In the context of our fuzzy inventory model, $\mathcal{X}$ could include elements such as the production rate, holding cost, deterioration rate, sales price, and other key factors listed in Table 1. These elements represent the various parameters that are subject to uncertainty.

To establish membership, a fuzzy set is defined as $\mathcal{F}\equiv \{(x,{\mu }_{\mathcal{F}}\left(x\right)):x\in \mathcal{X}\},$ where ${\mu }_{\mathcal{F}}\left(x\right)$ is the function that assigns a degree of membership in the interval of $[0,1]$ to each element $x\in \mathcal{X}$.

The support of a fuzzy set ($\mathcal{F}$) is the subset of elements with a membership greater than zero, stated as $\mathrm{support}\left(\mathcal{F}\right)\equiv \{x\in \mathcal{X}:{\mu }_{\mathcal{F}}\left(x\right)>0\}$. The core of a fuzzy set ($\mathcal{F}$) is the subset of elements with a membership equal to one, stated as $\mathrm{core}\left(\mathcal{F}\right)\equiv \{x\in \mathcal{X}:{\mu }_{\mathcal{F}}\left(x\right)=1\}$.

A fuzzy number is a specific type of fuzzy set that represents a quantity with uncertainty. Unlike general fuzzy sets, fuzzy numbers are often defined on the real number line and have specific shapes, such as hexagonal, pentagonal, and trapezoidal forms. These shapes are characterized by their membership functions, which define how the degree of membership varies over the range of possible values.

In the context of our inventory model, fuzzy numbers can represent various uncertain parameters, such as the holding cost, deterioration rate, and manufacturing rate. By using these specific forms of fuzzy numbers, we may model and manage the uncertainties in these parameters, leading to more robust decision making.

Different shapes of fuzzy numbers offer varying levels of detail and complexity in representing uncertainty as follows:

• TFNs provide a good approximation of uncertainty by allowing a plateau of full membership, which can model situations with a more defined range of most likely values.
• PFNs and HFNs offer even more flexibility, allowing for a more detailed representation of uncertainty, especially in cases where there are multiple levels of probability or asymmetric distributions.

Choosing the appropriate form of fuzzy numbers depends on the specific characteristics of the uncertainty being modeled and the available computational resources. More complex shapes can provide a more accurate representation of uncertainty but at the cost of increased computational complexity. We now proceed to define the $L,R$-type fuzzy numbers. Let $L,R:[0,1]\to [0,1]$ be two fixed functions that are both upper semicontinuous and decreasing such that $L\left(0\right)=R\left(0\right)=1$ and $L\left(1\right)=R\left(1\right)=0$. A fuzzy number (A) is said to be an $L,R$-type fuzzy number if the membership function (${\mu }_A$) is given by

$${\mu }_A\left(x\right)=\left\{\begin{array}{cc}L\left({\displaystyle \frac{C_L-x}{\eta }}\right),\hfill & \mathrm{if}{C}_L-\eta \le x<C_L;\hfill \\ 1,\hfill & \mathrm{if}{C}_L\le x\le C_U;\hfill \\ R\left({\displaystyle \frac{x-C_U}{\gamma }}\right),\hfill & \mathrm{if}{C}_U<x\le C_U+\gamma ;\hfill \\ 0,\hfill & \mathrm{otherwise};\hfill \end{array}\right.$$

where $\eta$ and $\gamma$ are the left and right spreads (widths) of ${\mu }_A$, respectively, and $[C_L,C_U]$ is the core of ${\mu }_A$. An $L,R$-type fuzzy number is denoted as $A\equiv \{(C_L,C_U;\eta ,\gamma \}$.

The membership function ${\mu }_A$ describes how each point (x) in the universe of discourse ($\mathcal{X}$) is mapped to a degree of membership ranging from zero to one. The left spread ($\eta$) determines how the membership value decreases from one to zero as x moves from $C_L$ to $C_L-\eta$, while the right spread ($\gamma$) determines how the membership value decreases from one to zero as x moves from $C_U$ to $C_U+\gamma$.

Next, we introduce the generalized TFNs. A quadruple $T\equiv \{(t_1,t_2,t_3,t_4);u,{\mu }_T$, $t_1\le t_2\le t_3\le t_4$, $t_1,t_2,t_3,t_4\in \mathbb{R}$, $u\in [0,1]\}$ is called a generalized TFN if its membership function is formulated as

$${\mu }_T\left(x\right)=\left\{\begin{array}{cc}{\mu }_{T_L}\left(x\right),\hfill & \mathrm{if}{t}_1\le x<t_2;\hfill \\ u,\hfill & \mathrm{for}{t}_2\le x<t_3;\hfill \\ {\mu }_{T_R}\left(x\right),\hfill & \mathrm{if}{t}_3\le x\le t_4;\hfill \\ 0,\hfill & \mathrm{otherwise};\hfill \end{array}\right.$$

where ${\mu }_{T_L}\left(x\right)=u(\left(x-t_1\right)/\left(t_2-t_1\right))$ and ${\mu }_{T_R}\left(x\right)=u(\left(t_4-x\right)/\left(t_4-t_3\right))$. If $t_2=t_3$, then the generalized TFN (T) becomes a generalized TFN and is denoted by the triplet $T=\{(t_1,t_2,t_3);u,{\mu }_T\}$.

Next, we define generalized PFNs. A fuzzy subset $P\equiv \{(p_1,p_2,p_3,p_4,p_5);$$u,u_L,u_R,{\mu }_P\}$, with $p_1\le p_2\le p_3\le p_4\le p_5$, and $p_i\in \mathbb{R}$, for $i\in \{1,\cdots ,5\}$, and $u,u_L,u_R\in [0,1]$, is called a generalized PFN is its membership function is presented as

$${\mu }_P\left(x\right)=\left\{\begin{array}{cc}{\mu }_L^1\left(x\right),\hfill & \mathrm{if}{p}_1\le x<p_2;\hfill \\ {\mu }_L^2\left(x\right),\hfill & \mathrm{if}{p}_2\le x<p_3;\hfill \\ u,\hfill & \mathrm{if}{p}_3\le x<p_4;\hfill \\ {\mu }_R^2\left(x\right),\hfill & \mathrm{if}{p}_4\le x\le p_5;\hfill \\ 0,\hfill & \mathrm{otherwise};\hfill \end{array}\right.$$

where ${\mu }_L^1\left(x\right)=u_L(\left(x-p_1\right)/\left(p_2-p_1\right))$, ${\mu }_L^2\left(x\right)=u_L+(u-u_L)(\left(x-p_2\right)/\left(p_3-p_2\right))$, ${\mu }_R^2\left(x\right)=u_R+(u-u_R)(\left(p_5-x\right)/\left(p_5-p_4\right)$, and ${\mu }_R^1\left(x\right)=u_R(\left(p_5-x\right)/\left(p_5-p_4\right))$ are the left lower, left upper, right upper, and right lower legs of P, respectively. A generalized PFN passes through $(p_1,0)$, $(p_2,u_L)$, $(p_3,u)$, $(p_4,u_R)$, and $(p_5,0)$. Hence, it is denoted by $P\equiv \{(p_1,p_2,p_3,p_4,p_5);u,u_L,u_R,{\mu }_P\}$.

Now, we define the generalized HFNs. A fuzzy subset $H\equiv \{(h_1,h_2,h_3,h_4,h_5,h_6)$; $u,u_L,u_R,{\mu }_H\}$, with $h_1\le h_2\le h_3\le h_4\le h_5\le h_6$, and $h_i\in \mathbb{R}$, for $i\in \{1,\cdots ,6\}$, where $u,u_L,u_R\in [0,1]$, is called a generalized HFN if its membership function is stated as

$${\mu }_H\left(x\right)=\left\{\begin{array}{cc}{\mu }_L^1\left(x\right),\hfill & \mathrm{if}{h}_1\le x<h_2;\hfill \\ {\mu }_L^2\left(x\right),\hfill & \mathrm{if}{h}_2\le x<h_3;\hfill \\ u,\hfill & \mathrm{if}{h}_3\le x<h_4;\hfill \\ {\mu }_R^2\left(x\right),\hfill & \mathrm{if}{h}_4\le x<h_5;\hfill \\ {\mu }_R^1\left(x\right),\hfill & \mathrm{if}{h}_5\le x\le h_6;\hfill \\ 0,\hfill & \mathrm{otherwise};\hfill \end{array}\right.$$

where ${\mu }_L^1\left(x\right)=u_L(\left(x-h_1\right)/\left(h_2-h_1\right))$, ${\mu }_L^2\left(x\right)=u_L+(u-u_L)(\left(x-h_2\right)/\left(h_3-h_2\right))$, ${\mu }_R^2\left(x\right)=u_R+(u-u_R)(\left(h_5-x\right)/\left(h_5-h_4\right))$, and ${\mu }_R^1\left(x\right)=u_R(\left(h_6-x\right)/\left(h_6-h_5\right))$ are the left lower, left upper, right upper, and right lower legs of H, respectively. A generalized HFN passes through $(h_1,0)$, $(h_2,u_L)$, $(h_3,u)$, $(h_4,u)$, $(h_5,u_R)$, and $(h_6,0)$. Hence, it is denoted by $H\equiv \{(h_1,h_2,h_3,h_4,h_5,h_6);u,u_L,u_R,{\mu }_H\}$.

2.3. Defuzzification

Defuzzification is the process of converting a fuzzy quantity into a precise (crisp) value. This process is crucial in practical applications where decisions need to be made based on fuzzy data. For instance, in a fuzzy inventory model, we may have fuzzy numbers representing uncertain parameters like holding cost, production rate, and deterioration rate. To make concrete decisions, such as determining optimal production levels or estimating costs, it is necessary to defuzzify these fuzzy numbers into specific values.

There are various methods of defuzzification, each with its own approach to handling the fuzziness of the data. Some common methods rely on the concept of h levels (also known as h cuts) for GMI. An h-level set of a fuzzy set ($\mathcal{F}$, where $h\in [0,1]$ is a real number representing the membership level used to define the cut) is a crisp set that contains all elements (x) for which ${\mu }_{\mathcal{F}}\left(x\right)\ge h$, and it is defined as ${\mathcal{F}}_h\equiv \{x\in X:{\mu }_{\mathcal{F}}\left(x\right)\ge h\}$.

For example, consider a TFN (D) representing the deterioration rate of a product defined as $D=(0.01,0.02,0.05,0.06)$. For a specific level h, say $h=0.5$, the cut h would be the set of values x where the membership function is ${\mu }_D\left(x\right)\ge 0.5$.

The GMI method deals with the weighted average of the membership values along the support of the fuzzy set using a specific weighted average of the cuts h. The two-step process is outlined as follows:

• Step 1—Determine the left and right boundaries of the fuzzy set at that level ($L^{-1}\left(h\right)$ and $R^{-1}\left(h\right)$, respectively), for each level h.
• Step 2—Compute a weighted average of the midpoint of the boundaries of Step 1, where each midpoint is weighted by its corresponding level h.

The continuous formula of the GMI method is given by

$$J\left(A\right)={\displaystyle \frac{{\int }_0^{w_A}h\left(\frac{L^{-1}\left(h\right)+R^{-1}\left(h\right)}{2}\right)\mathrm{d}h}{{\int }_0^{w_A}h\mathrm{d}h}},$$

(3)

where $L^{-1}\left(h\right)$ and $R^{-1}\left(h\right)$ are the inverse functions defining the left and right boundaries of the fuzzy set at level h, respectively, and $w_A$ is the width of the support of the fuzzy set (A), which can take any non-negative real value. The GMI method may be applied to various types of fuzzy numbers to obtain a crisp value. We now proceed with the application of the GMI method for TFNs, PFNs, and HFNs.

A TFN of $Z\equiv \left\{(y_1,y_2,y_3,y_4)\right\}$ is characterized by the following four points: $y_1$ (lower limit), $y_2$ (start of the plateau), $y_3$ (end of the plateau), and $y_4$ (upper limit). The membership function ${\mu }_Z$ increases linearly from $y_1$ to $y_2$, is constant from $y_2$ to $y_3$, and decreases linearly from $y_3$ to $y_4$. The GMI for a TFN is calculated as $J\left(Z\right)=\left(y_1+2y_2+2y_3+y_4\right)/6$ by considering the weighted average over the support of the fuzzy number. The weights correspond to the areas under the membership function, which can be divided into the following three regions: $y_1$ to $y_2$ with increasing membership, $y_2$ to $y_3$ with full membership, and $y_3$ to $y_4$ with decreasing membership. Points $y_2$ and $y_3$ are included twice in the weighted sum because they define the plateau where the membership is equal to one.

A PFN of $P\equiv \left\{(y_1,y_2,y_3,y_4,y_5)\right\}$ has the following five points: $y_1$ (lower limit), $y_2$ and $y_4$ (intermediate points), $y_3$ (peak), and $y_5$ (upper limit). The function ${\mu }_P$ increases linearly from $y_1$ to $y_2$, then from $y_2$ to $y_3$, whereas it decreases linearly from $y_3$ to $y_4$ and from $y_4$ to $y_5$. The GMI for a PFN is calculated as $J\left(P\right)=\left(y_1+3y_2+4y_3+3y_4+y_5\right)/12,$ where the weights correspond to the areas under the membership function as follows: $y_1$ to $y_2$ with increasing membership, $y_2$ to $y_3$ with increasing membership, $y_3$ to $y_4$ with decreasing membership, and $y_4$ to $y_5$ with decreasing membership. Points $y_2$ and $y_4$ have weights of three each, and the peak at $y_3$ corresponds to a weight of four because it is the maximum membership.

An HFN of $H\equiv \left\{(y_1,y_2,y_3,y_4,y_5,y_6)\right\}$ consists of the following six points: $y_1$ (lower limit), $y_2$ and $y_5$ (intermediate points), $y_3$ and $y_4$ (peak values), and $y_6$ (upper limit). The membership function ${\mu }_H$ is piecewise linear, with linear increases and decreases in membership. To derive the GMI for HFNs, we need to consider the weighted average over the support of the fuzzy number by applying the expression stated in (3). This involves calculating the areas under the membership function for different segments as follows. From $y_1$ to $y_2$ (increasing membership), we have ${\mu }_H\left(x\right)=\left(x-y_1\right)/\left(y_2-y_1\right)$ and therefore $L^{-1}\left(h\right)=y_1+h(y_2-y_1)$. From $y_2$ to $y_3$ (increasing membership), we have ${\mu }_H\left(x\right)=\left(x-y_2\right)/\left(y_3-y_2\right)$ and hence $L^{-1}\left(h\right)=y_2+h(y_3-y_2)$. From $y_3$ to $y_4$ (constant membership), we have ${\mu }_H\left(x\right)=1$$(\mathrm{maximum}\mathrm{membership})$, with $L^{-1}\left(h\right)$ and $R^{-1}\left(h\right)$ being the constants ($y_3$ and $y_4$). From $y_4$ to $y_5$ (decreasing membership), we have ${\mu }_H\left(x\right)=\left(y_5-x\right)/\left(y_5-y_4\right)$ and therefore $R^{-1}\left(h\right)=y_5-h(y_5-y_4)$. From $y_5$ to $y_6$ (decreasing membership), we have ${\mu }_H\left(x\right)=\left(y_6-x\right)/(y_6-y_5$ and thus $R^{-1}\left(h\right)=y_6-h(y_6-y_5)$.

By integrating these functions along h from zero to one, we obtain the weighted average as $J\left(H\right)=\left({\int }_0^{1}h(\left(L^{-1}\left(h\right)+R^{-1}\left(h\right)\right)/2)\mathrm{d}h\right)/{\int }_0^{1}h\mathrm{d}h.$ For each segment, the integral is calculated separately. From $y_1$ to $y_2$, we have ${\int }_0^{1}h(\left(y_1+h(y_2-y_1)+y_6-h(y_6-y_5)\right)/2)\mathrm{d}h.$ From $y_2$ to $y_3$, we have ${\int }_0^{1}h(\left(y_2+h(y_3-y_2)+y_5-h(y_5-y_4)\right)/2)\mathrm{d}h.$ From $y_3$ to $y_4$, we reach ${\int }_0^{1}h(\left(y_3+y_4\right)/2)\mathrm{d}h.$ From $y_4$ to $y_5$, we attain ${\int }_0^{1}h(\left(y_4+y_5-h(y_5-y_4)\right)/2)\mathrm{d}h.$ From $y_5$ to $y_6$, we obtain ${\int }_0^{1}h(\left(y_5+y_6-h(y_6-y_5)\right)/2)\mathrm{d}h.$ By combining these integrals, we derive the GMI formula for the HFN as $J\left(H\right)=\left(y_1+3y_2+2y_3+2y_4+3y_5+y_6\right)/12.$ Each weight corresponds to the area under the membership function, making points $y_2$ and $y_5$ (weights of three each) and peaks $y_3$ and $y_4$ (weights of two each) more relevant due to the maximum membership.

The earlier calculations illustrate how the GMI method is applied to different fuzzy numbers to obtain a crisp value that can be used for decision making in practical applications.

3. Cost Optimization and Solution Method

This section outlines the mathematical formulation of the TC and the optimization methods used to solve the inventory model. The Lagrangian method is applied to find the optimal solution, incorporating uncertainties through the use of fuzzy variables.

3.1. Mathematical Formulation of Inventory Total Cost

Consider a manufacturing inventory model with improved reliability and trade credit terms [50]. The producer fabricates non-defective items during the manufacturing process. The inventory starts with zero items. When customer’s demand is met, the manufacturer begins production at rate P at time $t=0$, with $0<t<T$ and T being the production cycle duration representing the total time over which production and inventory processes occur.

The TC of the goods includes ordering cost, carrying cost, deterioration cost, interest paid, and interest earned. The goal is to minimize the TC function expressed in (1). The production halt time ($t_1$), representing the time when production is stopped, is constrained by the expression defined in (2). We determine the optimal cycle period (L) by minimizing the TC while satisfying all constraints. To find the optimal cycle period (L), we differentiate $\mathrm{TC}\left(L\right)$ with respect to L and equate it to zero, yielding

$${\displaystyle \frac{R_0}{T^2}}+{\displaystyle \frac{f{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}^2}}-{\displaystyle \frac{(o^2-p^2(1-\vartheta ))Mh{J}_e}{2T^2}}={\displaystyle \frac{(d+⌀f+f{J}_c){\sigma }_{e}G}{⌀T^2}}\left(L{\displaystyle \frac{t_1}{T}}-t_1\right).$$

(4)

By substituting the expression for $t_1$ from (2) into (4), we reach

$$R_0+{\displaystyle \frac{f{J}_c(1-⌀o-e^{-⌀o})\left({\sigma }_{e}G\right)}{{\theta }^2}}-{\displaystyle \frac{(o^2-p^2(1-\vartheta ))Mh{J}_e}{2}}={\displaystyle \frac{(d+⌀f+f{J}_c)}{2⌀}}\left(h⌀-{\displaystyle \frac{h^2⌀}{{\sigma }_{e}G}}\right)L^2.$$

Thus, solving for L, we arrive at

$$L=\sqrt{{\displaystyle \frac{2\left(R_0+\frac{f{J}_c(1-⌀o-e^{-⌀o})\left({\sigma }_{e}G\right)-h}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))Mh{J}_e}{2}\right)}{\left(\frac{d+⌀f+f{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$

The manufacturing company halts production at time $t=T$, ensuring that all customer’s demand is satisfied by the goods produced during the cycle.

3.2. Application of the Kuhn–Tucker Conditions in Optimization

To solve the optimization problem, we employ the Kuhn–Tucker conditions, which ensure the optimality of solutions under both equality and inequality constraints. These conditions include the following:

• Primal feasibility—The solution must satisfy all the original constraints of the optimization problem.
• Dual feasibility—The Lagrange multipliers associated with inequality constraints must be non-negative.
• Complementary slackness—For each inequality constraint, either the constraint is active with a positive multiplier or it is inactive and its multiplier is zero.
• Stationarity—The gradient of the Lagrangian function, incorporating the objective function and constraints, must be zero at the optimal point.

When dealing with fuzzy variables, fuzzy arithmetic is applied to both the constraints and objective function. This arithmetic expands the feasible solution space and allows the optimization process to account for uncertainty, seeking robust solutions under varying fuzzy parameters.

3.3. Fuzzy Inventory Model for the Total Cost

In the fuzzy inventory model, we convert the variables d, $R_0$, f, T, and M from the model formulated in (1) into fuzzy variables, which can be represented as TFNs, PFNs, or HFNs. This fuzzification introduces uncertainty into key parameters such as demand (d), cost components ($R_0$, f, T, and M), and production reliability (${\sigma }_e$).

The Kuhn–Tucker conditions, as discussed in Section 3.2, are then applied to ensure that the solution satisfies both equality and inequality constraints. Using the GMI method, we defuzzify the fuzzy values and facilitate optimization, ensuring that the solution is robust under uncertainty. The TC function for the fuzzy inventory system is expressed as

$$\mathrm{TC}\left(L\right)={\displaystyle \frac{1}{6}}\left(\begin{array}{c}\hfill \left({\displaystyle \frac{R_{01}}{T}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T}}+{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^{2}T}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T}}\right)\\ \hfill +2\left({\displaystyle \frac{R_{02}}{T}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T}}+{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^{2}T}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T}}\right)\\ \hfill +2\left({\displaystyle \frac{R_{03}}{T}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T}}+{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^{2}T}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T}}\right)\\ \hfill +\left({\displaystyle \frac{R_{04}}{T}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T}}+{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^{2}T}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T}}\right)\end{array}\right).$$

(5)

To find the optimal value of T, we partially differentiate $\mathrm{TC}\left(L\right)$, as defined in (5), with respect to L, following the stationarity condition of the Kuhn–Tucker optimality criteria, reaching

$${\displaystyle \frac{\partial \mathrm{TC}\left(L\right)}{\partial L}}={\displaystyle \frac{1}{6}}\left(\begin{array}{c}\hfill \left(-{\displaystyle \frac{R_{01}}{T^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T^2}}{\sigma }_{e}G\left(L{\displaystyle \frac{t_1}{T}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T^2}}\right)\\ \hfill +2\left(-{\displaystyle \frac{R_{02}}{T^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T^2}}{\sigma }_{e}G\left(L{\displaystyle \frac{t_1}{T}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T^2}}\right)\\ \hfill +2\left(-{\displaystyle \frac{R_{03}}{T^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T^2}}{\sigma }_{e}G\left(L{\displaystyle \frac{t_1}{T}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T^2}}\right)\\ \hfill +\left(-{\displaystyle \frac{R_{04}}{T^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T^2}}{\sigma }_{e}G\left(L{\displaystyle \frac{t_1}{T}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T^2}}\right)\end{array}\right),$$

Hence, by solving the equation obtained when setting this derivative equal to zero, we determine the optimal value of T as given by

$$T=\sqrt{{\displaystyle \frac{2\left(\begin{array}{c}\hfill (R_{01}+2R_{02}+2R_{03}+R_{04})+\frac{(1-⌀o-e^{-⌀o})J_c(f_1+2f_2+2f_3+f_4)({\sigma }_{e}G-h)}{{⌀}^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))h{J}_e(M_1+2M_2+2M_3+M_4)}{2}\end{array}\right)}{\left(\frac{d_1+2d_2+2d_3+d_4+⌀(f_1+2f_2+2f_3+f_4)+J_c(f_1+2f_2+2f_3+f_4)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$

3.4. Lagrangian Method for Solving EOQ Models with Fuzzy Variables

We now explore the Lagrangian method commonly used for optimizing inventory models by finding the optimal solution that minimizes costs and maximizes efficiency [48,49,51]. In this case, the Kuhn–Tucker conditions are used to ensure that the optimal solution respects both the equality and inequality constraints. Specifically, the dual feasibility and complementary slackness conditions ensure that the Lagrange multipliers correctly reflect the impact of the fuzzy constraints in the optimization process.

By applying the Lagrangian method, we determine the average cost for the manufactured product under uncertainty conditions. The procedure involves the following steps:

• Step 1—Define the TC function stated in (5), including the ordering cost, carrying cost, deterioration cost, interest paid, and interest earned.
• Step 2—Convert crisp variables to fuzzy variables to account for uncertainty.
• Step 3—Apply the GMI method to calculate the average cost over subperiods, adjusting the number of subperiods according to the type of fuzzy numbers defined in the previous step (for example, four subperiods for TFNs, five for PFNs, and six for HFNs). The TC function presented in (5) is correspondingly modified for each case.
• Step 4—Use the Lagrangian function, which integrates the TC minimization objective with the problem constraints (such as production time limits and credit terms) while accounting for the division of the production cycle into subperiods.
• Step 5—Solve the system of equations obtained from the partial derivatives of the Lagrangian function across the subperiods, ensuring that the solution satisfies both the equality and inequality constraints to find the optimal production cycle period.
• Step 6—Verify that the solution satisfies all constraints, including the appropriate number of subperiods (4, 5, or 6) depending on the type of fuzzy numbers used.

By following these steps, as illustrated in Figure 1, we ensure that the optimal cycle period is found, satisfying all constraints and minimizing the TC.

3.5. Case of Trapezoidal Fuzzy Numbers

Next, we apply the Lagrangian method to the specific case of TFNs. Let $d=(d_1,d_2,d_3,d_4)$, $R_0=(R_{01},R_{02},R_{03},R_{04})$, $f=(f_1,f_2,f_3,f_4)$, $T=(T_1,T_2,T_3,T_4)$, and $M=(M_1,M_2,M_3,M_4)$.

In the fuzzy model formulated in (5), the cycle period (T) is divided into subperiods $T_1,T_2,T_3$, and $T_4$. The function $\mathrm{TC}\left(L\right)$ for the TFNs is expressed as

$$\mathrm{TC}\left(L\right)={\displaystyle \frac{1}{6}}\left(\begin{array}{c}\hfill \left({\displaystyle \frac{R_{01}}{T_4}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_4}}+{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_4}}\right)\\ \hfill +2\left({\displaystyle \frac{R_{02}}{T_3}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_3}}+{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_3}}\right)\\ \hfill +2\left({\displaystyle \frac{R_{03}}{T_2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_2}}+{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_2}}\right)\\ \hfill +\left({\displaystyle \frac{R_{04}}{T_1}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_1}}+{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_1}}\right)\end{array}\right),$$

where $0<T_1\le T_2\le T_3\le T_4$ and the cost components (ordering cost, carrying cost, deterioration cost, interest paid, and interest earned) are calculated for each subperiod. These inequalities can be rewritten as $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, and $T_1>0$. The optimal cycle length ($T^{*}$) is then calculated as follows.

To locate the minimum, we employ the Lagrangian method according to the following four steps:

• Step 1—Find the minimum of $\mathrm{TC}\left(L\right)$ from

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_1}}=0\Rightarrow T_1=\sqrt{{\displaystyle \frac{R_{04}+\frac{(1-⌀o-e^{-⌀o})f_4{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_4{J}_e}{2}}{\left(\frac{d_4+⌀f_4+f_4{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_2}}=0\Rightarrow T_2=\sqrt{2\left({\displaystyle \frac{R_{03}+\frac{(1-⌀o-e^{-⌀o})f_3{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_3{J}_e}{2}}{\left(\frac{d_3+⌀f_3+f_3{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}\right)},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_3}}=0\Rightarrow T_3=\sqrt{2\left({\displaystyle \frac{R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}\right)},$$

  and

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_4}}=0\Rightarrow T_4=\sqrt{{\displaystyle \frac{R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  In the above expressions, we observe that $T_1>T_2>T_3>T_4$, which contradicts the required condition of $0<T_1\le T_2\le T_3\le T_4$. This inconsistency indicates that the constraints or the method must be adjusted to ensure that the solution satisfies the required ordering of the cycle periods. Possible adjustments include revising the Lagrange multipliers or tightening the conditions imposed on the cycle lengths.
• Step 2—Reformulate the constraint of $T_2-T_1\ge 0$ as $T_2-T_1=0$ and next apply the Lagrangian function, $N(T_1,T_2,T_3,T_4,\tau )=\mathrm{TC}\left(L\right)-\tau (T_2-T_1)$ say, where $\tau$ is the Lagrange multiplier. By taking the partial derivatives of N with respect to the cycle periods $T_1,T_2,T_3,T_4$ and the Lagrange multiplier $\tau$, we get the system of equations given by

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_4^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{G}_c)}{\theta T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_4^2}}\right)+\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{02}}{T_3^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{G}_c)}{\theta T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_3^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_2^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{G}_c)}{\theta T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_2^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{04}}{T_1^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{G}_c)}{\theta T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_1^2}}\right)=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial \tau }}=-(T_2-T_1)=0,$$

  obtaining

  $$T_1=T_2=\sqrt{{\displaystyle \frac{2\left((R_{04}+2R_{03})+\frac{(1-⌀o-e^{-⌀o})(f_4{J}_c+2f_3{J}_c)({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h(M_4{J}_e+2M_3{J}_e)}{2}\right)}{\left(\frac{d_4+2d_3+⌀(f_4+2f_3)+(f_4{J}_c+2f_4{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_3=\sqrt{{\displaystyle \frac{2\left(R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}\right)}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $$T_4=\sqrt{{\displaystyle \frac{2\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}$$
  Based on this, we observe that $T_3>T_4$, which indicates that further adjustments to the constraints are required to ensure the correct order of $T_1\le T_2\le T_3\le T_4$. These adjustments may involve revising the Lagrange multipliers or modifying the conditions imposed on the system to guarantee compliance with the required sequence.
• Step 3—Convert the constraints of $T_2-T_1\ge 0$ and $T_3-T_2\ge 0$ into $T_2-T_1=0$ and $T_3-T_2=0$, respectively and then optimize $\mathrm{TC}\left(L\right)$ using the Lagrangian method as

  $$N(T_1,T_2,T_3,T_4,{\tau }_1,{\tau }_2)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2),$$

  where ${\tau }_1$ and ${\tau }_2$ are the Lagrange multipliers, with the partial derivatives of N being formulated as

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_4^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_4^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{02}}{T_3^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_3^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_2^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_2^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_1^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_1^2}}\right)=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2).$$
  From the above expressions, we obtain

  $$T_1=T_2=T_3=\sqrt{{\displaystyle \frac{2\left((R_{04}+2R_{03}+2R_{02})+\frac{(1-⌀o-e^{-⌀o})(f_4{J}_c+2f_3{J}_c+2f_2{J}_c)({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))(M_4{J}_e+2M_3{J}_e+2M_2{J}_e)}{2}\right)}{\left(\frac{d_2+2d_3+2d_2+⌀(f_4+2f_3+2f_2)+(f_4{J}_c+2f_3{J}_c+2f_2{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}$$

  and

  $$T_4=\sqrt{{\displaystyle \frac{2\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  indicating that $T_4>T_1$. Therefore, further adjustments to the constraints are needed to satisfy the required order of cycle periods.
• Step 4—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, and $T_4-T_3\ge 0$ into $T_2-T_1=0$, $T_3-T_2=0$, and $T_4-T_3=0$, respectively, and use the Lagrangian function given by

  $$N(T_1,T_2,T_3,T_4,{\tau }_1,{\tau }_2,{\tau }_3)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3),$$

  where ${\tau }_1$, ${\tau }_2$, and ${\tau }_3$ are the Lagrange multipliers, with the partial derivatives of N being presented as

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_4^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_4^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{02}}{T_3^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_3^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_2^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_2^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_1^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_1^2}}\right)-{\tau }_3=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3).$$
  Then, the cycle periods are given by

  $$T_1=T_2=T_3=T_4=\sqrt{{\displaystyle \frac{2\left((R_{04}+2R_{03}+2R_{02}+R_{01})+\frac{(1-⌀o-e^{-⌀o})J_c(f_4+2f_3+2f_2+f_1)({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{J}_e(M_4+2M_3+2M_2+M_1)}{2}\right)}{\left(\frac{d_4+2d_3+2d_2+d_1+⌀(f_4+2f_3+2f_2+f_1)+J_c(f_4+2f_3+2f_2+f_1)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  Therefore, the solution of $T^{*}=(T_1,T_2,T_3,T_4)$ satisfies all inequality constraints. Thus, we obtain the TC using the nonlinear Lagrangian method for the fuzzy inventory model defined in (5) for the case of TFNs. To summarize and visualize these steps, we present Algorithm 1 and the corresponding flow chart in Figure 2 for better visualization.

Algorithm 1: Lagrangian method for solving EOQ models.

3.6. Case of Pentagonal Fuzzy Numbers

We now extend the fuzzy inventory model TC to the case of PFNs. Unlike the trapezoidal case, where the cycle period (T) was divided into four subintervals, in the pentagonal case, the cycle length is divided into five subintervals. This affects the overall optimization process and allows for a more granular modeling of uncertainty. The relevant crisp variables are converted into fuzzy variables represented by PFNs, including the cycle length $T=(T_1,T_2,T_3,T_4,T_5)$, the demand rate $d=(d_1,d_2,d_3,d_4,d_5)$, the setup cost $R_0=(R_{01},R_{02},R_{03},R_{04},R_{05})$, the holding cost $\mathcal{F}=(f_1,f_2,f_3,f_4,f_5)$, and the manufacturing cost $M=(M_1,M_2,M_3,M_4,M_5)$.

Applying the GMI method to these variables, we reach the expression for the TC stated as

$$\mathrm{TC}\left(L\right)={\displaystyle \frac{1}{6}}\left(\begin{array}{c}\hfill \left({\displaystyle \frac{R_{01}}{T_5}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_5}}+{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_5}}\right)\\ \hfill +3\left({\displaystyle \frac{R_{02}}{T_4}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_4}}+{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_4}}\right)\\ \hfill +4\left({\displaystyle \frac{R_{03}}{T_3}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_3}}+{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_3}}\right)\\ \hfill +3\left({\displaystyle \frac{R_{04}}{T_2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_2}}+{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2}}\right)\\ \hfill +\left({\displaystyle \frac{R_{05}}{T_1}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)({\sigma }_{e}G{t}_1-hL)}{⌀T_1}}+{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_1}}\right)\end{array}\right),$$

where $0<T_1\le T_2\le T_3\le T_4\le T_5$. The inequality conditions $0<T_1\le T_2\le T_3\le T_4\le T_5$ can be expressed as $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, $T_5-T_4\ge 0$, and $T_1>0$. To determine the minimum, we employ the Lagrange method considering the following five steps:

• Step 1—Consider the results given by

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_1}}=0\Rightarrow T_1=\sqrt{{\displaystyle \frac{R_{05}+\frac{(1-⌀o-e^{-⌀o})f_5{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_5{J}_e}{2}}{\left(\frac{d_5+⌀f_5+f_5{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_2}}=0,\Rightarrow T_2=\sqrt{3{\displaystyle \frac{R_{04}+\frac{(1-⌀o-e^{-⌀o})f_4{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_4{J}_e}{2}}{\left(\frac{d_4+⌀f_4+f_4{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_3}}=0,\Rightarrow T_3=\sqrt{4{\displaystyle \frac{R_{03}+\frac{(1-⌀o-e^{-⌀o})f_3{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_3{J}_e}{2}}{\left(\frac{d_3+⌀f_3+f_3{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_4}}=0,\Rightarrow T_4=\sqrt{3{\displaystyle \frac{R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_5}}=0,\Rightarrow T_5=\sqrt{{\displaystyle \frac{R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above results, we obtain $T_1>T_2>T_3>T_4>T_5$, which does not satisfy the required condition of $0<T_1\le T_2\le T_3\le T_4\le T_5$. Therefore, we must adjust the constraints to ensure that the solution satisfies the required ordering. Possible adjustments include revising the Lagrange multipliers or introducing stricter conditions for the cycle lengths to guarantee the required sequence.
• Step 2—Convert the constraint of $T_2-T_1\ge 0$ into $T_2-T_1=0$ and use the Lagrangian function given by $N(T_1,T_2,T_3,T_4,T_5,\tau )=\mathrm{TC}\left(L\right)-\tau (T_2-T_1)$, where $\tau$ is the Lagrange multiplier, to obtain
  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_5^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{G}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_5^2}}\right)+\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_4^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{G}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_4^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{4}{6}}\left(-{\displaystyle \frac{R_{03}}{T_3^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{G}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_3^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{04}}{T_2^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{G}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{05}}{T_1^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{G}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_1^2}}\right)=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial \tau }}=-(T_2-T_1)=0,$$

  reaching

  $$T_1=T_2=\sqrt{{\displaystyle \frac{3\left((R_{05}+2R_{04})+\frac{(1-⌀o-e^{-⌀o})(f_5{J}_c+2f_4{J}_c)({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h(M_5{J}_e+2M_4{J}_e)}{2}\right)}{\left(\frac{d_5+2d_4+⌀(f_5+2f_4)+(f_5{J}_c+2f_4{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_3=\sqrt{{\displaystyle \frac{4\left(R_{03}+\frac{(1-⌀o-e^{-⌀o})f_3{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_3{J}_e}{2}\right)}{\left(\frac{d_3+⌀f_3+f_3{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_4=\sqrt{{\displaystyle \frac{3\left(R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}\right)}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_4>T_5$, which violates the required condition of $0<T_1\le T_2\le T_3\le T_4\le T_5$. Thus, it is necessary to further adjust the constraints to ensure that the solution meets the specified conditions.
• Step 3—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, and $T_5-T_4\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, $T_4-T_3=0$, and $T_5-T_4=0$, respectively, and apply the Lagrangian function defined by

  $$N(T_1,T_2,T_3,T_4,T_5,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4),$$

  where ${\tau }_1$, ${\tau }_2$, ${\tau }_3$, and ${\tau }_4$ are the Lagrange multipliers used, to obtain

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_5^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_5^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_4^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_4^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{4}{6}}\left(-{\displaystyle \frac{R_{03}}{T_3^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_3^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{04}}{T_2^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{05}}{T_1^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_1^2}}\right)-{\tau }_4=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4).\left)\right)$$
  Then, solving these equations, we get

  $$T_1=T_2=T_3=\sqrt{{\displaystyle \frac{3\left((R_{05}+3R_{04}+4R_{03})+\frac{(1-⌀o-e^{-⌀o})(f_5{J}_c+3f_4{J}_c+4f_3{J}_c)({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))(M_5{J}_e+3M_4{J}_e+4M_3{J}_e)}{2}\right)}{\left(\frac{d_5+3d_4+4d_3+⌀(f_5+3f_4+4f_3)+(f_5{J}_c+3f_4{J}_c+4f_3{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_4=\sqrt{{\displaystyle \frac{3\left(R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}\right)}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we find that $T_4>T_5$, which indicates a violation of the required condition of $0<T_1\le T_2\le T_3\le T_4\le T_5$. Therefore, further adjustments to the constraints are necessary to ensure that the solution meets the required conditions.
• Step 4—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, and $T_5-T_4\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, $T_4-T_3=0$, and $T_5-T_4=0$, respectively, and employ the Lagrangian function given by

  $$N(T_1,T_2,T_3,T_4,T_5,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4)$$

  where ${\tau }_1$, ${\tau }_2$, ${\tau }_3$, and ${\tau }_4$ are the Lagrange multipliers. Taking the partial derivatives of N with respect to the variables $T_1,T_2,T_3,T_4,T_5$ and the multipliers ${\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4$, we obtain the system of equations stated as

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_5^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_5^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_4^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_4^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{4}{6}}\left(-{\displaystyle \frac{R_{03}}{T_3^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_3^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{04}}{T_2^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{05}}{T_1^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_1^2}}\right)-{\tau }_4=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4).$$
  Then, we obtain
  $$T_1=T_2=T_3=T_4=\sqrt{{\displaystyle \frac{3\left((R_{05}+3R_{04}+4R_{03}+3R_{02})+\frac{(1-⌀o-e^{-⌀o})(f_5{J}_c+3f_4{J}_c+4f_3{J}_c+3f_2{J}_c)({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))(M_5{J}_e+3M_4{J}_e+4M_3{J}_e+3M_2{J}_e)}{2}\right)}{\left({\displaystyle \frac{d_5+3d_4+4d_3+3d_2+⌀(f_5+3f_4+4f_3+3f_2)+(f_5{J}_c+3f_4{J}_c+4f_3{J}_c+3f_2{J}_c)}{⌀}}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and
  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_5>T_1$, which contradicts the requirement of $0<T_1\le T_2\le T_3\le T_4\le T_5$. Therefore, we need to further adjust the constraints to ensure that the solution meets the required conditions.
• Step 5—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, and $T_4-T_3\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, and $T_4-T_3=0$, respectively, and optimize $\mathrm{TC}\left(L\right)$ using the Lagrangian method as

  $$N(T_1,T_2,T_3,T_4,T_5,{\tau }_1,{\tau }_2,{\tau }_3)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3).$$

  generating

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_5^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_5^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_4^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_4^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{4}{6}}\left(-{\displaystyle \frac{R_{03}}{T_3^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_3^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{04}}{T_2^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{05}}{T_1^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_1^2}}\right)=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),$$

  arriving at

  $$T_1=T_2=T_3=T_4=T_5=\sqrt{{\displaystyle \frac{3\left(\begin{array}{c}\hfill (R_{05}+3R_{04}+4R_{03}+3R_{02}+R_{01})+\frac{(1-⌀o-e^{-⌀o})J_c(f_5+3f_4+4f_3+3f_2+f_1)({\sigma }_{e}G-h)}{{⌀}^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))h{J}_e(M_5+3M_4+4M_3+3M_2+M_1)}{2}\end{array}\right)}{\left({\displaystyle \frac{d_5+3d_4+4d_3+3d_2+d_1+⌀(f_5+3f_4+4f_3+3f_2+f_1)+J_c(f_5+3f_4+4f_3+3f_2+f_1)}{⌀}}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}$$

  and

  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_5>T_1$, which violates the required condition of $0<T_1\le T_2\le T_3\le T_4\le T_5$. Consequently, it is necessary to further adjust the constraints to ensure that the solution meets the imposed conditions.
• Step 6—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, and $T_5-T_4\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, $T_4-T_3=0$, and $T_5-T_4=0$, respectively, with the corresponding Lagrangian function being expressed as

  $$N(T_1,T_2,T_3,T_4,T_5,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4),$$

  where ${\tau }_1$, ${\tau }_2$, ${\tau }_3$, and ${\tau }_4$ are the Lagrange multipliers. By taking the partial derivatives of N with respect to the cycle periods ($T_1,T_2,T_3,T_4,T_5$) and the Lagrange multipliers (${\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4$), we obtain
  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_5^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_5^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_4^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_4^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{4}{6}}\left(-{\displaystyle \frac{R_{03}}{T_3^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_3^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{04}}{T_2^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{05}}{T_1^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_1^2}}\right)-{\tau }_4=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4),$$

  reaching

  $$T_1=T_2=T_3=T_4=T_5=\sqrt{{\displaystyle \frac{3\left(\begin{array}{c}\hfill (R_{05}+3R_{04}+4R_{03}+3R_{02}+R_{01})+\frac{(1-⌀o-e^{-⌀o})J_c(f_5+3f_4+4f_3+3f_2+f_1)({\sigma }_{e}G-h)}{{⌀}^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))h{J}_e(M_5+3M_4+4M_3+3M_2+M_1)}{2}\end{array}\right)}{\left({\displaystyle \frac{d_5+3d_4+4d_3+3d_2+d_1+⌀(f_5+3f_4+4f_3+3f_2+f_1)+J_c(f_5+3f_4+4f_3+3f_2+f_1)}{⌀}}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$
  Thus, the solution of $T^{*}=(T_1,T_2,T_3,T_4,T_5)$ satisfies all the inequality constraints, ensuring an optimal cycle period.

3.7. Case of Hexagonal Fuzzy Numbers

Next, we extend the fuzzy inventory model to the case of HFNs. Similar to the trapezoidal and pentagonal cases, the relevant crisp variables are converted into fuzzy variables represented by HFNs. However, in this case, the cycle length (T) is divided into six subintervals, allowing for a finer and more detailed representation of uncertainty and variability. The variables involved in this division include the cycle length $T=(T_1,T_2,T_3,T_4,T_5,T_6)$, the demand rate $d=(d_1,d_2,d_3,d_4,d_5,d_6)$, the setup cost $R_0=(R_{01},R_{02},R_{03},R_{04},R_{05},R_{06})$, the holding cost $\mathcal{F}=(f_1,f_2,f_3,f_4,f_5,f_6)$, and the manufacturing cost $M=(M_1,M_2,M_3,M_4,M_5,M_6)$.

The TC for the manufactured product using HFNs and the GMI method is expressed as

$$\mathrm{TC}\left(L\right)={\displaystyle \frac{1}{6}}\left(\begin{array}{c}\left({\displaystyle \frac{R_{01}}{T_6}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)({\sigma }_{e}G{t}_1-h{L}_6)}{⌀T_6}}+{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_6}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_6}}\right)\hfill \\ +3\left({\displaystyle \frac{R_{02}}{T_5}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)({\sigma }_{e}G{t}_1-h{L}_5)}{⌀T_5}}+{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_5}}\right)\hfill \\ +2\left({\displaystyle \frac{R_{03}}{T_4}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)({\sigma }_{e}G{t}_1-h{L}_4)}{⌀T_4}}+{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_4}}\right)\hfill \\ +2\left({\displaystyle \frac{R_{04}}{T_3}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)({\sigma }_{e}G{t}_1-h{L}_3)}{⌀T_3}}+{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_3}}\right)\hfill \\ +3\left({\displaystyle \frac{R_{05}}{T_2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)({\sigma }_{e}G{t}_1-h{L}_2)}{⌀T_2}}+{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_2}}\right)\hfill \\ +\left({\displaystyle \frac{R_{06}}{T_1}}+{\displaystyle \frac{(d_6+⌀f_6+f_6{J}_c)({\sigma }_{e}G{t}_1-h{L}_1)}{⌀T_1}}+{\displaystyle \frac{f_6{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1}}-{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{6}h{J}_e}{2T_1}}\right)\hfill \end{array}\right),$$

where $0<T_1\le T_2\le T_3\le T_4\le T_5\le T_6$. These inequality conditions can be reformulated as $T_2-T_1\ge 0,T_3-T_2\ge 0,T_4-T_3\ge 0,T_5-T_4\ge 0,T_6-T_5\ge 0,T_1>0$. To determine the minimum value of the TC, we apply the Lagrangian method according to the following six steps:

• Step 1—Find the minimum of $\mathrm{TC}\left(\mathrm{L}\right)$ from

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_1}}=0\Rightarrow T_1=\sqrt{{\displaystyle \frac{R_{06}+\frac{(1-⌀o-e^{-⌀o})f_6{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_6{J}_e}{2}}{\left(\frac{d_6+⌀f_6+f_6{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_2}}=0\Rightarrow T_2=\sqrt{3{\displaystyle \frac{R_{05}+\frac{(1-⌀o-e^{-⌀o})f_5{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_5{J}_e}{2}}{\left(\frac{d_5+⌀f_5+f_5{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_3}}=0\Rightarrow T_3=\sqrt{2{\displaystyle \frac{R_{04}+\frac{(1-⌀o-e^{-⌀o})f_4{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_4{J}_e}{2}}{\left(\frac{d_4+⌀f_4+f_4{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_4}}=0\Rightarrow T_4=\sqrt{2{\displaystyle \frac{R_{03}+\frac{(1-⌀o-e^{-⌀o})f_3{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_3{J}_e}{2}}{\left(\frac{d_3+⌀f_3+f_3{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_5}}=0\Rightarrow T_5=\sqrt{3{\displaystyle \frac{R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $${\displaystyle \frac{\partial \mathrm{TC}}{\partial T_6}}=0\Rightarrow T_6=\sqrt{{\displaystyle \frac{R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_1>T_2>T_3>T_4>T_5>T_6$, which does not satisfy the required condition of $0<T_1\le T_2\le T_3\le T_4\le T_5\le T_6$.
• Step 2—Convert the constraint of $T_2-T_1\ge 0$ into the equality of $T_2-T_1=0$, with the corresponding Lagrangian function being defined as

  $$N(T_1,T_2,T_3,T_4,T_5,T_6,\tau )=\mathrm{TC}\left(L\right)-\tau (T_2-T_1),$$

  where $\tau$ is the Lagrange multiplier, to generate
  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_6^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{G}_c)}{\theta T_6^2}}{\sigma }_{e}G\left(T_6{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_6^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_6^2}}\right)+\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_5^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{G}_c)}{\theta T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_5^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_4^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{G}_c)}{\theta T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_4^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_3^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{G}_c)}{\theta T_3^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_3^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{04}}{T_2^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{G}_c)}{\theta T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_2^2}}\right)-\tau =0,$$

  $${\displaystyle \frac{\partial N}{\partial T_6}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{06}}{T_1^2}}+{\displaystyle \frac{(d_6+⌀f_6+f_6{G}_c)}{\theta T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_6{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{6}h{J}_e}{2T_1^2}}\right)=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial \tau }}=-(T_2-T_1)=0,$$

  obtaining

  $$T_1=T_2=\sqrt{{\displaystyle \frac{\left((R_{06}+2R_{05})+\frac{(1-⌀o-e^{-⌀o})(f_6{J}_c+2f_5{J}_c)({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h(M_6{J}_e+2M_5{J}_e)}{2}\right)}{\left(\frac{d_6+2d_5+⌀(f_6+2f_5)+(f_6{J}_c+2f_5{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_3=\sqrt{{\displaystyle \frac{3\left(R_{04}+\frac{(1-⌀o-e^{-⌀o})f_4{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_4{J}_e}{2}\right)}{\left(\frac{d_4+⌀f_4+f_4{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_4=\sqrt{{\displaystyle \frac{2\left(R_{03}+\frac{(1-⌀o-e^{-⌀o})f_3{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_3{J}_e}{2}\right)}{\left(\frac{d_3+⌀f_3+f_3{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}\right)}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $$T_6=\sqrt{{\displaystyle \frac{3\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_5>T_6$, which does not satisfy the condition of $0<T_1\le T_2\le T_3\le T_4\le T_5\le T_6$.
• Step 3—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, $T_5-T_4\ge 0$, and $T_6-T_5\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, $T_4-T_3=0$, $T_5-T_4=0$, and $T_6-T_5=0$, respectively, with the corresponding Lagrangian function being given by

  $$N(T_1,T_2,T_3,T_4,T_5,T_6,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4)-{\tau }_5(T_6-T_5),$$

  where ${\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4$, and ${\tau }_5$ are the Lagrange multipliers, attaining at
  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_6^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_6^2}}{\sigma }_{e}G\left(T_6{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_6^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_6^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_5^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_5^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_4^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_4^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{04}}{T_3^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_3^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{05}}{T_2^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_2^2}}\right)-{\tau }_4=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_6}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{06}}{T_1^2}}+{\displaystyle \frac{(d_1+⌀f_6+f_6{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_6{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{6}h{J}_e}{2T_1^2}}\right)-{\tau }_5=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4),$$

  generating

  $$T_1=T_2=T_3=\sqrt{{\displaystyle \frac{3\left(\begin{array}{c}\hfill (R_{06}+3R_{05}+2R_{04})+\frac{(1-⌀o-e^{-⌀o})(f_6{J}_c+3f_5{J}_c+2f_4{J}_c)({\sigma }_{e}G-h)}{{\theta }^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))(M_6+3M_5{J}_e+2M_4{J}_e)}{2}\end{array}\right)}{\left(\frac{d_6+3d_5+2d_4+⌀(f_6+3f_5+2f_4)+(f_6{J}_c+3f_5{J}_c+2f_4{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$

  $$T_4=\sqrt{{\displaystyle \frac{2\left(R_{03}+\frac{(1-⌀o-e^{-⌀o})f_3{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_3{J}_e}{2}\right)}{\left(\frac{d_3+⌀f_3+f_3{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}\right)}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $$T_6=\sqrt{{\displaystyle \frac{3\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{\theta }^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
• Step 4—Convert the constraints of $T_2-T_1\ge 0$ and $T_3-T_2\ge 0$ into equalities of $T_2-T_1=0$ and $T_3-T_2=0$, respectively, and next optimize $\mathrm{TC}\left(L\right)$ using the Lagrangian method defined as
  $$N(T_1,T_2,T_3,T_4,T_5,T_6,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4),$$

  where ${\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4$ are the Lagrange multipliers, attaining at
  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_6^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_6^2}}{\sigma }_{e}G\left(T_6{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_6^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_6^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_5^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_5^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_4^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_4^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{04}}{T_3^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_3^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{05}}{T_2^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_2^2}}\right)-{\tau }_4=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_6}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{06}}{T_1^2}}+{\displaystyle \frac{(d_6+⌀f_6+f_6{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_6{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}\right.\left.+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{6}h{J}_e}{2T_1^2}}\right)=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_3-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4).$$

  obtaining

  $$T_1=T_2=T_3=T_4=\sqrt{{\displaystyle \frac{3\left(\begin{array}{c}\hfill (R_{06}+3R_{05}+2R_{04}+2R_{03})+\frac{(1-⌀o-e^{-⌀o})(f_6{J}_c+3f_5{J}_c+2f_4{J}_c+2f_3{J}_c)({\sigma }_{e}G-h)}{{\theta }^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))(M_6+3M_5{J}_e+2M_4{J}_e+2M_3{J}_e)}{2}\end{array}\right)}{\left(\frac{d_6+3d_5+2d_4+2d_3⌀(f_6+3f_5+2f_4+2f_3)+(f_6{J}_c+3f_5{J}_c+2f_4{J}_c+2f_3{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  $$T_5=\sqrt{{\displaystyle \frac{\left(R_{02}+\frac{(1-⌀o-e^{-⌀o})f_2{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_2{J}_e}{2}\right)}{\left(\frac{d_2+⌀f_2+f_2{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}},$$

  and

  $$T_6=\sqrt{{\displaystyle \frac{\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_6>T_1$, which violates the condition of $0<T_1\le T_2\le T_3\le T_4\le T_5\le T_6$.
• Step 5—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, $T_5-T_4\ge 0$, and $T_6-T_5\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, $T_4-T_3=0$, $T_5-T_4=0$, and $T_6-T_5=0$, respectively, then optimize $\mathrm{TC}\left(L\right)$ using the Lagrangian method defined as
  $$N(T_1,T_2,T_3,T_4,T_5,T_6,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5)=\mathrm{TC}\left(L\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4)-{\tau }_5(T_6-T_5),$$

  where ${\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5$ are the Lagrange multipliers, to obtain
  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_6^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_6^2}}{\sigma }_{e}G\left(T_6{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_6^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_6^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_5^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_5^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_4^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_4^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{04}}{T_3^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_3^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{05}}{T_2^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_2^2}}\right)-{\tau }_4=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_6}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{06}}{T_1^2}}+{\displaystyle \frac{(d_6+⌀f_6+f_6{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_6{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{6}h{J}_e}{2T_1^2}}\right)-{\tau }_5=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4),{\displaystyle \frac{\partial N}{\partial {\tau }_5}}=-(T_6-T_5),$$

  generating

  $$T_1=T_2=T_3=T_4=T_5=\sqrt{{\displaystyle \frac{3\left(\begin{array}{c}\hfill (R_{06}+3R_{05}+2R_{04}+2R_{03}+3R_{02})+\frac{(1-⌀o-e^{-⌀o})(f_6{J}_c+3f_5{J}_c+2f_4{J}_c+2f_3{J}_c+3f_2{J}_c)({\sigma }_{e}G-h)}{{\theta }^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))(M_6+3M_5{J}_e+2M_4{J}_e+2M_3{J}_e+3M_2{J}_e)}{2}\end{array}\right)}{\left(\frac{d_6+3d_5+2d_4+2d_3+3d_2+⌀(f_6+3f_5+2f_4+2f_3+3f_2)+(f_6{J}_c+3f_5{J}_c+2f_4{J}_c+2f_3{J}_c+3f_2{J}_c)}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}$$

  and

  $$T_6=\sqrt{{\displaystyle \frac{\left(R_{01}+\frac{(1-⌀o-e^{-⌀o})f_1{J}_c({\sigma }_{e}G-h)}{{⌀}^2}-\frac{(o^2-p^2(1-\vartheta ))h{M}_1{J}_e}{2}\right)}{\left(\frac{d_1+⌀f_1+f_1{J}_c}{⌀}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  From the above expressions, we observe that $T_6>T_1$, which indicates the need for further adjustment of the constraints to ensure that the solution meets the required conditions.
• Step 6—Convert the constraints of $T_2-T_1\ge 0$, $T_3-T_2\ge 0$, $T_4-T_3\ge 0$, $T_5-T_4\ge 0$, and $T_6-T_5\ge 0$ into equalities of $T_2-T_1=0$, $T_3-T_2=0$, $T_4-T_3=0$, $T_5-T_4=0$, and $T_6-T_5=0$, respectively, and get the corresponding Lagrangian function expressed as
  $$N(T_1,T_2,T_3,T_4,T_5,T_6,{\tau }_1,{\tau }_2,{\tau }_3,{\tau }_4,{\tau }_5)=G\left(\mathrm{TC}\left(L\right)\right)-{\tau }_1(T_2-T_1)-{\tau }_2(T_3-T_2)-{\tau }_3(T_4-T_3)-{\tau }_4(T_5-T_4)-{\tau }_5(T_6-T_5).$$
  Hence, we attain

  $${\displaystyle \frac{\partial N}{\partial T_1}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{01}}{T_6^2}}+{\displaystyle \frac{(d_1+⌀f_1+f_1{J}_c)}{⌀T_6^2}}{\sigma }_{e}G\left(T_6{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_1{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_6^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{1}h{J}_e}{2T_6^2}}\right)+{\tau }_1=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_2}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{02}}{T_5^2}}+{\displaystyle \frac{(d_2+⌀f_2+f_2{J}_c)}{⌀T_5^2}}{\sigma }_{e}G\left(T_5{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_2{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_5^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{2}h{J}_e}{2T_5^2}}\right)-{\tau }_1+{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_3}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{03}}{T_4^2}}+{\displaystyle \frac{(d_3+⌀f_3+f_3{J}_c)}{⌀T_4^2}}{\sigma }_{e}G\left(T_4{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_3{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_4^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{3}h{J}_e}{2T_4^2}}\right)-{\tau }_2=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_4}}={\displaystyle \frac{2}{6}}\left(-{\displaystyle \frac{R_{04}}{T_3^2}}+{\displaystyle \frac{(d_4+⌀f_4+f_4{J}_c)}{⌀T_3^2}}{\sigma }_{e}G\left(T_3{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_4{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_3^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{4}h{J}_e}{2T_3^2}}\right)-{\tau }_3=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_5}}={\displaystyle \frac{3}{6}}\left(-{\displaystyle \frac{R_{05}}{T_2^2}}+{\displaystyle \frac{(d_5+⌀f_5+f_5{J}_c)}{⌀T_2^2}}{\sigma }_{e}G\left(T_2{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_5{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_2^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{5}h{J}_e}{2T_2^2}}\right)-{\tau }_4=0,$$

  $${\displaystyle \frac{\partial N}{\partial T_6}}={\displaystyle \frac{1}{6}}\left(-{\displaystyle \frac{R_{06}}{T_1^2}}+{\displaystyle \frac{(d_6+⌀f_6+f_6{J}_c)}{⌀T_1^2}}{\sigma }_{e}G\left(T_1{\displaystyle \frac{h{t}_1}{ht}}-t_1\right)-{\displaystyle \frac{f_6{J}_c(1-⌀o-e^{-⌀o})({\sigma }_{e}G-h)}{{⌀}^2{T}_1^2}}+{\displaystyle \frac{(o^2-(1-\vartheta )p^2)M_{6}h{J}_e}{2T_1^2}}\right)-{\tau }_5=0,$$

  and

  $${\displaystyle \frac{\partial N}{\partial {\tau }_1}}=-(T_2-T_1),{\displaystyle \frac{\partial N}{\partial {\tau }_2}}=-(T_3-T_2),{\displaystyle \frac{\partial N}{\partial {\tau }_3}}=-(T_4-T_3),{\displaystyle \frac{\partial N}{\partial {\tau }_4}}=-(T_5-T_4),{\displaystyle \frac{\partial N}{\partial {\tau }_5}}=-(T_6-T_5),$$

  obtaining $T_1=T_2=T_3=T_4=T_5=T_6$, where
  $$T_6=\sqrt{{\displaystyle \frac{3\left(\begin{array}{c}\hfill (R_{06}+3R_{05}+2R_{04}+2R_{03}+3R_{02}+R_{01})+\frac{(1-⌀o-e^{-⌀o})J_c(f_6+3f_5+2f_4+2f_3+3f_2+f_1)({\sigma }_{e}G-h)}{{⌀}^2}\\ \hfill -\frac{(o^2-p^2(1-\vartheta ))h{J}_e(M_6+3M_5+2M_4+2M_3+3M_2+M_1)}{2}\end{array}\right)}{\left({\displaystyle \frac{d_6+3d_5+2d_4+2d_3+3d_2+d_1+⌀(f_6+3f_5+2f_4+2f_3+3f_2+f_1)+J_c(f_6+3f_5+2f_4+2f_3+3f_2+f_1)}{⌀}}\right)\left(h⌀-\frac{h^2⌀}{{\sigma }_{e}G}\right)}}}.$$
  Thus, the solution of $T^{*}=(T_1,T_2,T_3,T_4,T_5,T_6)$ satisfies all the inequality constraints for the respective fuzzy number cases. By applying the nonlinear Lagrangian method, we calculate the TC for the model stated in (5). As fuzzy parameters are incorporated, the precision of the optimal solution improves. This calculation is systematically adapted for different fuzzy number types, ensuring that the solution respects the distinct subperiod divisions.

4. Results

In this section, we present the results of our study, demonstrating how advanced ML techniques were integrated with our novel inventory modeling strategy to evaluate various inventory management scenarios.

4.1. Machine Learning Methodology for Inventory Analysis

ML techniques offer powerful tools for analyzing and optimizing inventory. Weka, a comprehensive data mining platform, offers a wide range of algorithms that can be applied to practical inventory-related challenges. It supports data classification, clustering, preprocessing, regression, plotting, association rules, and other key functionalities.

In our analysis, we used a CSV file containing 12,291 rows with 18 parameters to calculate the values of TC and EOQ for both crisp and fuzzy inventory models. The TC was categorized as ‘profit’ and ‘non-profit’ based on the cycle length, where shorter cycle lengths increased profit and longer cycle lengths decreased profit.

To classify the dataset, LR was employed, which is a robust technique for classifying observations and identifying the most influential variables. It was chosen for its suitability in handling the probabilistic nature of our dataset, which encompasses various scenarios within the retailer supply chain. Recognized for its simplicity, efficiency, and effectiveness in managing large datasets with high variability and uncertainty, LR was ideal for our initial exploration.

While other classifiers such as support vector machines (SVMs), k-nearest neighbors (KNN), and random forest (RF) offer robust classification capabilities, their computational complexity depends on the data and tuning parameters. For this study, LR was chosen for its balance between predictive performance and computational efficiency. Its probabilistic and statistical outputs also provide straightforward interpretation, which is beneficial for decision makers.

Future studies could investigate the application of SVMs, KNN, RF, and other advanced techniques to our model, potentially offering deeper insights and improved predictive accuracy.

The methodological process adopted in our study is outlined in Figure 3 and includes the following steps:

• Step 1—Create a comprehensive dataset reflecting various scenarios within the retail supply chain, including demand data, production costs, and other parameters.
• Step 2—Fuzzify parameters, transforming deterministic values into fuzzy numbers to better represent uncertainties and variabilities in the supply chain.
• Step 3—Convert he dataset into attribute-relation file format (ARFF) format for compatibility with Weka.
• Step 4—Classify the dataset using the LR to categorize supply chain scenarios into ‘profitable’ and ‘non-profitable’.
• Step 5—Defuzzify the fuzzy results to crisp values for clearer interpretation.
• Step 6—Analyze the classification and defuzzification outcomes to evaluate model accuracy and applicability in real-world retail supply chain management.

Stratified cross-validation was employed to enhance the reliability of our findings, ensuring that each fold of the dataset accurately represented the overall distribution.

4.2. Classification Results for Trapezoidal Fuzzy Numbers Using Weka

Table 5 presents the results of classification using Weka and LR for TFNs. Table 6 reports the detailed accuracy by class, and Table 7 shows the confusion matrix.

From Table 5, note that the model correctly classified 11,924 instances (99.9748% of accuracy), and incorrectly classified 3 instances (0.0252% of error). The Kappa statistic of 0.9988 indicates almost perfect agreement between the predicted and observed classifications. The mean absolute error (MAE) is 0.0004, and the root mean squared error (RMSE) is 0.0125, indicating high precision and low error rates. Relative metrics such as the relative absolute error (0.2043%) and root relative squared error (3.9112%) further demonstrate the model’s accuracy.

Table 5. Results from LR for TFNs.

Metric	Value	Percentage
Correctly classified instances	11,924	99.9748%
Incorrectly classified instances	3	0.0252%
Kappa statistic	0.9988
Mean absolute error	0.0004
Root mean squared error	0.0125
Relative absolute error	0.2043%
Root relative squared error	3.9112%
Total number of instances	11,927

Table 5. Results from LR for TFNs.

Metric	Value	Percentage
Correctly classified instances	11,924	99.9748%
Incorrectly classified instances	3	0.0252%
Kappa statistic	0.9988
Mean absolute error	0.0004
Root mean squared error	0.0125
Relative absolute error	0.2043%
Root relative squared error	3.9112%
Total number of instances	11,927

Table 6 provides a detailed breakdown of the accuracy by class, showing that the TP rate is 1.000 for both ‘non-profit’ and ‘profit’ classes, which means that the model correctly classified every instance in each class. Similarly, the FP rate of 0.000 indicates that there were no misclassifications between the classes, that is, none of the ‘non-profit’ instances were mistakenly classified as ‘profit’ and vice versa. From Table 6, note that precision, recall, and F-measure scores are all exceptionally high (1.000 for ‘profit’ and 0.999 for ‘non-profit’), indicating that the model not only accurately classified the instances but that it also did so with minimal variance. The F-measure, a harmonic mean of precision and recall, confirms the model’s robustness in handling both classes without bias. Also from Table 6, observe that the Matthews correlation coefficient (MCC) of 0.999 further supports this, showing near-perfect correlation between predicted and observed classifications. The MCC is particularly helpful in situations with imbalanced datasets, as it takes into account TP, FP, and FN values. Here, it demonstrates the model’s ability to correctly classify instances in both classes equally well. The receiver operating characteristic (ROC) area is 1.000, indicating perfect sensitivity and specificity, meaning the model is able to distinguish between the ‘profit’ and ‘non-profit’ classes without error, even as classification thresholds vary. Similarly, the precision–recall curve (PRC) area is 1, showing that the model maintains high precision and recall, which is particularly helpful when there is class imbalance. These metrics underscore the model’s exceptional ability to balance TP and FP rates across different classification thresholds, confirming its robustness in a variety of scenarios.

Table 6. Detailed accuracy by class.

Class	TP Rate	FP Rate	Precision	Recall	F-Measure	MCC	ROC Area	PRC Area
Non-profit	1.000	0.000	0.998	1.000	0.999	0.999	1.000	1.000
Profit	1.000	0.000	1.000	1.000	1.000	0.999	1.000	1.000
Weighted average	1.000	0.000	1.000	1.000	1.000	0.999	1.000	1.000

Table 6. Detailed accuracy by class.

Class	TP Rate	FP Rate	Precision	Recall	F-Measure	MCC	ROC Area	PRC Area
Non-profit	1.000	0.000	0.998	1.000	0.999	0.999	1.000	1.000
Profit	1.000	0.000	1.000	1.000	1.000	0.999	1.000	1.000
Weighted average	1.000	0.000	1.000	1.000	1.000	0.999	1.000	1.000

The confusion matrix presented in Table 7 further illustrates the model’s performance. Out of 1378 instances classified as ‘non-profit’, all were correctly classified, with no misclassifications. Similarly, out of 10,549 instances classified as ‘profit’, 10,546 were correctly classified, with only 3 instances misclassified as ‘non-profit’. This strong performance indicates the LR model’s effectiveness in distinguishing between ‘profit’ and ‘non-profit’.

Table 7. Confusion matrix.

Classified as	Non-Profit (a)	Profit (b)
Non-Profit (a)	1378	0
Profit (b)	3	10,546

Table 7. Confusion matrix.

Classified as	Non-Profit (a)	Profit (b)
Non-Profit (a)	1378	0
Profit (b)	3	10,546

By analyzing the results reported in Table 5, Table 6 and Table 7, we can conclude that the LR model with TFNs performs exceptionally well in distinguishing between ‘profit’ and ‘non-profit’ classes. The high accuracy, low error rates, and strong performance metrics across both classes highlight the model’s effectiveness for this classification task.

4.3. Classification Results for Pentagonal Fuzzy Numbers Using Weka

Table 8 presents the results of classification using Weka and LR for PFNs. Table 9 reports the detailed accuracy by class, and Table 10 shows the confusion matrix.

Table 8 reports the model correctly classified 11,759 instances, corresponding to an accuracy of 98.5914%, while incorrectly classifying 168 instances, leading to an error rate of 1.4086%. The Kappa statistic of 0.9294 indicates strong agreement between the predicted and observed classifications, although it is less than the Kappa observed for TFNs. This may suggest that the increased complexity introduced by PFNs impacts the model’s ability to classify with the same level of precision as with TFNs.

Also from Table 8, note that the MAE of 0.0359 and the RMSE of 0.1202 remain acceptable, but the relative absolute error (17.5687%) and root relative squared error (37.5965%) are notably greater than for TFNs, indicating high variability and reduced precision in the model’s predictions when dealing with PFNs.

Table 8. Results from LR for PFNs.

Metric	Value	Percentage
Correctly classified instances	11,759	98.5914%
Incorrectly classified instances	168	1.4086%
Kappa statistic	0.9294
Mean absolute error	0.0359
Root mean squared error	0.1202
Relative absolute error	17.5687%
Root relative squared error	37.5965%
Total number of instances	11,927

Table 8. Results from LR for PFNs.

Metric	Value	Percentage
Correctly classified instances	11,759	98.5914%
Incorrectly classified instances	168	1.4086%
Kappa statistic	0.9294
Mean absolute error	0.0359
Root mean squared error	0.1202
Relative absolute error	17.5687%
Root relative squared error	37.5965%
Total number of instances	11,927

Table 9 details the accuracy by class, highlighting TP rates of 0.911 for the ‘non-profit’ class and 0.996 for the ‘profit’ class. The lower TP rate for the ‘non-profit’ class suggests that the model struggles more to correctly classify instances in this class, which is further evidenced by the FP rate of 0.089 for ‘profit’. Precision, recall, and F-measure remain high for both classes, reflecting a generally robust performance, although improvements could be made for the ‘non-profit’ class.

Table 9. Detailed accuracy by class.

Class	TP Rate	FP Rate	Precision	Recall	F-Measure	MCC	ROC Area	PRC Area
Non-profit	0.911	0.004	0.965	0.911	0.937	0.930	0.987	0.962
Profit	0.996	0.089	0.989	0.996	0.992	0.930	0.987	0.997
Weighted average	0.986	0.079	0.986	0.986	0.986	0.930	0.987	0.993

Table 9. Detailed accuracy by class.

Class	TP Rate	FP Rate	Precision	Recall	F-Measure	MCC	ROC Area	PRC Area
Non-profit	0.911	0.004	0.965	0.911	0.937	0.930	0.987	0.962
Profit	0.996	0.089	0.989	0.996	0.992	0.930	0.987	0.997
Weighted average	0.986	0.079	0.986	0.986	0.986	0.930	0.987	0.993

The confusion matrix presented in Table 10 further illustrates the classification performance. Out of 1378 instances classified as ‘non-profit’, 1256 were correctly classified, while 122 were misclassified as ‘profit’. In contrast, out of 10,549 instances classified as ‘profit’, 10,503 were correctly classified, with 46 misclassifications as ‘non-profit’. This difference in misclassification rates between the two classes should be addressed in future research, as the misclassification of ‘non-profit’ instances may have practical implications for inventory management decisions.

Table 10. Confusion matrix.

Classified as	Non-Profit (a)	Profit (b)
Non-Profit (a)	1256	122
Profit (b)	46	10,503

Table 10. Confusion matrix.

Classified as	Non-Profit (a)	Profit (b)
Non-Profit (a)	1256	122
Profit (b)	46	10,503

In summary, from Table 8, Table 9 and Table 10, the performance metrics suggest that while the LR model is effective for PFNs, it does face challenges with the ‘non-profit’ class, which may warrant further refinement of the model or exploration of alternative ML techniques to enhance classification accuracy for more complex fuzzy numbers. By analyzing these results, we can conclude that the LR model with PFNs performs adequately in distinguishing between ‘profit’ and ‘non-profit’ classes. While the model demonstrates high accuracy and strong performance metrics for the ‘profit’ class, the lower TP rate and higher misclassification rate for the ‘non-profit’ class suggest that there is room for improvement in handling more complex scenarios. Overall, the LR model remains effective for this classification task, but further refinement may enhance its performance in future studies.

4.4. Classification Results for Hexagonal Fuzzy Numbers Using Weka

Table 11 presents the results of classification using Weka and LR for HFNs. Table 12 reports the detailed accuracy by class, and Table 13 shows the confusion matrix.

From Table 11, observe that the model correctly classified 11,608 instances, which corresponds to 97.3254% accuracy, and incorrectly classified 319 instances, corresponding to 2.6746% error. The Kappa statistic of 0.8739 indicates substantial agreement between the predicted and observed classifications. The MAE is 0.0946, and the RMSE is 0.1817, indicating good precision and acceptable error rates. Relative metrics such as the relative absolute error (46.2914%) and root relative squared error (56.8264%) further demonstrate the model’s accuracy.

Table 11. Results from LR for HFNs.

Metric	Value	Percentage
Correctly classified Instances	11,608	97.3254%
Incorrectly classified Instances	319	2.6746%
Kappa statistic	0.8739
Mean absolute error	0.0946
Root mean squared error	0.1817
Relative absolute error	46.2914%
Root relative squared error	56.8264%
Total number of instances	11,927

Table 11. Results from LR for HFNs.

Metric	Value	Percentage
Correctly classified Instances	11,608	97.3254%
Incorrectly classified Instances	319	2.6746%
Kappa statistic	0.8739
Mean absolute error	0.0946
Root mean squared error	0.1817
Relative absolute error	46.2914%
Root relative squared error	56.8264%
Total number of instances	11,927

Table 12 details the accuracy by class, showing TP rates of 0.928 for the ‘non-profit’ class and 0.979 for the ‘profit’ class, with corresponding FP rates of 0.021 and 0.072, respectively. Precision, recall, and F-measure are also high, demonstrating the model’s robust performance across different classes. The weighted averages reflect the overall performance metrics, which are strong and consistent.

Table 12. Detailed accuracy by class.

Class	TP Rate	FP Rate	Precision	Recall	F-Measure	MCC	ROC Area	PRC Area
Non-profit	0.928	0.021	0.853	0.928	0.889	0.875	0.978	0.825
Profit	0.979	0.072	0.991	0.979	0.985	0.875	0.978	0.997
Weighted average	0.973	0.066	0.975	0.973	0.974	0.875	0.978	0.977

Table 12. Detailed accuracy by class.

Class	TP Rate	FP Rate	Precision	Recall	F-Measure	MCC	ROC Area	PRC Area
Non-profit	0.928	0.021	0.853	0.928	0.889	0.875	0.978	0.825
Profit	0.979	0.072	0.991	0.979	0.985	0.875	0.978	0.997
Weighted average	0.973	0.066	0.975	0.973	0.974	0.875	0.978	0.977

The confusion matrix presented in Table 13 further illustrates the model’s performance, showing that out of 1378 instances classified as ‘non-profit’, 1279 were correctly classified and 99 were misclassified as ‘profit’. Similarly, out of 10,549 instances classified as ‘profit’, 10,329 were correctly classified, with 220 instances misclassified as ‘non-profit’.

Table 13. Confusion matrix.

Classified as	Non-Profit (a)	Profit (b)
Non-Profit (a)	1279	99
Profit (b)	220	10,329

Table 13. Confusion matrix.

Classified as	Non-Profit (a)	Profit (b)
Non-Profit (a)	1279	99
Profit (b)	220	10,329

In summary, from Table 11, Table 12 and Table 13, we can conclude that the LR model with HFNs performs well in distinguishing between ‘profit’ and ‘non-profit’ classes. However, the slightly lower accuracy, higher misclassification rates, and decreased performance metrics for the ‘non-profit’ class, as reflected by the confusion matrix and class-specific metrics, indicate that the model faces limitations in correctly classifying certain instances. Despite these limitations, the model remains effective overall, but further refinements could improve its ability to handle the variability inherent in the ‘non-profit’ class.

4.5. Results of Cross-Validation Techniques

Next, we present the results of the classification task using three types of fuzzy numbers, namely TFNs, PFNs, and HFNs. For each case, cross-validation techniques were applied, and the performance of the LR classifier was evaluated. Figure 4 shows the box plots of the resupply order cost for the ‘profit’ and ‘non-profit’ categories, comparing the three fuzzy number types (TFNs, PFNs, and HFNs). The plot indicates that resupply order costs are much higher for ‘profit’ instances across all fuzzy number types. In particular, PFNs exhibit the highest median resupply order cost, at 2061, compared to 689 for ‘non-profit’ instances. Moreover, TFNs display more variability within the ‘profit’ category, with an interquartile range (IQR) of 1006, suggesting greater cost dispersion.

Figure 5 illustrates the box plots of manufacturing costs by ‘profit’ and ‘non-profit’ categories. Production costs are notably higher for ‘profit’ instances, with PFNs showing the highest median manufacturing cost of 517, while TFNs and HFNs have comparable medians around 507 and 499, respectively, in the ‘profit’ category. For ‘non-profit’ instances, the range of costs is lower and more consistent across all fuzzy number types, with a minimum value of 50 for each type.

Figure 6 shows the production cost by ‘profit’ and ‘non-profit’ category. The results reveal that ‘profit’ instances have substantially higher production costs across all fuzzy number types. PFNs show the highest median production cost of 411, compared to 63 for ‘non-profit’ instances. Additionally, the variability in production costs is more pronounced in ‘profit’ instances, with TFNs showing an IQR of 273, while PFNs and HFNs have an IQR of 272, indicating wider cost dispersion in profitable scenarios.

Overall, Figure 4, Figure 5 and Figure 6 demonstrate clear differences in resupply order, manufacturing, and production costs between ‘profit’ and ‘non-profit’ categories. The comparison between fuzzy number types (TFNs, PFNs, and HFNs) highlights how each method captures uncertainty and affects cost variability. PFNs generally exhibit the highest costs across all categories, while TFNs show the most variability in profit-making scenarios.

The second part of this analysis presents the results of a classification task using three types of cross-validation, namely TFNs, PFNs, and HFNs, along with their respective performance metrics. Trapezoidal cross-validation yields the highest accuracy and the lowest error rates compared to pentagonal and hexagonal cross-validations. Specifically, trapezoidal cross-validation achieved an accuracy of 99.97%, with an MAE of 0.0004 and an RMSE of 0.0125. These metrics demonstrate that the model predictions are highly accurate and precise. Moreover, the detailed accuracy by class shows perfect performance in classifying both ‘non-profit’ and ‘profit’ instances, with high precision, recall, and F-measure values for both classes. Pentagonal and hexagonal cross-validations exhibit slightly lower accuracies and higher error rates compared to trapezoidal cross-validation. Pentagonal cross-validation achieved an accuracy of 98.59%, with an MAE of 0.0359 and an RMSE of 0.1202, while hexagonal cross-validation yielded an accuracy of 97.33%, with an MAE of 0.0946 and an RMSE of 0.1817. Although these performance metrics indicate good results, they are not as strong as those of trapezoidal cross-validation. The superior performance of trapezoidal cross-validation can be attributed to its ability to delineate decision boundaries between classes more precisely, possibly due to the specific characteristics of TFNs, which may better capture the underlying patterns in the data than PFNs and HFNs. Consequently, trapezoidal cross-validation outperforms pentagonal and hexagonal cross-validations in terms of accuracy and error rates, making it the preferred choice for this classification task. The near-perfect classification results highlight its effectiveness in distinguishing between ‘non-profit’ and ‘profit’ instances, making it an optimal choice for similar tasks.

In summary, while the LR model performed well for all fuzzy number types, the TFN method achieved the highest performance, likely due to its ability to represent uncertainties in a more controlled and precise manner. The comparison between TFNs, PFNs, and HFNs underscores how different fuzzification approaches affect model performance, with the trapezoidal approach proving the most advantageous in this classification task.

5. Discussion and Conclusions

This study examined a manufacturing inventory model for deteriorating products, focusing on price-dependent demand and production rates associated with reliability under partial trade credits. The analysis highlighted the relevance of product reliability and trade credit policies in optimizing total production costs in competitive industrial environments.

The results of numerical simulations revealed that when the credit period extended to customers surpasses both the business cycle and the supplier’s credit term, the system’s total cost reaches its minimum. This is particularly relevant for industries dealing with perishable or short-lived products, such as pharmaceuticals, food products, and high-tech electronics, where efficient inventory management is critical for maintaining product quality and minimizing costs.

The proposed model assumes constant values for the deterioration rate and the production reliability factor. Although this assumption is common in inventory models, it may not fully capture the fluctuations observed in dynamic industrial settings. However, it provides a practical solution for scenarios where these rate and factor exhibit moderate variability. Additionally, the model focused on a single stock-keeping unit, which is appropriate for industries where high turnover or deterioration-sensitive products are central. Extending this model to accommodate multiple stock-keeping units would introduce substantial complexity and was intentionally excluded to maintain the clarity of the current analysis.

In terms of trade credit policies, a simplified approach was used, assuming a uniform credit period throughout the production cycle. While this reflects common industrial practices, future studies could explore more complex credit terms to account for variability in specific contexts. The assumptions and simplifications of the model make it well-suited for a broad range of industrial scenarios, particularly those prioritizing efficiency, predictability, and cost control in systems with deteriorating products and price-sensitive demand. These scenarios were rigorously addressed, and the proposed extensions aim to broaden the model’s applicability in more complex and dynamic environments.

The model is particularly relevant for industries dealing with perishable or deteriorating goods, such as pharmaceuticals and food production, where managing spoilage and maintaining product quality are critical. Minimizing production costs while balancing demand and ensuring reliability is essential in these sectors. The model is also applicable to industries that utilize partial trade credit policies, such as retail and consumer electronics, where fluctuating demand and cost control present important challenges. The inclusion of fuzzy logic to manage uncertainties in demand, reliability, and credit policies enhances the model’s practical applicability in such scenarios.

Future research could extend the model to handle multiple stock-keeping units, necessitating adjustments to the cost function to account for interactions between products. Managing multiple products introduces added complexity, as each stock-keeping unit may have distinct deterioration rates, demand patterns, and reliability factors. Optimizing both individual and interdependent costs, as well as shared resources, would require the development of sophisticated algorithms. This extension would enhance the model’s relevance for industries with diverse product lines, such as retail chains and manufacturing sectors.

Additionally, exploring more advanced fuzzy numbers, such as the spherical case [52,53], could better model uncertainties in complex scenarios. While trapezoidal and pentagonal fuzzy numbers balance simplicity and computational efficiency in the current model, more intricate uncertainty patterns might benefit from advanced techniques. Incorporating varying deterioration rates or fluctuating production reliability would further enhance the model’s realism, making it more applicable to industries with varying environmental conditions or production capabilities. Moreover, integrating advanced machine learning techniques, such as deep neural networks or ensemble models [54], could improve the model’s predictive accuracy and its ability to handle highly uncertain inventory scenarios. In some cases, quantile regression could serve as an alternative to median-based methods to better capture the distribution of key variables [55].

Exploring reliability models beyond the inventory domain, such as epidemic models or hybrid systems, could also broaden the model’s applicability, providing deeper insights into the role of reliability in diverse contexts.

In conclusion, this study presented a novel approach to managing inventory systems with deteriorating goods under partial trade credit policies. By incorporating fuzzy numbers to address uncertainties in production costs, reliability, and demand, the model demonstrated high accuracy in decision making under uncertainty. The primary contributions of this research lie in its integration of fuzzy logic for uncertainty management and its applicability to real-world scenarios, particularly in industries dealing with perishable goods and fluctuating demand.