A numerical method for parabolic complementarity problem

Haiyan Song; Fei Sun; Haiyan Song; Fei Sun

doi:10.3934/era.2023052

Electronic Research Archive

2023, Volume 31, Issue 2: 1048-1064. doi: 10.3934/era.2023052

Previous Article Next Article

Theory article Special Issues

A numerical method for parabolic complementarity problem

Haiyan Song ¹,
Fei Sun ^{2
,
,}

1.
School of Computer and Data Engineering, NingboTech University, Ningbo 315100, China
2.
School of Mathematics and Computational Science, Wuyi University, Jiangmen 529020, China

Received: 06 October 2022 Revised: 24 November 2022 Accepted: 30 November 2022 Published: 13 December 2022

In this paper, we study the numerical solution of a parabolic complementarity problem which is a widely used model in many fields, such as option pricing, risk measures, etc. Using a power penalty method we represent the complementarity problem as a nonlinear parabolic partial differential equation (PDE). Then, we use the trapezoidal rule as the time discretization, for which we have to solve a nonlinear equation at each time step. We solve such a nonlinear equation by the fixed-point iteration and in this methodology solving a tridiagonal linear system is the major computation. We present an efficient backward substitution algorithm to handle this linear system. Numerical results are given to illustrate the advantage of the proposed algorithm (compared to the built-in command backslash in Matlab) in terms of CPU time.

Keywords:

Citation: Haiyan Song, Fei Sun. A numerical method for parabolic complementarity problem[J]. Electronic Research Archive, 2023, 31(2): 1048-1064. doi: 10.3934/era.2023052

Related Papers:

[1]	Massoumeh Nazari, Mahmoud Dehghan Nayeri, Kiamars Fathi Hafshjani . Correction: Developing mathematical models and intelligent sustainable supply chains by uncertain parameters and algorithms. AIMS Mathematics, 2024, 9(9): 25223-25231. doi: 10.3934/math.20241230
[2]	Zhimin Liu, Ripeng Huang, Songtao Shao . Data-driven two-stage fuzzy random mixed integer optimization model for facility location problems under uncertain environment. AIMS Mathematics, 2022, 7(7): 13292-13312. doi: 10.3934/math.2022734
[3]	Zhimin Liu . Data-driven two-stage sparse distributionally robust risk optimization model for location allocation problems under uncertain environment. AIMS Mathematics, 2023, 8(2): 2910-2939. doi: 10.3934/math.2023152
[4]	Shuai Huang, Youwu Lin, Jing Zhang, Pei Wang . Chance-constrained approach for decentralized supply chain network under uncertain cost. AIMS Mathematics, 2023, 8(5): 12217-12238. doi: 10.3934/math.2023616
[5]	Yuejia Dang . Intelligent optimization algorithm for strategic planning in economics with multi-factors assessment: A MEREC-driven Heronian mean framework. AIMS Mathematics, 2025, 10(5): 10866-10897. doi: 10.3934/math.2025494
[6]	Ibrahim M. Hezam, Pratibha Rani, Arunodaya Raj Mishra, Ahmad Alshamrani . An intuitionistic fuzzy entropy-based gained and lost dominance score decision-making method to select and assess sustainable supplier selection. AIMS Mathematics, 2023, 8(5): 12009-12039. doi: 10.3934/math.2023606
[7]	Xiaodie Lv, Yi Liu, Yihua Zhong . A novel method to solve the optimization problem of uncertain network system based on uncertainty theory. AIMS Mathematics, 2023, 8(3): 5445-5461. doi: 10.3934/math.2023274
[8]	Jun Wang, Dan Wang, Yuan Yuan . Research on low-carbon closed-loop supply chain strategy based on differential games-dynamic optimization analysis of new and remanufactured products. AIMS Mathematics, 2024, 9(11): 32076-32101. doi: 10.3934/math.20241540
[9]	Muhammad Ihsan, Muhammad Saeed, Atiqe Ur Rahman, Mazin Abed Mohammed, Karrar Hameed Abdulkaree, Abed Saif Alghawli, Mohammed AA Al-qaness . An innovative decision-making framework for supplier selection based on a hybrid interval-valued neutrosophic soft expert set. AIMS Mathematics, 2023, 8(9): 22127-22161. doi: 10.3934/math.20231128
[10]	Muhammad Saqlain, Xiao Long Xin, Rana Muhammad Zulqarnain, Imran Siddique, Sameh Askar, Ahmad M. Alshamrani . Energy supplier selection using Einstein aggregation operators in an interval-valued q-rung orthopair fuzzy hypersoft structure. AIMS Mathematics, 2024, 9(11): 31317-31365. doi: 10.3934/math.20241510

Abstract

1. Introduction

Sustainable supply chain management involves integrating economic, social, and environmental considerations into intra-organizational systems developed for managing materials and costs of production, distribution, and procurement of services/products ^[1]. Limited studies have addressed environmental issues in supply chain management, mainly focusing on reverse logistics, closed-loop, and green supply chain management. This may be due to the complexity and ambiguity of sustainability research, the blurred boundaries between interdisciplinary fields, and the vague definitions provided for sustainable development and sustainability. Social sustainability entails meeting people's needs, economic sustainability requires maximizing profits, and environmental sustainability calls for minimizing waste and non-renewable material consumption while avoiding long-term environmental harm ^[2]. Deteriorating environmental conditions and associated uncertainties make adopting green production processes in supply chain management crucial ^[3].

Artificial Intelligence (AI) is an effective approach to handling uncertainty ^[4]. Various environmental fluctuations and operational challenges in institutions complicate the collection of planning information. Moreover, the complexity of market systems makes predictions challenging. However, forecasting market price and demand is highly beneficial for investors.

Supply chain problems in the automotive industry can significantly negatively impact the environment and annual costs. These issues can result in disruptions in production and delays in the delivery of components and materials, leading to increased energy consumption, waste generation, and greenhouse gas emissions. They can also cause higher expenses. Furthermore, such problems can damage the reputation of automotive companies, result in customer dissatisfaction, reduce customer loyalty, and hinder progress towards sustainability goals, such as reducing pollutants and utilizing renewable materials. To address these challenges and improve supply chain performance, it is essential to implement efficient and sustainable practices.

Using traditional plastic materials in automotive manufacturing can significantly impact environmental pollution, such as extracting fossil fuels that contribute to greenhouse gas emissions. Additionally, plastic waste can take hundreds of years to decompose, leading to pollution in landfills and of marine environments. Environmental pollution can occur through various activities such as emissions, waste disposal, resource consumption, and habitat destruction. Companies can contribute to air pollution by releasing harmful gases and particles into the atmosphere. These can be from industrial processes, burning fossil fuels, or improper waste management practices. Air pollution can have serious health effects on both humans and wildlife. Companies often discharge untreated or improperly treated wastewater into rivers, lakes, or oceans. This can contaminate water bodies with toxic chemicals, heavy metals, and other pollutants, affecting aquatic life and posing risks to human health. Improper disposal or mishandling of hazardous waste can contaminate soil, making it unsuitable for agriculture or other purposes. Chemical spills, improper storage and disposal of industrial waste, and landfills are common sources of soil pollution caused by companies. Industrial activities, machinery, and transportation associated with companies can generate excessive noise levels, causing disturbances to both humans and wildlife. Prolonged exposure to high noise levels can lead to stress, hearing problems, and other health issues. Companies that emit greenhouse gases, such as carbon dioxide (CO₂), methane (CH₄), and nitrous oxide (N₂O) contribute to climate change. Industries like power generation, manufacturing, and transportation are principal sources of greenhouse gas emissions, which can lead to global warming and its associated impacts.

The purpose of a sustainable automotive supply chains with an AI approach in conditions of uncertainty is to address challenges and mitigate risks associated with operating in unpredictable environments. By leveraging AI, the supply chain can better adapt to changing circumstances and make informed decisions in real-time. Most models addressing uncertainty in sustainable supply chains consider uncertain parameters as a singular type of uncertainty, often posing challenges in real-world scenarios. This model can increase reliability, remain stable against potential events, consider various uncertainties, and accurately find the optimal solution. This perspective on uncertainty in mathematical models brings them closer to reality, produces more reliable solutions, reduces computation time and problem complexity, and enhances the robustness and resilience of sustainable supply chains against risks and uncertainties. In this study, a robust heuristic-mathematical programming model is used to make the obtained solutions reliable because of their robustness. Accordingly, dealing with the uncertainty in the sustainable supply chain and obtaining results with sufficient certainty require robust programming under uncertain conditions (e.g., random, fuzzy) using AI to provide reliable results to managers. The present study is of particular importance due to its innovation, diversity, and dynamics in the face of the uncertainties affecting business, especially in active supply chains in the automotive industry. In addition to the unprecedented research and innovative aspects of AI capabilities in predicting uncertain parameters of sustainable supply chains, this study developed a sustainable supply chain model by employing robust heuristic optimization techniques under uncertain conditions. This model is helpful in rapidly foreseeing uncertainty.

Also, this study, by utilizing hybrid model AI, the mathematical optimization algorithms NSGA-Ⅱ and HBA can improve the quality of decisions while reducing costs in a sustainable supply chain. These algorithms possess exceptional efficacy in minimizing environmental pollutants and system expenses, ultimately leading to improved decision-making outcomes. It also attempts to answer the following research questions:

● What are the parameters causing uncertainty in the mathematical model of the sustainable supply chain?

● What are the robust parameters in the sustainable supply chain model?

● Under uncertain conditions, which parameters sustainable supply chain model can be predicted using the AI approach?

2. Literature review

Generally, this literature review covers the sustainable supply chain (SSC), emphasizing Artificial Intelligence (AI), robust mathematical models, NSGA-Ⅱ, HBA algorithms, uncertainty, and robust optimization.

After conducting a thorough review of both local and international research, it has been determined that, despite the numerous studies conducted on sustainable supply chains, there are shortcomings that require improvement. In a recent study, researchers used intelligent methods to improve the efficiency and sustainability of the supply chain ^[5]. In another study, Wang et al. utilized intelligentization in the supply chain to enhance environmental sustainability. They assessed the ecological repercussions of the supply chain by using AI-based techniques and data analysis, and proposed measures to mitigate these repercussions. By utilizing AI and optimization techniques, they significantly enhanced energy consumption and reduced greenhouse gas emissions throughout the supply chain ^[6]. Vali-Siar, and Roghanian ^[7], developed a stable, flexible, and responsive mixed SC network under uncertainty due to the disruption of the COVID-19 pandemic. The researchers developed a multi-objective mixed integer linear programming (MOMILP) model for solving the flexible, responsive, and stable open-loop and closed-loop SC network design problem. The authors also developed a Lagrangian relaxation method and a constructive heuristic algorithm to overcome the complexity of the problem and solve large-scale problems. Hashim et al. ^[8] investigated a new meta-heuristic algorithm for solving optimization problems called the Honey Badger Algorithm, inspired by intelligent search behavior to develop an efficient mathematical strategy for solving optimization problems. They compared several well-known meta-heuristic algorithms with the solutions obtained from the HBA algorithm, including simulated annealing (SA), particle swarm optimization (PSO), covariance matrix adaptive evolution strategy (CMA-ES), and adaptive differential of evolution with linear population size reduction (L-SHADE), Moth-flame (MFO), Elephant Herding (EHO), the Whale (WOA), the Grasshopper (GOA), Thermal Exchange (TEO), and the Harris optimization (HHO). Yang et al. ^[9] have investigated the impact of manufacturing intelligence on green innovation performance in China. Their study has examined the effect of manufacturing intelligence on the effectiveness of green innovation and its internal mechanisms from a theoretical and empirical perspective. From the view of static and dynamic productivity, a dynamic spatial delay model, a mediating effect model, and a moderate effect model have been used for analysis. Homayouni et al. ^[10] investigated a robust optimization approach for carbon outflow control under uncertainty, which was solved using multiple-choice programming.

Based on the review of domestic and foreign research literature, it is evident that while many studies have been conducted on sustainable supply chains, nevertheless gaps need to be addressed. Foroozesh et al. ^[11] mainly focused on demand and fuzzy uncertainty while ignoring costs, carbon emissions, and economic dimensions. Hashim et al. ^[8] emphasized significant costs, demand and sustainable supply chains, and parameter prediction using neural networks. Instead, they introduced a meta-heuristic HBA algorithm. Yang et al. ^[9] mainly focused on demand and fuzzy uncertainty while ignoring carbon emissions and economic aspects. Homayouni et al. ^[10] considered pollutants, cost, demand, and economic dimension. But, similar to the above, they did not predict the neural network parameters and did not use HBA in their model. Govindan et al. ^[12] mainly focused on the demand and economic dimension and did not pay much attention to carbon emissions and uncertainty. Dormaz and Bilgen ^[13] only considered costs and did not pay much attention to carbon pollutants, demand, economic dimension, and uncertainty. Santibanes et al. ^[14] focused on the parameters related to carbon pollutants and costs but did not consider the economic dimension and the demand parameter. Instead, they used a stochastic uncertainty technique. In Franco and Alfonso Lizarazzo's ^[15] study, the aspect of uncertainty is not mentioned in any of the main parameters of the sustainable supply chain model. Gholizadeh et al. ^[16] did not consider carbon pollutants and only focused on demand, economic dimension, and costs.

Based on the investigated gaps, the present study, showing the capability of the sustainable SC model, considers demand, costs, and capacity reduction due to disruption (factory, supplier, recycling, and distributor) as non-deterministic parameters. It uses AI to predict shipping cost parameters. In addition, the NSGA-Ⅱ and HBA algorithms are used to optimize the model in the NP-hard problem to reduce the computation time. In addition to speeding up managers' problem-solving and decision-making processes, this model minimizes environmental risks and maximizes profits.

3. Methodology

The present study adopted a descriptive-analytical approach. A hybrid strategy was employed to tackle the inherent uncertainty of the problem. This strategy consisted of AI, the robust optimization model, and the meta-heuristic Honey Badger and NSGA-Ⅱ algorithms. This model produced greater resilience than basic models and was designed to manage problem complexity and solve large-scale instances. Additionally, this innovative heuristic model could significantly reduce computation time for large-scale problems. Data was analyzed using GAMS software with a simplex solver and Minitab 16, MATLAB 2015, and Microsoft Visio software. The non-deterministic parameters of shipping costs were predicted using neural networks and integrated into the model. However, due to the absence of disruption-induced capacity reduction rate and demand parameters, these parameters were entered non-deterministically and addressed through the robust model proposed by Mulvey. In general, the Honey Badger and NSGA-Ⅱ algorithms contribute to improving the performance and sustainability of supply chains by providing unique optimization methods. These algorithms, considering uncertain conditions, sudden changes, and multiple objectives, enable resource optimization, waste reduction, and improved supply chain performance. The combination of AI, HBA, and NSGA-Ⅱ and the mathematical robustness model provides a new approach. NSGA-Ⅱ helps balance conflicting goals. Utilizing AI and advanced optimization techniques can mitigate risks associated with uncertainty. Companies can achieve a more sustainable supply chain by using AI-based optimization algorithms like HBA and NSGA-Ⅱ, enabling informed decisions, considering uncertainties, and balancing economic, environmental, and social objectives. This approach improves efficiency, reduces environmental impact, and enhances overall sustainability performance. The aim is to provide accurate, robust, and adaptive models that effectively handle uncertainty, leading to improved decision-making and outcomes. This study presents a precise, adaptable, and resilient model using AI, the HBA, and NSGA-Ⅱ in the face of uncertainty. By utilizing mathematical models in the context of a sustainable supply chain under uncertainty, with the approach of AI and multi-objective algorithms desirable and fruitful results can be achieved. These methods, employing AI algorithms and optimization, allow for the examination and evaluation of various variables in the supply chain, leading to better and more optimal decisions in uncertain conditions. This study applied various innovations that maintain a specific priority:

1) Implementing a closed-loop sustainable chain model.

2) Using two objective functions in the robust model.

3) Including the disruption-induced capacity reduction parameter under uncertainty conditions.

4) Applying the ε-constraint technique predominantly used in cases where the values of objective functions differ and a priority exists.

5) Comparing the newly developed HBA algorithm (2022) with the older NSGA-Ⅱ algorithm.

6) Utilizing AI to predict non-deterministic parameters of shipping costs, ranging from the supply of polymer materials to the factory, from the factory to distributors of automotive plastic accessories, from distributors to customers of plastic car accessories, and from customers to the recycling and back to the factory.

7) Employing a hybrid approach consisting of AI, the robust model, the HBA, and NSGA-Ⅱ.

8) Precisely adjusting the parameters of meta-heuristic algorithms due to the inherently complex nature of design variables, which are largely non-linear and non-convex.

9) Highlighting the high speed of convergence in the HBA for NP-hard problems compared to other techniques.

Also, the conceptual model of the study is shown in Figure 1.

Figure 1. Conceptual models.

DownLoad: Full-Size Img PowerPoint

According to Azar and Farrokh ^[17], robust stochastic optimization, first presented by Mulvey ^[18], is a stochastic programming technique to manage the impacts of uncertain parameters within optimization models. This technique combines ideal programming formulations with scenario-based data related to the problem for robust modeling. Robust modeling is comprised of two decision-making variables: control and structural (design). For structural variables, decisions are made before the realization of random parameters, and once this realization is completed, these variables cannot be modified or adjusted. Control variables, on the other hand, are modified or adjusted when uncertainty parameters occur. The solution obtained with this model has lower sensitivity to variations in input data. In the current theory, a solution to the optimization problem is considered optimally (i.e., the robustness of the key) and justifiably (i.e., model robustness) robust. For a solution to be justifiably robust, it should maintain its justifiability within the problem's constraints for nearly all possible scenarios. The technique proposed by Mulvey et al. ^[18] uses a penalty function to measure the unjustifiability level of the solution to maintain its justifiability. The penalty function must assume a low value. Optimization robustness can be defined as a state where the objective function value for the solution retains its near-optimal value or its unfavorable deviations from the optimal value are minimal for every scenario. In Mulvey's model, the objective function has a penalty function for solving and the robustness of the model, each receiving weight per the modeler's defined preferences ^[17]. In Mulvey's model ^[18] (Eq (2.1)), the vector y represents control variables, and the vector x describes design variables.

A and B denote vectors of parameter coefficients, and b represents vectors of parameters (values on the right). A and B are definite values, but parameters B, C, and e are non-deterministic. The specific state of the parameters is considered a scenario, denoted by s. The probability of each scenario is determined by the parameter Ps. Ω represents a set of scenarios. As a result, for each scenario ( $s\in \mathit{\Omega}$ ), Bs and Cs are the coefficients that indicate uncertainty. After the scenario, y (the control variable) is adjusted. Due to the uncertainty of parameters, the model may be unjustified for some non-deterministic plans. Consequently, ηs can be referred to as the model when subjected to unreasonable methods. ηs is zero when the model is justified. Otherwise, it will assume a positive value. The objective function is given by the following:

$\begin{array}{l} \ \ \ \mathit{\boldsymbol{M}}\mathit{\boldsymbol{i}}\mathit{\boldsymbol{n}}\mathit{\boldsymbol{\sigma }}(\mathit{\boldsymbol{x}}, {\mathit{\boldsymbol{y}}}_{{\bf{1}}}, {\mathit{\boldsymbol{y}}}_{{\bf{2}}}, ....., {\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{s}}})+\mathit{\boldsymbol{\omega }}\mathit{\boldsymbol{\rho }}({\mathit{\boldsymbol{\eta }}}_{{\bf{1}}}, {\mathit{\boldsymbol{\eta }}}_{{\bf{2}}}, ......, {\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{s}}})\\ \ \ \ \mathit{\boldsymbol{s}}.\mathit{\boldsymbol{t}}.\\ \ \ \ \mathit{\boldsymbol{A}}\mathit{\boldsymbol{x}} = \mathit{\boldsymbol{b}}, \\ {\mathit{\boldsymbol{B}}}_{\mathit{\boldsymbol{s}}}\mathit{\boldsymbol{x}}+{\mathit{\boldsymbol{c}}}_{\mathit{\boldsymbol{s}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{s}}}+{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{s}}} = {\mathit{\boldsymbol{e}}}_{\mathit{\boldsymbol{s}}}, \mathit{\boldsymbol{x}}\ge 0, {\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{s}}}\ge 0, {\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{s}}}\ge 0, \forall \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{\mathit{ \boldsymbol{\varOmega}} }} . \end{array}$

(2.1)

The solution robustness is calculated by the first term, which indicates the level of risk aversion of decision-makers and their willingness to reduce costs. The second function measures the model's robustness by penalizing solution violations. It is combined using a coefficient (weight) γ to strike a balance between the model and solution robustness. The effect of γ can be understood further by inserting a small value into this parameter. It then minimizes the first objective function and increases the probability of achieving impossible solutions. For large γ values, the solution feasibility will increase, but the first part of objective function σ(x, y1, y2, ...) takes on higher values. ξ = f(x, y) and ξ are recognized as the utility and cost functions, respectively. High values of variance ξs = f(x, ys) for each scenario indicate risk decisions. In other words, minor variations in uncertain parameters may cause significant variations in the function f values. Mulvey et al. ^[19] confirmed the robustness of the solution by employing the terms in Eq (2.2). δ represents the weight indicating the variance of the solution.

$\sigma \left(0\right) = \sum\limits_{s\in \mathit{\Omega} }{p}_{s}{\xi }_{s}+\gamma {\sum\limits_{s\in \mathit{\Omega} }{p}_{s}\left({\xi }_{s}-\sum\limits_{s\in \mathit{\Omega} }{p}_{s}^{'}{\xi }_{s}^{'}\right)}^{2}.$

(2.2)

In Eq (2.2), the squared term ${{\sum }_{\mathit{\boldsymbol{s}}\in \mathit{\Omega} }{\mathit{\boldsymbol{p}}}_{\mathit{\boldsymbol{s}}}\left({\mathit{\boldsymbol{\xi }}}_{\mathit{\boldsymbol{s}}}-{\sum }_{\mathit{\boldsymbol{s}}\in \mathit{\Omega} }{\mathit{\boldsymbol{p}}}_{\mathit{\boldsymbol{s}}}^{'}{\mathit{\boldsymbol{\xi }}}_{\mathit{\boldsymbol{s}}}^{'}\right)}^{2}$ leads to increased computation time required to solve the model. Thus, the absolute value of the term was provided ^[19] to reduce the operations associated time with total computations:

$\sigma \left(o\right) = {\sum }_{s\in \mathit{\Omega} }{p}_{s}{\xi }_{s}+\gamma {\sum }_{s\in \mathit{\Omega} }{p}_{s}\left|{\xi }_{s}-{\sum }_{s\in \mathit{\Omega} }{p}_{s}^{'}{\xi }_{s}^{'}\right| .$

(2.3)

Two additional variables ${\mathcal{Q}}_{\mathrm{s}}^{+}$ and ${\mathcal{Q}}_{\mathrm{s}}^{-}$ are defined in Eq (2.3), which denotes absolute value and nonlinear function, so that the obtained objective function is linearized. If ${\sum }_{\mathrm{s}\in \mathrm{\mathit{\Omega} }}{\mathrm{p}}_{\mathrm{s}}^{'}{\mathrm{\xi }}_{\mathrm{s}}^{'}$ is greater than ${\mathrm{\xi }}_{\mathrm{s}}$ , then ${\mathcal{Q}}_{\mathrm{s}}^{-}$ is returned; otherwise, ${\mathcal{Q}}_{\mathrm{s}}^{+}$ is not returned. Consequently, Eq (2.3) can be reformulated as Eq (2.4):

$\begin{array}{l} Min\sum\limits_{s\in \mathit{\Omega} }{p}_{s}{\xi }_{s}+\gamma \sum\limits_{s\in \mathit{\Omega} }{p}_{s}({Q}_{s}^{+}+{Q}_{s}^{-})\\ s.t.\\ {\xi }_{s}-\sum\limits_{s\in \mathit{\Omega} }{p}_{s}^{'}{\xi }_{s}^{'} = {Q}_{s}^{+}+{Q}_{s}^{-}, s\in \mathit{\Omega} , \\ {Q}^{+}, {Q}^{-}\ge 0, s\in \mathit{\Omega} . \end{array}$

(2.4)

In line with the constraints of Eq (2.4), it should be noted that one of the value $Q_s^ -, Q_s^ +$ is always equal to zero for every $\gamma \geqslant 0$ . Based on Eq (2.1), the objective function for the robust model can be reformulated as Eq (2.5):

$\begin{array}{l} Min\sum\limits_{s\in \mathit{\Omega} }{p}_{s}{\xi }_{s}+\gamma \sum\limits_{s\in \mathit{\Omega} }{p}_{s}\left[\right({\xi }_{s}-\sum\nolimits_{s\in \mathit{\Omega} }{p}_{s}^{'}{\xi }_{s}^{'})+2{\theta }_{s}], \\ s.t. \\ \ \ \ \ \ \ {\xi }_{s}-\sum\limits_{s\in \mathit{\Omega} }{p}_{s}^{'}{\xi }_{s}^{'}+{\theta }_{s}\ge 0, {\theta }_{s}\ge 0 . \end{array}$

(2.5)

The expected value of the error vector is considered as a penalty function ${\sum }_{s\in \mathit{\Omega} }{p}_{s}{\eta }_{s}$ , and, finally, the robust programming model can be presented as Eq (2.6):

$\begin{array}{l} Min\ \ \ \mathop {}\limits_{} \mathop {}\limits_{} \sum\limits_{s \in \mathit{\Omega} } {{p_s}{\xi _s} + \gamma } \sum\limits_{s \in \mathit{\Omega} } {{p_s}\left[ {\left( {{\xi _s} - \sum\limits_{s \in \mathit{\Omega} } {{p_{s'}}{\xi _{s'}}} } \right) + 2{\theta _s}} \right]} + \omega \sum\limits_{s \in \mathit{\Omega} } {{p_s}{\eta _s}} \hfill \\ s.t.\ \ \ \mathop {}\limits_{} \mathop {}\limits_{} Ax = b, \hfill \\ \ \ \ \ \ \ \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} {B_s}x + {c_s}{y_s} + {\eta _s} = {e_s}, \mathop {}\limits_{} \mathop {}\limits_{} \forall s \in \mathit{\Omega} \hfill \\ \ \ \ \ \ \ \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} {\xi _s} - \sum\limits_{s \in \mathit{\Omega} } {{p_{s'}}{\xi _{s'}}} + {\theta _s} \geqslant 0\mathop {}\limits_{} s \in \mathit{\Omega} \hfill \\ \ \ \ \ \ \ \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} \mathop {}\limits_{} {\theta _s} \geqslant 0, x \geqslant 0, {y_s} \geqslant 0, {\eta _s} \geqslant 0, \forall s \in \mathit{\Omega} \hfill \end{array}.$

(2.6)

3.1. Mathematical model design

3.1.1. Robust indexes

Robust indexes are presented in Table 1.

Table 1. Robust indexes and their corresponding groups.

Index	Group
S, s^′∈S	Set of scenarios
i∈I	Suppliers of polymer material
j∈J	Factories
k∈K	Distributors of automotive plastic accessories
c∈C	Customers of plastic car accessories
v∈V	Vehicles
h∈H	Recycling center
t∈T	Period series
p∈P	Products of automotive plastic pieces
r∈R	Polymer materials

| Show Table

DownLoad: CSV

3.1.2. AI parameters

$\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{x}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : Cost of transportation of polymer materials r from the supplier of polymer materials i to the factory j within the period t using the vehicle v after determination based on the neural network (NN).

$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : Cost of shipping polymer materials r from the recycling center h to the factory j within the period t using the vehicle, v after determination based on the neural network (NN).

$\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : The shipping cost of returned products of automotive plastic pieces p from the customers of plastic car accessories c to the recycling center h within the period t using the vehicle v after determination based on the neural network (NN).

$\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}^{'}$ : The shipping cost of the products of automotive plastic pieces p from the distribution center k to the customers of plastic car accessories c within the period t using the vehicle v after determination based on the neural network (NN).

$\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : The shipping cost of products of automotive plastic pieces p from the factory center j to the distributors of automotive plastic accessories k within the period t using the vehicle v after determination based on the neural network (NN).

3.1.3. Conceptual model parameters in condition certain

${\mathit{\boldsymbol{x}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : Shipment volume of polymer materials r from the supplier of polymer materials i to the factory j during period t using the vehicle v.

$\mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : Shipment volume of the polymer materials r from recycling center h to the factory j during period t using the vehicle v.

${\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The returned amount of products of automotive plastic pieces p from customers of plastic car accessories c to recycling center h during period t using the vehicle v.

${\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}^{'}$ : Amount of the products of automotive plastic pieces p from distribution center k to customers of plastic car accessories c during period t using vehicle v.

${\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}}$ : Amount of the products of automotive plastic pieces p by factory center j to distributors of automotive plastic accessories k during period t using vehicle v.

3.1.4. Robust parameters

${\mathit{\boldsymbol{\lambda }}}_{\mathit{\boldsymbol{r}}\mathit{\boldsymbol{p}}}$ : The rate at which polymer material r is converted into the products of automotive plastic pieces p.

$\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The inventory cost of the products of automotive plastic pieces p at the factory j within the period t.

$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The inventory cost of the products of automotive plastic pieces p at the distributors of automotive plastic accessories k within the period t.

$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}}$ : The inventory cost of the polymer material r at the recycling center h within the period t.

${\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The non-deterministic demand for products of automotive plastic pieces p for the customers of plastic car accessories c in period t per scenario s.

$\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The factory capacity j for the products of automotive plastic pieces p within the period t.

$Y{J}_{j}$ : If a factory is selected, one; otherwise, zero.

$\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}}$ : The polymer material supplier i, for supplying polymer material r within the period t.

$Y{I}_{is}$ : If the supplier is selected, one; otherwise, zero.

$\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The capacity of the distributors of automotive plastic accessories k for the products of automotive plastic pieces p within the period t.

$Y{K}_{k}$ : If the distributor is selected, one; otherwise, zero.

$\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The capacity of the recycling center h for the products of automotive plastic pieces p within the period t.

$Y{H_h}$ : If the recycling centers are selected, one, otherwise zero.

$y{y}_{kcpt}$ : If the quantity of goods is shipped from the distribution center to the customer in the period, one; otherwise, zero.

${Z}_{kcpt}$ : If the number of goods from the distribution center to the customer in the period was lower than the minimum quality in terms of quality, one, otherwise zero.

${\mathit{\boldsymbol{\varphi }}}_{\mathit{\boldsymbol{r}}\mathit{\boldsymbol{p}}}$ : Percentage of the polymer materials r that is recyclable from the products of automotive plastic pieces p.

${\mathit{\boldsymbol{\alpha }}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}}$ : Percentage of the demand products of automotive plastic pieces p from the customers of plastic car accessories c for the purpose of recycling.

${\mathit{\boldsymbol{\gamma }}}_{\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}}$ : The amount of pollution from the vehicle v per unit of the products of automotive plastic pieces p.

${\mathit{\boldsymbol{\gamma }}}_{\mathit{\boldsymbol{v}}\mathit{\boldsymbol{r}}}^{'}$ : The amount of pollution caused by the vehicle v per unit of the polymer materials r.

$\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The capacity reduction rate caused by disruption for the factory j and products of automotive plastic pieces p within the period t per scenario s.

$\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The capacity reduction rate caused by disruption for the supplier of polymer materials i, for the polymer material r, within the period t per scenario s.

$\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The capacity reduction rate caused by disruption for the distributors of automotive plastic accessories k for products of automotive plastic pieces p within the period t per scenario s.

$\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The capacity reduction rate caused by disruption for the recycling center h for products of automotive plastic pieces p within the period t per scenario s.

${\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{t}}}$ : The inflation rate during period t.

${\mathit{\boldsymbol{P}}}_{\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The selling price of the products of automotive plastic pieces p, within period t.

$\mathit{\boldsymbol{C}}{\mathit{\boldsymbol{S}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The shortage cost of the products of automotive plastic pieces p at the factory j during period t.

$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The shortage cost of the product at the distributors of automotive plastic accessories k during period t.

$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}}$ : The shortage cost of the polymer material r at the recycling center h during period t.

$\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{d}}}_{\mathit{\boldsymbol{t}}}$ : The sum of the budget for each period t.

${\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}^{'}$ : The quality of the products of automotive plastic pieces p for distributors of automotive plastic accessories k during period t.

${\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ : The minimum demanded quality of products of automotive plastic pieces p for customers of plastic car accessories c within the period t.

3.1.5. Robust variables

${\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : Inventory of the products of automotive plastic pieces p at the factory j during period t per scenario s.

${\mathit{\boldsymbol{S}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The shortage of the products of automotive plastic pieces p at the factory j during period t per scenario s.

$\mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The inventory of the products of automotive plastic pieces p at the distributors of automotive plastic accessories k during period t per scenario s.

$\mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The inventory of the polymer material r at the recycling center h during period t per scenario s.

$\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The shortage of the products of automotive plastic pieces p at the distributors of automotive plastic accessories k during period t per scenario s.

$\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : The shortage of polymer materials r in the recycling center h during period t per scenario s.

${\mathit{\boldsymbol{x}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : Shipment volume of polymer materials r from the supplier of polymer materials i to the factory j during period t using the vehicle v per scenario s.

$\mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : Shipment volume of the polymer materials r from recycling center h to the factory j during period t using the vehicle v per scenario s.

${\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : Returned amount of products of automotive plastic pieces p from customers of plastic car accessories c to recycling center h during period t using the vehicle v per scenario s.

${\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}$ : Amount of the products of automotive plastic pieces p from distribution center k to customers of plastic car accessories c during period using vehicle v per scenario s.

${\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}$ : Amount of the products of automotive plastic pieces p by factory center j to distributors of automotive plastic accessories k during period t using vehicle v per scenario s.

${\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}$ : The demand of the customers of plastic car accessories c during period t based on the inflation rate per scenario s.

${\mathit{\boldsymbol{\xi }}}_{\mathit{\boldsymbol{s}}}$ : The total value of the objective functions for every single scenario irrespective of its probability.

3.1.6. Binary variables

$Y{I_i}, Y{K_k}, Y{H_h}, Y{J_j}, y{y_{kcpts}}, {z_{kcpts}}.$

3.1.7. Robust decision variables

${I_{jpts}} , {S_{jpts}} , I{K_{kpts}} , I{H_{hrts}} , S{K_{kpts}} , S{H_{hrts}} , {x_{ijrvts}} , z{z_{hjrvts}} , {q_{chvpts}} , {y'_{kcpvts}} , {D'_{cpts}} , {y_{jkpvts}} \geqslant 0 ,$

$y{y_{kcpts}} , {Z_{kcpts}} , Y{I_i} , Y{K_k} , Y{J_j} , Y{H_h} \in \left\{ {0, 1} \right\} .$

3.1.8. Modeling the objective functions

Optimum maintenance is the first part of the objective function and calculates the costs of inventory and shortage in the factory, distributors, and recycling centers and the shipping cost of the product to the factory, distributors, customers, and recycling centers during each period.

$F{1}_{ts} = \sum\limits_{j\in J}\sum\limits_{p\in P}\left({I}_{jpts}.C{I}_{jpt}+{S}_{jpts}.C{S}_{jpt}\right)$

$+\sum\limits_{k\in K}\sum\limits_{p\in P}\left(I{K}_{kpts}.CI{K}_{kpt}+S{K}_{kpts}.CS{K}_{kpt}\right)$

$+\sum\limits_{h\in H}\sum\limits_{r\in R}\left(I{H}_{hrts}.CI{H}_{hrt}+S{H}_{hrts}.CS{H}_{hrt}\right)$

$+\sum\limits_{j\in J}\sum\limits_{i\in I}\sum\limits_{v\in V}\sum\limits_{r\in R}{x}_{ijrvts}.C{x}_{ijrvt}$

$+\sum\limits_{j\in J}\sum\limits_{h\in H}\sum\limits_{v\in V}\sum\limits_{r\in R}z{z}_{hjrvts}.Cz{z}_{hjrvt}$

$+\sum\limits_{c\in C}\sum\limits_{h\in H}\sum\limits_{v\in V}\sum\limits_{p\in P}{q}_{chvpts}.C{q}_{chvpt}$

$+\sum\limits_{j\in J}\sum\limits_{k\in K}\sum\limits_{v\in V}\sum\limits_{p\in P}{y}_{jkpvts}.C{y}_{jkpvt}$

$+{\sum }_{c\in C}{\sum }_{k\in K}{\sum }_{v\in V}{\sum }_{p\in P}{y}_{kcpvts}^{'}.C{y}_{kcpvt}^{'} , \forall s\in S, t\in T .$

The maintainability of the models is the second part of the objective function and calculates the volume of environmental pollution generated during the shipping process under each scenario.

$F{2}_{s} = \sum\limits_{j\in J}\sum\limits_{h\in H}\sum\limits_{v\in V}\sum\limits_{r\in R}\sum\limits_{t\in T}{\gamma }_{vr}^{'}.z{z}_{hjrvts}$

$+\sum\limits_{c\in C}\sum\limits_{h\in H}\sum\limits_{v\in V}\sum\limits_{p\in P}\sum\limits_{t\in T}{\gamma }_{vp}.{q}_{chvpts}$

$+\sum\limits_{j\in J}\sum\limits_{i\in I}\sum\limits_{v\in V}\sum\limits_{r\in R}\sum\limits_{t\in T}{\gamma }_{vr}^{'}.{x}_{ijrvts}$

$+\sum\limits_{j\in J}\sum\limits_{k\in K}\sum\limits_{v\in V}\sum\limits_{p\in P}\sum\limits_{t\in T}{y}_{jkpvts}.{\gamma }_{vp}$

$+{\sum }_{c\in C}{\sum }_{k\in K}{\sum }_{v\in V}{\sum }_{p\in P}{\sum }_{t\in T}{y}_{kcpvts}^{'}.{\gamma }_{vp} , \forall s \in S, t \in T.$

Overall objective function 1 minimizes the amount of system costs for each scenario, as below:

$Min\left(A{1}_{s}\right) = {\sum }_{t\in T}\left(F{1}_{ts}\right), \forall s\in S .$

Objective function 2 minimizes the amount of environmental pollution per scenario, as below:

$Min\left(A{2}_{s}\right) = F{2}_{s} .$

If the weighted summation technique is used, the sum of the first and second objective functions is as follows:

${\xi }_{s} = \left({w}_{1}.A{1}_{s}+{w}_{2}.A{2}_{s}\right)\forall s\in S .$

3.1.9. Modeling the constraints

Constraint 1 calculates the shortage, inventory, shipment, and receipt at the factory for each scenario.

$\begin{array}{l} {I}_{jpts}-{S}_{jpts} = {I}_{jpt-1s}-{S}_{jpt-1s}+\sum\limits_{r\in R}\sum\limits_{v\in V}\sum\limits_{i\in I}{x}_{ijrvts}.{\lambda }_{rp} \\ -{\sum }_{k\in K}{\sum }_{v\in V}{y}_{jkpvts}+{\sum }_{r\in R}{\sum }_{v\in V}{\sum }_{h\in H}z{z}_{hjrvt-1s}.{\lambda }_{rp} , \forall \mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. \end{array}$

(3.1)

Constraint 2 calculates the capacity of the supplier of polymer material in terms of disruption for each scenario.

${\sum }_{v\in V}{\sum }_{j\in J}{x}_{ijrvts}-{\eta }_{irts}^{1}\le (1-\beta {I}_{irts}).cap{I}_{irt}.Y{I}_{is} , \forall \mathit{\boldsymbol{i}}\in \mathit{\boldsymbol{I}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{r}}\in \mathit{\boldsymbol{R}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} .$

(3.2)

Constraint 3 calculates the factory capacity in terms of disruption for each scenario.

${\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{2}}}\le ({\bf{1}}-\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}).\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}.\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}} , \forall \mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.3)

Constraint 4 calculates the capacity of plastic accessory distributors in terms of disruption for each scenario.

$\begin{array}{c} {\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}-{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{3}}}\le ({\bf{1}}-\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}).\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}.\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}} , \\ \forall \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. \end{array}$

(3.4)

Constraint 5 calculates the maximum shipment to the recycling centers of plastic auto pieces in terms of disruption for each scenario.

$\begin{array}{l} \sum\limits_{v\in V}\sum\limits_{c\in C}{q}_{chpvts}-{\eta }_{hpts}^{4}\le \left(1-\beta {H}_{hpts}\right).cap{H}_{hpt}.Y{H}_{h}, \\ \ \ \ \ \ \ \ \ \ \forall \mathit{\boldsymbol{h}}\in \mathit{\boldsymbol{H}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. \end{array}$

(3.5)

Constraint 6 calculates the shortage, inventory, shipment, and receipt at the distributors of automotive plastic accessories for every scenario.

$\begin{array}{l} \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} = \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}} \\ \ \ \ \ \ \ \ \ \ \ \ \ \ +\sum\limits_{\mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\sum\limits_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}, \forall \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. \end{array}$

(3.6)

Constraint 7 calculates the amount of shipments to customers for every scenario.

${\sum }_{\mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}}{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} = {\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} , \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} .$

(3.7)

Constraint 8 calculates the shortage, inventory, and shipment of the polymer materials at the recycling center for each scenario.

$\begin{array}{l} \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} = \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}\sum\limits_{\mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}}\mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}+\sum\limits_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}\sum\limits_{\mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{\varphi }}}_{\mathit{\boldsymbol{r}}\mathit{\boldsymbol{p}}}.{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}, \\ \ \ \ \ \ \ \ \ \ \forall h\in \mathit{\boldsymbol{H}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{r}}\in \mathit{\boldsymbol{R}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} . \end{array}$

(3.8)

Constraint 9 determines the quantity demanded in accordance with the price and inflation rate for every scenario.

${\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}+{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{5}}} = {\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-({\bf{1}}+{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{t}}}).{\mathit{\boldsymbol{P}}}_{\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}} , \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} .$

(3.9)

Constraint 10 determines the costs for each scenario based on the budget.

$F{1}_{ts}\le B{d}_{t} .$

(3.10)

Constraints 11 and 12 calculate the binary variable for every scenario.

$\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\le {\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} , \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.11)

$\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}.\mathit{\boldsymbol{M}}\mathit{\boldsymbol{M}}\ge {\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} , \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} .$

(3.12)

For each scenario, constraints 13 and 14 determine, that if the quality of the shipped product is lower than the minimum quality, it is allowed to be sent to the recycling location.

$({\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}-{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}^{'}).\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\ge ({\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\bf{1}}).\mathit{\boldsymbol{M}}\mathit{\boldsymbol{M}} , \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.13)

$({\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}-{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}^{'}).\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\le {\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}.\mathit{\boldsymbol{M}}\mathit{\boldsymbol{M}} , \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} .$

(3.14)

Constraints 15 and 16 determine the returns of plastic accessories to the recycling center based on low-quality products for every scenario.

(3.15)

$\begin{array}{c} {\sum }_{\mathit{\boldsymbol{h}}\in \mathit{\boldsymbol{h}}}{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\le {\mathit{\boldsymbol{\alpha }}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}}.{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}+\left({\bf{1}}-{\mathit{\boldsymbol{Z}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\right).\mathit{\boldsymbol{M}}\mathit{\boldsymbol{M}} , \\ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} . \\ {\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , {\mathit{\boldsymbol{S}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}, \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , \mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , {\mathit{\boldsymbol{x}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , \mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , {\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , {\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} , {\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} , {\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\ge {\bf{0}} , \\ \mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , {\mathit{\boldsymbol{Z}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{i}}} , \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}} , \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}} , \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}}\in \left\{{\bf{0}}, {\bf{1}}\right\} . \end{array}$

(3.16)

Consequently, according to the technique suggested by Mulvey ^[18] to deal whit this type of uncertainty, the above model is reformulated in the next section using the following two new parameters and two positive variables.

Parameters:

pp_s: The probability of each scenario.

ω: The unjustifiability weight of the model.

Positive variables:

θ_s′: The difference between the objective function from its mean for each scenario.

${\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{1}}}, {\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{2}}}, {\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{3}}}, {\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{4}}}, {\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{5}}}$ : The constraint violation is due to uncertainty in the parameter (indicating that the model is not justified).

3.2. Finalized model

Mulvey's ^[18] objective function can be rewritten as follows:

$\begin{array}{l} \ \ \ \ Min\left(Z\right) = \sum\limits_{{s}^{'}\in S}p{p}_{{s}^{'}}{\xi }_{{s}^{'}}+\gamma .\sum\limits_{s\in S}p{p}_{{s}^{'}}\left[\right({\xi }_{{s}^{'}}-p{p}_{s}{\xi }_{s})+2.{\theta }_{{s}^{'}}] \\ \ \ \ \ \ \ \ \ \ \ \ \ +\omega .\sum\limits_{{s}^{'}\in S}p{p}_{{s}^{'}}.(\sum\limits_{r\in R}\sum\limits_{i\in I}\sum\limits_{t\in T}{\eta }_{irts}^{1}+\sum\limits_{j\in J}\sum\limits_{p\in P}\sum\limits_{t\in T}{\eta }_{jpts}^{2}+\sum\limits_{k\in K}\sum\limits_{p\in P}\sum\limits_{t\in T}{\eta }_{kpts}^{3}\\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ +{\sum }_{h\in H}{\sum }_{p\in P}{\sum }_{t\in T}{\eta }_{hpts}^{4}+{\sum }_{c\in C}{\sum }_{p\in P}{\sum }_{t\in T}{\eta }_{cpts}^{5}) .\\ {\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}}{\mathit{\boldsymbol{x}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{1}\le (1-\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}).\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}}.\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{s}}}\ \ \forall \mathit{\boldsymbol{i}}\in \mathit{\boldsymbol{I}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{r}}\in \mathit{\boldsymbol{R}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. , \end{array}$

(3.17)

${\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{2}}}\le ({\bf{1}}-\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}).\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}.\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{J}}}_{\mathit{\boldsymbol{j}}} , \ \ \forall \mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.18)

${\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}-{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{3}}}\le ({\bf{1}}-\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}).\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}.\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}} , \ \ \forall \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.19)

$\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}\sum\limits_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{4}}}\le \left({\bf{1}}-\mathit{\boldsymbol{\beta }}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\right).\mathit{\boldsymbol{c}}\mathit{\boldsymbol{a}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}.\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}}, \ \ \forall \mathit{\boldsymbol{h}}\in \mathit{\boldsymbol{H}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.20)

${\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}+{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{{\bf{5}}} = {\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\left({\bf{1}}+{\mathit{\boldsymbol{\eta }}}_{\mathit{\boldsymbol{t}}}\right).{\mathit{\boldsymbol{P}}}_{\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}, \ \ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.21)

${\mathit{\boldsymbol{\xi }}}_{{\mathit{\boldsymbol{s}}}^{'}}-\sum\limits_{\mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}}\mathit{\boldsymbol{p}}{\mathit{\boldsymbol{p}}}_{\mathit{\boldsymbol{s}}}{\mathit{\boldsymbol{\xi }}}_{\mathit{\boldsymbol{s}}}+{\mathit{\boldsymbol{\theta }}}_{{\mathit{\boldsymbol{s}}}^{'}}\ge {\bf{0}}, \ \ \forall {\mathit{\boldsymbol{s}}}^{'}\in \mathit{\boldsymbol{S}}.$

(3.22)

$\begin{array}{l} {\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\mathit{\boldsymbol{S}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} = {\mathit{\boldsymbol{I}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-{\mathit{\boldsymbol{S}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}+\sum\limits_{\mathit{\boldsymbol{r}}\in \mathit{\boldsymbol{R}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}\sum\limits_{\mathit{\boldsymbol{i}}\in \mathit{\boldsymbol{I}}}{\mathit{\boldsymbol{x}}}_{\mathit{\boldsymbol{i}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}{.\mathit{\boldsymbol{\lambda }}}_{\mathit{\boldsymbol{r}}\mathit{\boldsymbol{p}}} \\ \ \ -{\sum }_{\mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}}{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}+{\sum }_{\mathit{\boldsymbol{r}}\in \mathit{\boldsymbol{R}}}{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{h}}\in \mathit{\boldsymbol{H}}}\mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}{.\mathit{\boldsymbol{\lambda }}}_{\mathit{\boldsymbol{r}}\mathit{\boldsymbol{p}}} , \\ \ \ \ \ \ \ \ \ \ \ \ \ \forall \mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. \end{array}$

(3.23)

$\mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} = \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}+\sum\limits_{\mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{j}}\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\sum\limits_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}, \ \ \ \forall \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.24)

${\sum }_{\mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}}{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} = {\mathit{\boldsymbol{D}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'} , \ \ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.25)

$\begin{array}{c} \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} = \mathit{\boldsymbol{I}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{H}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{t}}-{\bf{1}}\mathit{\boldsymbol{s}}}-{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{j}}\in \mathit{\boldsymbol{J}}}\mathit{\boldsymbol{z}}{\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{h}}\mathit{\boldsymbol{j}}\mathit{\boldsymbol{r}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}+{\sum }_{\mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}}{\sum }_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\sum }_{\mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}}{\mathit{\boldsymbol{\varphi }}}_{\mathit{\boldsymbol{r}}\mathit{\boldsymbol{p}}}.{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}} , \\ \forall \mathit{\boldsymbol{h}}\in \mathit{\boldsymbol{H}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{r}}\in \mathit{\boldsymbol{R}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}} . \end{array}$

(3.26)

$\mathit{\boldsymbol{F}}{{\bf{1}}}_{\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\le \mathit{\boldsymbol{B}}{\mathit{\boldsymbol{d}}}_{\mathit{\boldsymbol{t}}} .$

(3.27)

$\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\le \sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}, \ \ \ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.28)

$y{y}_{kcpts}.MM\ge \sum\limits_{v\in V}{y}_{kcpvts}^{'}, \ \ \ \forall c\in C, k\in K, p\in P, t\in T, s\in S.$

(3.29)

$\left({\mathit{\boldsymbol{Q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}-{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}^{'}\right).\mathit{\boldsymbol{y}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\ge \left({\mathit{\boldsymbol{z}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}-{\bf{1}}\right).\mathit{\boldsymbol{M}}\mathit{\boldsymbol{M}}, \ \ \ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.30)

$\left({Q}_{cpt}-{q}_{kpt}^{'}\right).y{y}_{kcpts}\le {z}_{kcpts}.MM, \ \ \ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.31)

$\sum\limits_{\mathit{\boldsymbol{h}}\in \mathit{\boldsymbol{H}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{q}}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{h}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\ge {\mathit{\boldsymbol{\alpha }}}_{\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}}\sum\limits_{\mathit{\boldsymbol{v}}\in \mathit{\boldsymbol{V}}}{\mathit{\boldsymbol{y}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{v}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}^{'}-\left({\bf{1}}-{\mathit{\boldsymbol{Z}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{c}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}\mathit{\boldsymbol{s}}}\right).\mathit{\boldsymbol{M}}\mathit{\boldsymbol{M}}, \ \ \ \forall \mathit{\boldsymbol{c}}\in \mathit{\boldsymbol{C}}, \mathit{\boldsymbol{k}}\in \mathit{\boldsymbol{K}}, \mathit{\boldsymbol{p}}\in \mathit{\boldsymbol{P}}, \mathit{\boldsymbol{t}}\in \mathit{\boldsymbol{T}}, \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}.$

(3.32)

${\sum }_{h\in h}{\sum }_{v\in V}{q}_{chpvts}\le {\alpha }_{cp}.{\sum }_{v\in V}{y}_{kcpvts}^{'}+\left(1-{Z}_{kcpts}\right).MM , \ \ \ \forall c\in C, k\in K, p\in P, t\in T, s\in S.$

(3.33)

$\begin{array}{c} {\mathit{\boldsymbol{\xi }}}_{\mathit{\boldsymbol{s}}} = \left({\mathit{\boldsymbol{w}}}_{{\bf{1}}}.\mathit{\boldsymbol{A}}{{\bf{1}}}_{\mathit{\boldsymbol{s}}}+{\mathit{\boldsymbol{w}}}_{{\bf{2}}}.\mathit{\boldsymbol{A}}{{\bf{2}}}_{\mathit{\boldsymbol{s}}}\right), \ \ \ \forall \mathit{\boldsymbol{s}}\in \mathit{\boldsymbol{S}}. \\ {I}_{jpts} , {S}_{jpts} , I{K}_{kpts} , I{H}_{hrts} , S{K}_{kpts} , S{H}_{hrts} , {x}_{ijrvts} , z{z}_{hjrvts} , {q}_{chvpts} , {y}_{kcpvts}^{'} , {D}_{cpts}^{'} , {y}_{jkpvts}\ge 0 , \\ y{y}_{kcpts} , z{z}_{kcpts} , Y{I}_{i} , Y{K}_{k} , Y{J}_{j} , Y{H}_{h}\in \left\{\mathrm{0, 1}\right\} , \\ {\theta }_{{s}^{'}}, {\eta }_{irts}^{1}, {\eta }_{jpts}^{2}, {\eta }_{kpts}^{3}, {\eta }_{hpts}^{4}, {\eta }_{cpts}^{5}\ge 0 . \end{array}$

(3.34)

4. Validation and analysis

Initially, the deterministic model was validated using examples and random data (Figure 2) to ensure that the model aligned with the system behavior.

Figure 2. Validation of the random example set.

DownLoad: Full-Size Img PowerPoint

The factory received a total of 987 tons of polymer materials from suppliers 1, 3, and 4. It produced 848,000 plastic auto pieces but delivered 817,000 pieces, to the distributors of automotive plastic accessories by vehicle 2 to deliver to customers of plastic car accessories. The distributor of automotive plastic accessories sent the received pieces by vehicle 3 to the customers of plastic car accessories and stored the rest in the warehouse in the following way: The distributor delivered 235,963 pieces to customer 1 in period 3,203,963 pieces to customer 2 in period 3, and 230,963 pieces to customer 3 in period 3, considering that the distributors of automotive plastic accessories already had a shortage of 528,000 pieces. Out of the 817,000 pieces received, the distributors of automotive plastic accessories still had 146,011 pieces, which were added to their warehouse. It should be noted that the distributor's previous inventory was 178,000 pieces, indicating the model's validity.

A sensitivity analysis was conducted based on the shortage costs at the automotive plastic accessories distributor center $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ due to the role of customer losses and the credibility of the distributors. The sensitivity of parameters $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{\mathit{\boldsymbol{k}}\mathit{\boldsymbol{p}}\mathit{\boldsymbol{t}}}$ for period 3 (t = 3), distributor ${\mathit{\boldsymbol{k}}}_{1}$ , and product ${\mathit{\boldsymbol{p}}}_{1}$ were examined, and the results are presented in Table 2.

Table 2. The sensitivity analysis based on randomness data.

Parameter	Variations in shortage cost $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{1}}}{\mathit{\boldsymbol{p}}}_{{\bf{1}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$	Shortage rate $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{1}}}{\mathit{\boldsymbol{p}}}_{{\bf{1}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$	A1	A2	Weighted sum of the objective function
	+50%	0	9675400	5409,983	4835000
$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{1}}}{\mathit{\boldsymbol{p}}}_{{\bf{1}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$	0	426,962	9663000	5392,468	4828800
	-50%	652,943	9516200	5314,817	4755500
	-70%	1074,537	9423900	4930,272	4709500

| Show Table

DownLoad: CSV

According to the decided parameters and with the increment in deficiency costs due to the increment in costs, the deficiency desire for that center ought to be decreased. Also, due to the reduction of shortage costs and, due to the state of zero changes, a shortage value has been considered for it, so the most optimal state is to increase the shortage for that period and the distributor due to the reduction of the shortage cost for the distribution center ${\mathit{\boldsymbol{k}}}_{{\bf{1}}}$ . As a result, according to the outputs of the above table, this amount of change in the variables and objective functions is consistent with the described and expected behavior, which indicates that it is correct.

Also, according , the sensitivity of the model was analyzed in real conditions to ensure the reliability and performance of the system. A sensitivity analysis was conducted based on the shortage costs at the automotive plastic accessories distributor center $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{2}}}{\mathit{\boldsymbol{p}}}_{{\bf{9}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$ due to the role of customer losses and the credibility of the distributors. The sensitivity of parameter $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{2}}}{\mathit{\boldsymbol{p}}}_{{\bf{9}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$ for period 3 (t = 3), distributor k₂, and product p₉ was examined and the results are presented in Table 3.

Table 3. The sensitivity analysis of the uncertainty model subject to real-world conditions.

Parameter	Variations in shortage cost $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{2}}}{\mathit{\boldsymbol{p}}}_{{\bf{9}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$	Shortage rate $\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{2}}}{\mathit{\boldsymbol{p}}}_{{\bf{9}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$	A1	A2	Weighted sum of the objective function
$\mathit{\boldsymbol{C}}\mathit{\boldsymbol{S}}{\mathit{\boldsymbol{K}}}_{{\mathit{\boldsymbol{k}}}_{{\bf{2}}}{\mathit{\boldsymbol{p}}}_{{\bf{9}}}{\mathit{\boldsymbol{t}}}_{{\bf{3}}}}$	+50%	0	2.3242e+8	41300.696	1.1623e+8
	-70%	683.020	2.3174e+8	41303.445	1.1589e+8
	-80%	737.190	2.3018e+8	41239.445	1.1511e+8

| Show Table

DownLoad: CSV

The increased shortage costs and the described parameters necessitate decreasing the expected shortage at the distribution centre to prevent high costs. A decreasing trend is observed in the second objective function (41300.696) as the prices increase by 50%, but the value of the first objective function (2.3242e+8) is improved, and the number of shortages becomes zero.

Nonetheless, as the costs increase by 50%, the weighted sum of the objective functions increases to a reliable high value, and reducing shortage costs by 70% leads to an increase in the shortage rate (683.020). However, compared to when the prices are increased, a decrease is observed in the value of the first objective function (2.3174e+8) and the value of the second objective function (41303.445). In addition, the weighted sum of the second and first objective functions declines to 1.1589E+8. Decreasing the shortage costs by 80% leads to an increased shortage rate of 737.190, and the dual (41239.445) and first (2.3018e+8) objective functions. Also, the weighted sum of the same objective functions (1.1511e+8) indicates a declining trend compared to increased costs. Therefore, the optimal action is to increase the shortage for that period and for the distributor to reduce the shortage costs. Therefore, based on the outputs presented in Table 2, the variations observed in the objective functions and variables agree with the predicted and defined behavior, reflecting the model's accuracy. As a result, eco-friendly automotive companies can modify their strategies according to the mathematical model developed in this investigation to be prepared for quick reactions in risky and challenging conditions to prevent bankruptcy and, avoid the purchase of low-quality and non-eco-friendly materials with low monetary values to prevent endangering the environment.

In this research, the Taguchi method, a method for designing experiments, is used to reduce the number of experiments and adjust the input parameters in the design algorithm. Utilizing the Taguchi method, estimating various response functions based on the specified factors is possible. These estimated results aid in identifying the factors that yield the optimal results for the given experiment. Overall, this method decreases the number of tests needed for optimization and enhances the accuracy of the results. Mulvey's model is a multi-product, multi-period sustainable supply chain model. Initially, it was solved on a small scale using the GAMS software and CPLEX to demonstrate its correct functionality. Subsequently, MATLAB software and the NSGA-Ⅱ and HB algorithms were employed to solve the model on a larger scale due to their high efficiency. Before solving it, the solution parameters were determined using the Taguchi parameter setting method. The model was then solved on a large scale, and the results were recorded. The impacts of adjusting the parameters with the Taguchi method are given in Figures 3–6.

Figure 3. HBA mean of means.

DownLoad: Full-Size Img PowerPoint

Figure 4. Signal noise of HBA.

DownLoad: Full-Size Img PowerPoint

Figure 5. NSGA-Ⅱ mean of means.

DownLoad: Full-Size Img PowerPoint

Figure 6. Signal noise diagram for NSGA-Ⅱ.

DownLoad: Full-Size Img PowerPoint

After validating the parameters using the Taguchi method in Minitab 16, the model was implemented in natural conditions at Gibor Sepehr Automotive Company, an auto parts supplier. This company works in the field of the design and production of all kinds of plastic and metal parts, the design and manufacture of all types of industrial models, the design and production of all kinds of car parts, the production of all types of features, the agricultural industry, the import and export of all raw materials related to the purchase of manufactured products, and the sale of all kinds of parts and industrial machines. The shipping cost parameter from the polymer material supplier to the factory was analyzed after data normalization and removal of disruptive outliers. As shown in Figure 7, this parameter consists of an intermediate layer with 25 neurons. As shown in Figure 7, this parameter consists of an intermediate layer with 25 neurons. Also, shallow neural network perceptron and the Levenberg-Marquardt machine learning algorithm, and three-time periods and MSE performance (i.e., mean squared errors of the performance function) are used. In this parameter used, the train-validation-test data division ratio is 20/10/75 and NARNET neural network and sigmoid activation function. After six iterations, the NN model stabilized with MSE = 4.05*10⁻²¹ and with very good qualitative performance.

Figure 7. Epochs of shipping costs from the supplier of polymer materials to the factory.

DownLoad: Full-Size Img PowerPoint

As shown in Figure 8, the gradient reached a minimum value of 1.08*10⁻¹¹. The average shipping costs from the supplier of polymer materials to the factory after six iterations (1*10⁻⁰⁹) achieved its minimum value and showed a downward trend.

Figure 8. Gradient and mean of shipping costs from the supplier of polymer materials to the factory.

DownLoad: Full-Size Img PowerPoint

Next the parameter of shipping costs from the factory to the distributors of automotive plastic accessories was analyzed. This parameter (Figure 9) employed the perceptron of a shallow neural network, a middle layer with 20 neurons, three time periods, the Levenberg-Marquardt machine learning algorithm, train-validation-test data split ratio of 85/5/10, the MSE quality performance, a sigmoid activation function, and a NARNET neural network. The neural network model stabilized after 12 iterations with MSE = 4.24*10⁻¹⁷, and has excellent qualitative performance.

Figure 9. Epochs of shipping costs from the distributors of automotive plastic accessories to the customers of plastic car accessories.

DownLoad: Full-Size Img PowerPoint

As shown in Figure 10, the gradient reached a minimum of 4.40*10⁻⁰⁹. The average shipping costs from the factory to the distributors of automotive plastic accessories (1*10⁻⁰⁹) in 12 repetitions passed zero and had a downward trend.

Figure 10. Mean and gradient of shipping costs from the distributors of automotive plastic accessories to the customers of plastic car accessories.

DownLoad: Full-Size Img PowerPoint

Next, the parameter of shipping costs from the distributors of automotive plastic accessories to the customers of plastic car accessories was analyzed. As shown in Figure 11, this parameter employed the perceptron of a shallow neural network, a middle layer with 27 neurons, three time periods, the Levenberg-Marquardt machine learning algorithm, a train-validation-test data split ratio of 85/10/5, the MSE quality performance, a sigmoid activation function, and a NARNET neural network. The neural network model stabilized in eight iterations with MSE = (3.32*10²⁵) qualitatively performed very well.

Figure 11. Epochs of shipping costs from the distributors of automotive plastic accessories to the customers of plastic car accessories.

DownLoad: Full-Size Img PowerPoint

According to Figure 12, the gradient reached its minimum of 6.34*10⁻¹³. Mean shipping costs from the factory to the distributors of automotive plastic accessories (1.00*10⁻¹¹) in eight repetitions passed zero and had a downward trend.

Figure 12. Mean and gradient of shipping costs from the distributors of automotive plastic accessories to the customers of plastic car accessories.

DownLoad: Full-Size Img PowerPoint

Next, the parameter of shipping costs from the customers of plastic car accessories to the recycling center was analyzed. As shown in Figure 13, this parameter employed the perceptron of a shallow neural network, a middle layer with 42 neurons, three time periods, the Levenberg-Marquardt machine learning algorithm, a train-validation-test data split ratio of 75/15/10, the MSE quality performance, a sigmoid activation function, and a NARNET neural network. The neural network model stabilized in four iterations with MSE = (1.08*10⁻²²), and qualitatively performed very well.

Figure 13. Epochs of shipping costs from the customers of plastic car accessories to the recycling center.

DownLoad: Full-Size Img PowerPoint

As shown in Figure 14, the gradient reached a minimum of 3.148*10⁻¹². The mean of shipping costs from the customers of plastic car accessories to the recycling center (1.00*10^-7) in four repetitions passed 0.

Figure 14. Mean and gradient of shipping costs from the customers of plastic car accessories to the recycling center.

DownLoad: Full-Size Img PowerPoint

The parameter of shipping costs from the recycling center to the factory was analyzed next. As shown in Figure 15, this parameter employed the perceptron of a shallow neural network, a middle layer with 38 neurons, three time periods, the Levenberg-Marquardt machine learning algorithm, a train-validation-test data split ratio of 45/35/20, the MSE quality performance, a sigmoid activation function, and a NARNET neural network. The neural network model stabilized in three iterations with MSE = (5.32*10⁻¹⁶) and qualitatively performed very well.

Figure 15. Epochs of shipping costs from the plastic pieces recycling center to the factory the auto plastic pieces.

DownLoad: Full-Size Img PowerPoint

As shown in Figure 16, the gradient reached a minimum of 3.04*10⁻⁰⁸. Mean shipping costs from the recycling center to the factory (1.00*10⁻⁰⁶), in three repetitions passed 0.

Figure 16. Mean and gradient of shipping costs from the plastic pieces recycling center to the factory for the auto plastic pieces.

DownLoad: Full-Size Img PowerPoint

After predicting the AI parameters, the data were entered into The GAMS software, and the outputs of the robust optimization model and meta-heuristic HBA and NSGA-Ⅱ algorithms were obtained. The Pareto front of the robust ε-constraint model is shown in Figure 17.

Figure 17. The Pareto front of the robust ε-constraint model.

DownLoad: Full-Size Img PowerPoint

Table 4 presents the optimal, best, and worst solutions of the GAMS objective functions compared to the solutions obtained from the HBA and NSGA-Ⅱ algorithms, highlighting their respective best and worst solutions. Table 4 compares the performance of the employed algorithms mentioned above.

Table 4. The Pareto front of the objective functions.

	ε-constraint		HBA		NSGA-Ⅱ
	ff1	ff2	ff1	ff2	ff1	ff2
1	4.37e+04	1.81e+08	1627223532.2	20082.333	1605120815.82	20145.0415
2	4.83e+04	1.35e+08	1678270448.7	20001.9491	1840651942.59	19126.8662
3	5436.73	9.03e+07	1787959645.8	19941.6493	1667663152.29	19653.1567
4	162115.564	451e+07	1789224333.8	19864.7403	1778997220.47	19217.8833
5	73394.876	0			1704746697.55	19423.0377
6	40031.743	2.26e+08			1721771354.47	19263.8599
Best	0	40031.748	1.6272e+09	1.9865e+0	1.6051e+09	1.9127e+04
Worst	2.257432e+8	73394.876	1.7892e+09	2.0082e+0	1.8407e+09	2.0145e+04

| Show Table

DownLoad: CSV

Overall, this study indicated the effectiveness and feasibility of the presented model in achieving high-quality solutions within shorter processing times. After ensuring model validity and tuning the parameters of the meta-heuristic algorithm, AI was used. The NAR neural network was employed for time series forecasting and successfully predicted the parameters of shipping costs, surpassing other techniques in performance. Next, the mathematical robust optimization model (presented by Mulvey ^[19]) was used. As shown in Table 4, the worst point in the first objective function is 2.257432e+8, and the worst point in the second objective function is 73394.876. In the HBA, the worst point of the first objective function was 1.7892e+09, and the worst point of the second objective function was 2.0082e+04. Based on the results, the best solution obtained with the HBA indicated that, in the best case, the minimums of the second (1.9865e+04) and first (1.6272e+09) objective functions dominate other points of the Pareto front.

NSGA-Ⅱ algorithm results (see Table 4) indicate that the worst point of the first objective function was 1.8407e+09, and the worst point of the second objective function was 2.0145e+04. According to the best solution obtained with the NSGA-Ⅱ algorithm, in the best case, the minimums of the second (1.9127e+04) and first (1.6051e+09) objective functions are dominant over other points of the Pareto front. In the HBA, the worst solution on the Pareto front contains the maximums of the second (2.0082e+04) and first (1.7892e+09) objective functions, and in the NSGA- Ⅱ algorithm, the worst solution on the Pareto front contains the maximums of the second (2.0145e+04) and first (1.8407e+09) objective functions, which are the non-dominated points (see Table 4).

Comparing the Pareto fronts of the robust optimization model and the HBA and NSGA-Ⅱ algorithms showed that, in the robust optimization model, point (40031.748) dominates other points because the goal is minimization. In the HBA and NSGA-Ⅱ algorithms, the Pareto front points 1.9865e+04 and 1.9127e+04 dominate other points, respectively. It should be noted that, since the goal of the model is minimization, the lowest point dominates the higher points and is selected.

Then, Mulvey's ε-constraint model, the HBA algorithm, and the NSGA-Ⅱ algorithm were compared based on performance indicators (see Table 5). The findings indicated that the MID of the HBA algorithm (0.0604) had a shorter distance than the ε-constraint model (0.789) and the NSGA-Ⅱ algorithm (0.0822). In terms of performance of the indicator of computation time, Mulvey's model took a shorter time (0.327) to solve the problems than the NSGA-Ⅱ (203.0374) and HB (140.7682) algorithms.

Table 5. Performance indicators of the algorithms.

	DM	SNS	NPS	MID	Time
Robust ε-constraint, GAMS	2.257433e+08	1.378*	6*	0.789	0.327*
HBA	1.6200e+08	0.0457	4	0.0604*	140.7682
NSGA-Ⅱ	2.3553e+08*	0.0382	6*	0.0822	203.0374
Note: * indicates optimal.

| Show Table

DownLoad: CSV

With regard to the number of Pareto solutions (NPS), Mulvey's model (6) and the NSGA-Ⅱ algorithm showed better performance than the HBA (4). The SNS was higher in Mulvey's model (1.378) than the NSGA-Ⅱ (0.0382) and HB (0.0457) algorithms. This finding demonstrates the more robust performance of this indicator. Additionally, a higher DM value was obtained for the NSGA-Ⅱ algorithm (2.3553e+08) than Mulvey's model (2.257432e+8) and the HBA (1.6200e+08). In general, the mentioned comparisons indicate that the hybrid model features a suitable structure for problem-solving in uncertain conditions, and eco-friendly companies can use this model to modify their necessary strategies.

Model and algorithm performance indicators: Spread of non-dominated solutions (SNS), Diverty Metric (DM), Mean Ideal Distance (MID), and Number of Pareto Solutions (NPS), are given below:

Finally, the following questions of the study were answered:

1) What are the parameters causing uncertainty in the mathematical model of the sustainable supply chain? Uncertainty parameters in this model are shipping costs estimated by the neural network and the demand and capacity reduction due to disruption.

2) What are the robust parameters in the sustainable supply chain model? The non-deterministic demand, the capacity reduction caused by disruption in the factory, the capacity reduction caused by disruption in the supplier, the capacity reduction caused by disruption in the distributors, the capacity reduction caused by disruption in the recycling.

3) Under uncertain conditions, which of the parameters of the sustainable supply chain model can be predicted using the AI approach? Parameters: shipping costs from supplier to factory, from factory to distributor, from distributor to customer, from customer to recycling.

5. Conclusions

The present study designed a sustainable supply chain model to minimize system costs and environmental pollutants. In the proposed two-objective model, a new hybrid approach composed of AI, a robust optimization model, and the HBA and NSGA-Ⅱ meta-heuristic algorithms were used to address the uncertain nature of the problem and deal with the uncertainties. The hybrid model of the study can significantly improve the economic and environmental aspects. Additionally, the proposed model can be an efficient tool for designing a stable and flexible supply chain network and making related decisions. Experimental results and analysis demonstrate the effectiveness of Mulvey's model, the HBA, and the NSGA-Ⅱ algorithm in solving optimization problems with complex and challenging search spaces. This hybrid model proved effective in overcoming problem complexity and solving large-scale examples with high processing quality. Heuristic solutions also demonstrated a significant reduction in computation time for large problems. The results of the AI approach showed its superiority over other techniques, particularly in predicting shipping cost fluctuations in conditions of uncertainty. The robustness-based model effectively dealt with demand uncertainties, transportation costs, and disruptions that reduced capacity. Even though the robust model exhibited limitations and could not solve problems with larger sets the HBA and NSGA-Ⅱ algorithms were used to address these larger-scale problems. Experimental results showed the effectiveness of the HBA and NSGA-Ⅱ algorithms in solving NP-hard optimization problems, surpassing other methods in convergence speed. The analysis of results indicated that the presented hybrid model features a suitable structure for problem-solving in uncertain conditions, and eco-friendly companies can use this model to modify their necessary strategies. These results also highlight their superiority in terms of convergence speed compared to other methods. Furthermore, the results indicate the effectiveness and capability of implementing the presented model, leading to high-quality solutions in less processing time. This research helps decision-makers in making quick predictions and developing appropriate strategies in uncertain situations. Moreover, our findings affirm that the proposed model holds great significance in ensuring the sustainability of the supply chain. Based on sensitivity analysis, the level of accountability greatly influences the goals of a sustainable supply chain, and managers should consider the trade-off between accountability and their objectives. We provided practical suggestions for improving sustainable supply chains:

● To moderate environmental pollution caused by companies, different measures can be taken which include implementing sustainable practices, adopting cleaner technologies, reducing waste generation, recycling, treating wastewater before disposal, and adhering to ecological regulations and standards. Additionally, companies can promote environmental awareness, engage in responsible corporate practices, and invest in renewable energy sources to reduce their overall environmental impact.

● To have significant impact on the world community and improve air quality, the ultimate goal should be to reduce greenhouse gas emissions and promote sustainable practices across industries. This can involve a combination of measures, including transitioning to renewable energy sources, implementing stricter environmental regulations, and promoting the use of green materials in manufacturing processes.

● By transitioning to green materials, automotive companies can reduce their carbon footprint and minimize the negative environmental impact of plastic production and disposal. Green materials, such as recycled plastics or bio-based alternatives, have a lower carbon footprint and can be more easily recycled or biodegraded. While there may be upfront costs associated with adopting green materials, automotive companies can achieve profitability by complying with these standards in several ways.

● Cost savings in the long run: Green materials can be more cost-effective in the long term, as they often require less energy and resources. Additionally, using recycled materials can reduce the need for raw material extraction.

● Market differentiation: Companies that prioritize sustainability, and utilize green materials can pull in environmentally-conscious shoppers, driving them to make strides in deals and brand dependability.

● Access to green markets: Some regions or countries may have specific requirements or preferences for eco-friendly products. By complying with these standards, automotive companies can gain access to these markets and expand their customer base.

Also, governments can implement stricter rules and standards on using green materials in automotive manufacturing, which can include requirements for the percentage of recycled or bio-based materials used in plastic parts. Governments can provide financial incentives and support to automotive companies that adopt green materials in their manufacturing processes, which can offset the additional costs associated with using environmentally friendly materials. Increasing consumer awareness and demand for eco-friendly products can push automotive companies to prioritize green materials. Companies are more likely to adopt sustainable practices if they see a market for such products. Further suggestions for future research are provided:

● It is suggested that, in future research, all costs (production, inventory, warehouse, salaries, marketing, advertising, inventory control, etc.) should be considered as uncertainty parameters in the model, and only certain costs should not be considered.

● It is suggested that, according to the current research model, in AI, all uncertain parameters should be used for forecasting. This is because, in uncertain conditions, rapidly forecasting uncertainties will speed up problem-solving and timely strategy formulation by managers.

● It is suggested that future researchers study several larger automobile companies that play a prominent role in environmental pollution. This will help take better and more important steps to minimize environmental pollution. Additionally, this supply chain sustainability problem cannot be solved simply by studying one supplier or one company.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Acknowledgments

The authors would like to acknowledge the valuable contributions of the respected chief editor, anonymous referees, and editors. Their insightful comments, feedback and suggestions have significantly improved the quality of this work.

Conflict of interest

The authors declare there is no conflict of interest.

References

[1]	L. Angermann, S. Wang, Convergence of a fitted finite volume method for European and American option valuation, Numer. Math., 106 (2007), 1–40. https://doi.org/10.1007/s00211-006-0057-7 doi: 10.1007/s00211-006-0057-7
[2]	A. Bensoussan, J. L. Lions, Applications of Variational Inequalities in Stochastic Control, North-Holland Amsterdam, New York, Oxford, 1978.
[3]	T. B. Gyulov, M. N. Koleva, Penalty method for indifference pricing of American option in a liquidity switching market, Appl. Numer. Math., 172 (2022), 525–545. https://doi.org/10.1016/j.apnum.2021.11.002 doi: 10.1016/j.apnum.2021.11.002
[4]	J. Haslinger, M. Miettinen, Finite Element Method for Hemivariational Inequalities, Kluwer Academic Publisher, 1999. https://doi.org/10.1007/978-1-4757-5233-5
[5]	L. Scurria, D. Fauconnier, P. Jiranek, T. Tamarozzi, A Galerkin/hyper-reduction technique to reduce steady-state elastohydrodynamic line contact problems, Comput. Methods Appl. Mech. Eng., 386 (2021), 114132. https://doi.org/10.1016/j.cma.2021.114132 doi: 10.1016/j.cma.2021.114132
[6]	W. D. Zhao, Numerical methods for forward backward stochastic differential equations, Math. Numer. Sin., 37 (2015), 337–373. https://doi.org/10.12286/jssx.2015.4.337 doi: 10.12286/jssx.2015.4.337
[7]	L. Jiang, Convexity, translation invariance and subadditivity for g-expectations and related risk measures, Ann. Appl. Probab., 18 (2008), 245–258. https://doi.org/10.1214/105051607000000294 doi: 10.1214/105051607000000294
[8]	I. Penner, A. Reveillac, Risk measures for processes and BSDEs, Finance Stochastics, 19 (2015), 23–66. https://doi.org/10.1007/s00780-014-0243-x doi: 10.1007/s00780-014-0243-x
[9]	R. Glowinski, Numerical Methods for Nonlinear Variational Problems, Springer-Verlag, 1984. https://doi.org/10.1007/978-3-662-12613-4
[10]	S. Wang, X. Q. Yang, K. L. Teo, Power penalty method for a linear complementarity problem arising from American option valuation, J. Optim. Theory Appl., 129 (2006), 227–254. https://doi.org/10.1007/s10957-006-9062-3 doi: 10.1007/s10957-006-9062-3
[11]	S. Wang, C. S. Huang, A power penalty method for solving a nonlinear parabolic complementarity problem, Nonlinear Anal. Theory Methods Appl., 69 (2008), 1125–1137. https://doi.org/10.1016/j.na.2007.06.014 doi: 10.1016/j.na.2007.06.014
[12]	M. Chen, C. Huang, A power penalty method for a class of linearly constrained variational inequality, J. Ind. Manage. Optim., 14 (2018), 1381–1396. https://doi.org/10.3934/jimo.2018012 doi: 10.3934/jimo.2018012
[13]	A. M. Rubinov, X. Q. Yang, Lagrange-type Functions in Constrained Non-convex Optimization, Kluwer Academic Publishers, Dordrecht, Holland, 2003.
[14]	X. Q. Yang, X. X. Huang, Nonlinear lagrangian approach to constrained optimization problems, SIAM J. Optim., 11 (2001), 1119–1144. https://doi.org/10.1137/S1052623400371806 doi: 10.1137/S1052623400371806
[15]	M. H. Holmes, Introduction to Mumerical Methods in Differential Equations, Springer New York, 2009.
[16]	O. T. Hanna, New explicit and implicit "improved Euler" methods for the integration of ordinary differential equations, Comput. Chem. Eng., 12 (1988), 1083–1086. https://doi.org/10.1016/0098-1354(88)87030-3 doi: 10.1016/0098-1354(88)87030-3
[17]	R. Santos, L. Alves, A comparative analysis of explicit, IMEX and implicit strong stability preserving Runge-Kutta schemes, Appl. Numer. Math., 159 (2021), 204–220. https://doi.org/10.1016/j.apnum.2020.09.007 doi: 10.1016/j.apnum.2020.09.007
[18]	E. Zeidler, Nonlinear Functional Analysis and Its Applications Ⅱ/B: Nonlinear Monotone Operators, Springer-Verlag, New York, 1990.
[19]	Y. L. Zhao, P. Y. Zhu, X. M. Gu, X. L. Zhao, H. Y. Jian, A preconditioning technique for all-at-once system from the nonlinear tempered fractional diffusion equation, J. Sci. Comput., 83 (2020), 10. https://doi.org/10.1007/s10915-020-01193-1 doi: 10.1007/s10915-020-01193-1
[20]	G. H. Golub, C. F. Van Loan, Matrix Computations, Johns Hopkins University Press, 2012.
[21]	C. M. da Fonseca, On the eigenvalues of some tridiagonal matrices, J. Comput. Appl. Math., 200 (2007), 283–286. https://doi.org/10.1016/j.cam.2005.08.047 doi: 10.1016/j.cam.2005.08.047
[22]	A. R. Willms, Analytic results for the eigenvalues of certain tridiagonal matrices, SIAM J. Matrix Anal. Appl., 30 (2008), 639–656. https://doi.org/10.1137/070695411 doi: 10.1137/070695411
[23]	W. Luo, X. M. Gu, B. Carpentieri, A hybrid triangulation method for banded linear systems, Math. Comput. Simul., 194 (2022), 97–108. https://doi.org/10.1016/j.matcom.2021.11.012 doi: 10.1016/j.matcom.2021.11.012

This article has been cited by:

Massoumeh Nazari, Mahmoud Dehghan Nayeri, Kiamars Fathi Hafshjani, Correction: Developing mathematical models and intelligent sustainable supply chains by uncertain parameters and algorithms, 2024, 9, 2473-6988, 25223, 10.3934/math.20241230

Reader Comments

Your name:*

Email:*
© 2023 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)

通讯作者: 陈斌, bchen63@163.com

1.
沈阳化工大学材料科学与工程学院沈阳 110142

Electronic Research Archive

1.1 1.7

Metrics

Article views(1931) PDF downloads(76) Cited by(0)

Preview PDF

Download XML

Export Citation

Article outline

Show full outline

Figures and Tables

Figures(8)

Electronic Research Archive

A numerical method for parabolic complementarity problem

Related Papers:

Abstract

1. Introduction

2. Literature review

3. Methodology

3.1. Mathematical model design

3.1.1. Robust indexes

3.1.2. AI parameters

3.1.3. Conceptual model parameters in condition certain

3.1.4. Robust parameters

3.1.5. Robust variables

3.1.6. Binary variables

3.1.7. Robust decision variables

3.1.8. Modeling the objective functions

3.1.9. Modeling the constraints

3.2. Finalized model

4. Validation and analysis

5. Conclusions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Catalog

Electronic Research Archive

A numerical method for parabolic complementarity problem

Related Papers:

Abstract

1. Introduction

2. Literature review

3. Methodology

3.1. Mathematical model design

3.1.1. Robust indexes

3.1.2. AI parameters

3.1.3. Conceptual model parameters in condition certain

3.1.4. Robust parameters

3.1.5. Robust variables

3.1.6. Binary variables

3.1.7. Robust decision variables

3.1.8. Modeling the objective functions

3.1.9. Modeling the constraints

3.2. Finalized model

4. Validation and analysis

5. Conclusions

Use of AI tools declaration

Acknowledgments

Conflict of interest

References

This article has been cited by:

Reader Comments

通讯作者: 陈斌, bchen63@163.com

Metrics

Figures and Tables

Other Articles By Authors

Related pages

Tools

Export File

Citation

Format

Content

Catalog