Solving the Sustainable Automobile Production-Distribution Joint Optimization in the Physical Internet-Enabled Hyperconnected Order-to-Delivery System by I-NSGAIII

The Physical Internet (PI)-enabled hyperconnected order-to-delivery system (OTD) provides new solutions for sustainable supply chains from production perspectives. In this system, a PI-enabled hyperconnected manufacturing system is more closely tied with other functions through Internet-of-Things (IoT)-enabled machines for communication. In the OTD, the PI-enabled hyperconnected production–distribution system (PI-H) is modelled by multi-objective mixed-integer-nonlinear programming to evaluate sustainability. We develop an improved reference-point based non-dominated sorting genetic algorithm (I-NSGAIII) to solve practical-scale PI-enabled hyperconnected production-distribution scheduling problems, with the problem-specific solution expression and dynamic programming, subproblem-guided crossover and mutation strategies, and adaptive evolution mechanisms. I-NSGAIII’s performance advantages and PI-H’s sustainable advantages are validated through extensive experiments.


I. INTRODUCTION
For auto enterprises, implementing sustainable supply chain management (SSCM) within the framework of government regulations on environmental policies has the highest priority [4]. The order-to-delivery system (OTD) serves as a supply chain platform, connecting the market and manufacturer. This platform helps manage the production status and delivery progress of orders in production and transportation systems [2]. Hence, sustainably producing and delivering customised vehicles to customers in shorter delivery times is the primary purpose of OTD. Meanwhile, the Internet-of-Things (IoT) revolution promotes the Physical Internet (PI) paradigm in OTD, where the end-to-end visibility of the PI objects, operations and systems is realised [3]. A Physical The associate editor coordinating the review of this manuscript and approving it for publication was Shih-Wei Lin .
Internet-enabled hyperconnected OTD brings a new application scenario of a digital-driven manufacturing system that is more conducive to achieving sustainability in jointly scheduling production and distribution.
Auto supply chains with highly modular production processes have huge potential in digital supply chain implementation [2]. Ji et al. [4] studied the application of PI and found that a PI-enabled auto supply chain (PI-S) can provide a better service level with less cost compared to a traditional auto supply chain (TS), without considering sustainability. However, their study only investigated the advantages of transparency from the tracking function of PI, owing to existing limitations in the infrastructure and ignoring the hyperconnected distribution brought by PI. PI can digitally connect and share information in real-time among physical machines, assets and components, generating the operation scenario with an IoT-enabled machine-to-machine (IoT-M2M) [5]. VOLUME 11, 2023 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Our article formulates a model to evaluate the sustainability performance of the operation scenario in IoT-M2M, including the M2M-based automated modular roll spraying system and the hyperconnected distribution. In detailed, the M2M-based automated modular roll spraying system with multiple nozzles inside the robot arm is adopted to form the digital-driven painting workshop(PT). Colours in the nozzle constitute a colour package that massively changes the production logic and human participation compared with the machine-driven PT. Additionally, the hyperconnected distribution of resource flows between plants provides more flexibility for the auto enterprise, requiring unique production facility configurations and costly large-scale equipment [6]. By implementing the PI-enabled hyperconnected production-distribution system (PI-H), these factors promote the formation of hyperconnected OTD for auto enterprises. Moreover, many scholars have conducted extensive qualitative analysis, arguing that PI-M2M will fully and optimally leverage machines while liberating the workforce and improving resource utilisationwhich are two key sustainability indicators [7]. Hence, the potential positive influence brought by PI-H on sustainability is quite worthy of study. As defined in the Johannesburg Summit, sustainable development should balance not just economy and environmental protection but also social development [8]. These gaps between practical requirements of sustainability management in the auto supply chain and existing research allow us to explore implementable PI-H and further investigate its sustainability advantages from economic, environmental and social dimensions. The integrated production-distribution scheduling (IPDS) has been extensively studied in pursuing a more economical supply chain. Moreover, the interconnection, dynamics and openness brought by PI make IPDS has the characteristic of hyperconnection, w.r.t, PI-HPDS. In this article, the PI-HPDS is formulated as a multiple-objective model to comprehensively evaluate its sustainability from the economic, environmental, and social aspects. Moreover, an I-NSGAIII is proposed to solve the practical-scale problem based on the characteristics of production and distribution in PI-HPDS. The PI-H, PI-S, and TS are investigated under different make-to-order market (MTO) scenarios, combining different customer sizes and market structures. Notably, PI-H has superior performance on sustainability in all instances compared to PI-S and TS. Furthermore, PI-H can intensify these advantages over PI-S as market structures' complexity increases. Finally, the practical value of adopting PI-H for auto-enterprises is discussed. The contributions of this article are summarised as follows: • The PI-enabled hyperconnected production-distribution system is modelled as MOMINLP to assess its sustainability performance. The M2M-based modular roll spraying system and hyperconnected distribution of resource flow provide a new production scenario for the PI-HPDS. • We develop the I-NSGAIII to solve practical-scale problems more effectively. In detail, problem-specific solution expression and dynamic programming, subproblem-guided crossover and mutation operators are proposed. In addition, its performance is verified via statistical comparisons.
• Under the different market structures and demand scales which are the practical situation in MTO, the sustainability performance of PI-H is quantitatively investigated, and some managerial implications are established.
The remaining contents are structured as follows. Section II summarises relevant literature. Section III provides a detailed description of the model. Section IV describes the I-NSGAIII algorithm. Section V performs extensive computational experiments to discuss the efficiency of I-NSGAIII and the performance of the PI-enabled hyperconnected supply chain. Finally, Section VI presents our conclusions and some directions for future work.

II. LITERATURE REVIEW
Current research on jointly scheduling production and distribution prioritises the sustainable performance of the economy, environment and society instead of previous research that has focused on the economy alone. Chandra and Fisher [9] were the first to discuss the huge economic advantage brought by IPDS. Subsequently, IPDS only targeting economic optimisation has been extensively studied under different production configurations. Some researchers assumed that production is performed by a single machine [10]. IPDS with a parallel machine configuration, wherein a job is assigned to multiple machines, has also been studied in [11]. Moreover, Yagmur and Kesen [12] and Homayouni, et al. [13] studied IPDS with a more complex production configuration (e.g. single-stage flow-shop and job-shop). Ji et al. [4] further investigated the cost optimisation of IPDS with a multi-stage manufacturing system (MMS). Through a limited buffer system, jobs are then transferred sequentially in two workshops. Jobs cannot be produced in subsystem i + 1 until subsystem i completes the production and transfer of related jobs. The above studies only focus on the economy. Subsequently, Aghajani and Al-E-Hashem [14], Al-E-Hashem et al. [15], Solina and Mirabelli [16] study how IPDS can consider both environmental and economic goals. Particularly, Al-E-Hashem et al. [15] only focused on carbon emissions in transportation, while Aghajani and Al-E-Hashem [14] and Solina and Mirabelli [16] further studied energy consumption on a single production line. Peng et al. [22] evaluated the sustainability performance of PI-enabled production-inventory-distribution system from economic, environmental and social dimensions. They found that PI will bring better sustainability performance for the supply chain. However, they did not extensively study the social impact of production (e.g. the worker's rights from worker safety and fair production), which is receiving increasing attention. The aforementioned research gaps drive us to explore the sustainability performance of PI-enabled IPDS in more practical production scenarios.
Montreuil [17] developed the concept of PI, a physical metaphor for the digital internet. PI promotes the interconnection and openness of independent supply networks to obtain significant advantages [17]. Moreover, Crainic and Montreuil [18] emphasised hyperconnectivity in logistics and supply chains. Ben Mohamed et al. [19] optimised interconnected urban logistics with economic and ecological attributes. Moreover, Chargui et al. [20] and Hu et al. [21] studied the interconnected distribution system. Peng et al. [22] and Ji et al. [23] successively evaluated PI-enabled production-inventory-distribution system based on cost and sustainability performance. Regardless, the above research simplified the production process as the order allocation problem. Ji et al. [4] further studied IPDS with PIenabled MMS production configuration, which considers job scheduling in sequentially connected two workshops. However, in this article, the two plants independently perform production, and the significance of production for sustainable management is not explored. Furthermore, Marcotte, Montreuil and Coelho [6] firstly introduced the hyperconnected mobile production system. A PI-enabled hyperconnected production system not only relocates processing modules based on demand but also dynamically transfers resources between open plants. Notably, digitally driven manufacturing brought by IoT-M2M in Industry 4.0 considerably changes production scenes in existing MMS and open plants, including the M2M-based modular roll spraying system and the hyperconnected semi-finished delivering system. The research gaps between current research and practical requirements for auto enterprises motivated us to study PI-HPDS based on the PIenabled hyperconnected production-distribution system and evaluate its sustainability.
In recent years, a lot of approaches have been developed to address IPDS, such as the exact method [24] and heuristic algorithms (e.g. adaptive genetic algorithm [10], memetic algorithm [12] and NSGAIII [25]). Exact methods are not undeniably feasible when the problem is scaled up to a certain level. Moreover, most practical problems for enterprises involve multi-objective optimisation. Hence, multiobjective evolutionary algorithms (MOEA) have become suitable solution methods. MOEA based on the Pareto concept has been considered more effective compared to the other methods, including the weighted sum of objective functions and approaches based on population [26]. For instance, strength Pareto evolutionary algorithms (SPEA and SPEA2) have been exposed [27], and non-dominated sorting genetic algorithms (NSGA and NSGA-II) are proposed [28]. Deb and Jain [29] propose the reference-point-based NSGAIII following the NSGA-II framework, which is much better than a classical generating method in many-objective problems. However, NSGAIII does not perform well in largescale optimisation problems. Wang [30] combines several information feedback models with NSGAIII to improve its ability to solve large-scale optimisation problems. However, the above works disregard the mostly nonlinear nature of real problems. In nonlinear problems, the minimal evolutionary fluctuations in the decision space will cause great distribution changes in the target space. The information entropy theory is suited for solving the PI-HPDS. Table 1 summarises the research gaps in the related literature.

III. MODEL FORMULATION A. PROBLEM DEFINITION
The proposed PI-HPDS can realise hyperconnected production between two PI-enabled plants and then deliver finishedvehicles to retailers through a multi-level hyperconnected distribution network (Figure 1). The green section describes the hyperconnected production system. Every PI-enabled plant has multiple workshops. Semi-finished products can be delivered between two plants to complete the remaining steps of the production process. Hereafter, only finished-vehicles, which are off the assembly line, will be delivered to customers through a two-level PI-enabled hyperconnected distribution system with RFID, WSN and GPS to perform tracking and tracing. After the finished vehicle leaves the plant, it can be transported to multiple PI-hubs or directly to the retailer. The blue section illustrates the first layer in the transportation system, where trains provide service between plants and hubs. The second layer, marked in yellow, covers routing between PI-hubs and retailers. The PI-HPDS can be defined as a directed network G = (V , E). While V is the vertex set, and E is the edge set, represent arcs between nodes in the vertex set V . Vertex set V includes the plant set P, virtual stock point set P ′ , hub set H and retailer set R. The virtual stock point set P ′ and hub set H constitute N 1 = P ′ ∪ H . The time horizon is T (t ∈ {1, 2, . . . , |T |}), each t in T is the uniformly sized discrete period, and the length of each unit time period can be user-defined when doing time discretisation.
In PI-HPDS, Original equipment manufacturers (OEMs) perform the hyperconnected semi-finished transport for car-shells before the painting shop (PT). Orders are initially pre-allocated to different plants based on the practical constraints of the press workshop only. Hence, the body-inwhite in Plant A will be sent to Plant B for painting and assembly under the whole process of tracking and accurate identification. Meanwhile, the digital painting workshop adopts an M2M-based modular roll spraying system, which has multiple nozzles inside the robot arm. In this study, the total number of spray nozzles κ n = 3 and colours in spray nozzles constitute a colour package. If the colour required at the current tact is contained by the colour package, changing the nozzle becomes unnecessary. Notably, spray nozzles can be replaced one or two at a time if they do not work in the current tact. When two spray nozzles are replaced in the same tact, pollution emission becomes less than that of two nozzles replaced separately in a different tack. Otherwise, at least one of the spray nozzles will be replaced in the previous production tact. A special organic solvent is needed to clean the paint material, which would generate a lot of VOC, which is harmful to the environment and human beings. Hence, minimising the times of colour switching is required, as spray gun cleaning is performed whenever the colour needs to be changed [31]. Moreover, two types of parts can be found in the M2M-AS: standard and key parts. This article only considers the key parts (G) which are affected by the assembly sequence. Smoothing the part demand [32] is obtained in the first function to reduce demand fluctuations in each period. Meanwhile, the assembly load balancing is formulated as a constraint satisfaction problem [33].
To balance the sustainability of economic, environmental and social dimensions, a MOMINLP model is proposed to formulate this problem. Notably, the production also has a major influence on the environment and society. Consequently, not only sustainability from transportation is considered but also the environmental and social effects of production. In the environmental objective, VOC and CO 2 emission costs are counted. The social objective considers the social welfare contribution of carbon mitigation and workers' rights from workers' safe and fair production.
The sets, parameters and decision variables are defined as the following notation list: The economic objective function (1) minimises the sum of fixed setup cost and variable cost for production, inventory cost, loading cost, transportation cost and delivery tardiness cost. The variable production cost includes the colour switching cost for PT and the supply cost for AS.

2) ENVIRONMENTAL OBJECTIVE
In the environmental objective, the cost of GHG emissions caused by painting and transportation are considered, as there is no energy consumption reduction during the smoothing part demand for AS. The first and the second terms are GHG consumption for different transportation modes ( [34], [35]). Moreover, the third term calculates the process cost of VOC [36].

3) SOCIAL OBJECTIVE
Regarding social objective, the objective function (3), as shown at the bottom of the page, is evaluated from three aspects: the social contribution of carbon mitigation, workers' safety concerns [37] and production fairness [38]. One part of the social contribution of carbon mitigation comes from PI-enabled transportation. The other part comes from the recycling system of VOC in PT, which is generated during the cleaning of nozzles. Moreover, the workers' safety concern is contributed by the cost-saving of colour switching between the traditional PT and M2M-PT [39]. Fair production is calculated based on the Coefficient of Variation for the production period between plants.

C. GENERIC MODEL
The sustainability of PI-HPDS can be presented as min f Sus = (f Eoc , f Env , −f Soc ). The constraints are as (4)- (18), shown at the bottom of page 7. The assembly schedule must satisfy the limitations of production capacity, production sequence, and production bottlenecks. Constraints (4) indicate that the order a must not be assembled in the plant i at the period t when the plant i does not perform the production at the period t. Constraints (4) also restrict the production capacity of AS for every plant in each period. Constraints (5) further stipulate the transfer of semi-finished products between the plants. Constraints (6) define the Coefficient of Variation for fair production. Constraints (7) specify that each order must be produced at most once. The assembly sequence is represented by constraints (8) - (15), where virtual orders 0 and |O| + 1 are set as the head and tail orders for each time period in constraints (9) and (10). Constricts (11) and (12) define the cross-period connections for the assembly sequence. Constraints (13) and t∈T i∈P VOLUME 11, 2023 specify that the maximum number of key parts G selected in a continuously produced Q g order is q g . Constraints (16) make sure that the installation of each workstation is completed within the tact. Constraints (17) and (18) (24), as shown at the bottom of the next page, illustrate the logical relationship of colour switching in M2M-PT. Constraints (19) specify the replacement of colour e at the time slot nt is performed only when it is in the spray package during nt-1. Constraints (20) specify that if the colour e is replaced by colour d at the time slot nt, then colour d must be included in the spray package at nt+1. Constraints (21) stipulate that if not being switched, the spray package's colours are still used at the next time slot. Constraints (22) ensure the maximum number of colours in the spray package, and constraints (23) restrict that at least one colour should not be switched. Constraints (24) state that order a cannot be painted on time slot nt if colour e of order a is not currently in the spray package. Constraints (25), as shown at the bottom of the next page, indicate that an order must be painted and assembled in the same plant. Constraints (26) - (29), as shown at the bottom of the next page, describe the production capacity for PT in each time period t, where the total time period T can also be counted by time slot nt. The production sequence constraints in PT are stated as constraints (30) - (37), as shown at the bottom of the next page. Constraints (38), as shown at the bottom of the next page, specify that order a cannot start to be assembled before the time that painting operations are performed and the transferring time between CBS. Constraints (39), as shown at the bottom of the next page, describe the logic of colour switching in typical PT.
Transportation-related constraints are described above. Constraints (40), as shown at the bottom of page 10, specify that each vehicle must be delivered to a virtual plant or hub after being assembled. Constraints (41), as shown at the bottom of page 10, restrict the number of finished-vehicles that can be delivered. Constraints (42), as shown at the bottom of page 10, forbid orders to revisit any arcs. Constraints (43) - (45), as shown at the bottom of page 10, are flow balance constraints. Constraints (46), as shown at the bottom of page 10, specify the departure time of orders. Constraints (47) -(51), as shown at the bottom of page 10, define the shipping time and inventory. Notability, the stock of plants is placed in the virtual node P ′ , so it has zero stock.

IV. SOLVING PROCEDURES
As stated above, PI-HPDS is formulated as a MOMINLP model, simultaneously optimising production in the hyperconnected production system and vehicle routing problem in the hyperconnected distribution system. This model would also be NP-hard. The solution representation and the operators of crossover and mutation based on the characteristics of PI-H are proposed for a more effective progeny population generation. NSGAIII is used as the basic framework of the developed algorithm. Moreover, dynamic programming is proposed to get the optimum scheme for the painting subproblem, improving the quality of the solution and the performance of the algorithm. The PI-HPDS proposed in this article is a nonlinear problem, and even if small evolutionary fluctuations in the decision space may lead to significant changes in the target space. Therefore, we simultaneously analyse the population from both decision space and target space and propose adaptive adjustment mechanisms. The adaptive information feedback mechanism and the adaptive reference point-based selection mechanism are adopted into I-NSGAIII, balancing the convergence and the diversity of the algorithm.

A. SOLUTION REPRESENTATION
Solution representation contributes to mapping PI-HPDS's solution space to I-NSGAIII's search space. Since the intermediate production and transportation need to be determined simultaneously, they partition the incidence matrix A = λ ij , i = 1, . . . , 2 · |T | + 2, j ∈ N into submatrices A 1 , A 2 , and A 3 . Thus we can take A = (A 1 , A 2 , A 3 ), in which A 1 , A 2 , and A 3 represent painting sequencing, assembly sequencing, and transportation scheme. Obviously, λ 1j and λ 2j denotes order at index j of the painting and assembly sequencing, respectively. If the order at index j of A 1 is the same as the order at index k of A 2 , then j and k should be subject to a relation that j < k − ℓ h + ϕ s . The visited route and transportation mode of each order partition the matrix A 3 into T-paired two-rowed submatrices , j ∈ N describes the transportation for the order at index j of A 2 in n th period, including the visited route i = 2n + 1 and transportation state i = 2n + 2. A 2 is set as the order index because only the assembly scheme directly affects subsequent transport scheduling. A simple example with the size T × N = 4 × 5 is as follows: which describes that order 4 is painted in first and assembly in first. Subsequently, order 4 is delivered through routes 3 and 8.

B. SUBPROBLEM GUIDED DYNAMIC PROGRAMMING
According to the practical production process of PT, spray nozzles can be changed two colours at a time or one colour at a time and the cost involved varies, thus switching scheme affects the incurred costs. Moreover, the digital-driven PT can meet any production requirements of assembly, while the subsequent transport scheduling is only affected by assembly sequencing. We propose dynamic programming to get the optimum change scheme of spray nozzles in PT, which has no negative influence on assembly and transportation. It is beneficial for the quality of the solution and the performance of I-NSGAIII.
First, we reconstruct the painting sequence in each individual to obtain the compact painting sequence. We replace the consecutive identical colours in the sequence as a colour block. Then, let arr be the array of colours in the compact production sequence, and P be the set of the index of colour in the compact production sequence, and S be the set of spray nozzles, i.e., {1, 2, 3}, and E i is the set of colours in the spray colour package when the painting process reaches the colour on the index i. j 1 , j 2 , k respectively represent the vector of available sprays position, idle sprays position and the underworking spray position ∀j 1 , j 2 , k ∈ S. t ∈ P is the time when the sprays nozzle is used again. To simplify the recurrence function, we introduce a function δ + (i, j 1 , j 2 , k, E, t) to calculate the preceding state when the painting process reaches the i th colour in compact sequence, with available and unavailable spray position (j 1 , j 2 , k), colour package (E), and |T |·ϕ s nt=1 i∈P STP a = nt · λ nt ai ∀i ∈ P, ∀a ∈ O, ∀nt ∈ NT (35) |e · ξ ae · λ nt−1 ai − e · ξ ae · λ nt ai | ≤ M · ϒ nt-1 i ∀nt ∈ NT, ∀i ∈ P (39) VOLUME 11, 2023 the time of next use (t). The function is defined as (52), shown at the bottom of the next page. The piecewise linear function g i, E ′ describes the pollution emission, σ 1 and σ 2 are the pollution emission of changing one or two colours at a time, where σ 2 > σ 1 ∨ σ 2 < 2 · σ 1 as (53), shown at the bottom of the next page.
Then the recurrence function, used to calculate the total emission, can be written as: In most metaheuristic algorithms, the updating process does not utilise the historical information of individuals in the previous iterations, leading to the loss of some information conducive to population evolution. We perform the information feedback mechanism to select the previous individuals to update the current offspring, regardless of good or bad individuals. The weight of influence is determined based on the fitness function. We further propose the self-adaptively selection strategy for avoiding the premature convergence problem, based on the evolution state. The main idea is to evaluate the distribution of the current decision space using information entropy, which is then adopted to turn the selection strategy.
We adopt the affinity theory of the immune system to describe the distribution of the population. The population is composed of N individuals having M genes. According to the problem-specific solution representation, we set the set of S ={S 1 , S 2 , S 3 } to describe all the elements that occur in the corresponding gene, where S 1 is the set of order, S 2 is the set of the route, and S 3 is the set of transportation state. The population's entropy is calculated as Eq.(55). Eq.(56) calculates the information entropy of the j th gene of all individuals, where p ij = n j /N is the probability of i th element at j th gene, and n j is the size of the set associated with the j th gene. Further, the affinity of i th and j th individuals is calculated as Eq.(57), where H (2) is the entropy of i th and j th individuals. So the affinity of the population is calculated as Eq.(58), in which Q P ∈ [0, 1]. During the evolution, the population changes from disorder to order and gradually converges. At the early stage of the evolution, the value of Q P is small and with the process proceeding, the value becomes bigger.
The adaptive information feedback mechanism is given in algorithm 2.

E. SUBPROBLEM GUIDED CROSSOVER AND MUTATION OPERATORS
To realise a more effective progeny population generation, we propose a strategy to identify poor variants and improve Inversion_Mutation and Subtour_Exchange_Crossover based on the subproblem's production logic. Considering the online time of painting and assembly, this article adopts the synchronous crossover and mutation to simultaneously update the relative sequence of orders in painting and assembly. Especially for PT, randomly select two positions in the production sequence i, i ′ to perform a mutation, if the three colour blocks before and after i th position have the same colour of arr[i], while do not include colour arr[i], the mutation is not allowed. In the crossover process, the intersection segment and the three colour blocks before and after it are updated into a new colour blocks sequence. If there is the same colour in the original or the new local colours blocks sequence, we remove one of them. Next, we compare the length of the two sequences, if the new sequence is longer than the original sequence, no crossover is performed. The strategy is described in Figure 2. Due to the directionality of transportation routes, we carry out Cross_Exchange_Crossover, Or_Opt Mutation, and 2_Opt Mutation while keeping local nodes unchanged. Meanwhile,   else // calculate similarity of individuals within this reference point 10: if Q t p ≤ Q 0 then 11: if D i , j ≤ ϖ then 16: end if 20: ρ j = ρ j + 1, F l = f l \s 21: k = k + 1 22: else Z r = Z r /J 23: end if 24: while kϵK the crossover and mutation operations are performed only between paths but not within the path.
The overview of crossover and mutation operators is given in algorithm 3.

V. EXPERIMENTAL SETUP A. TEST INSTANCES AND ALGORITHMS
A large dataset is extracted in consultation with industrial managers. Expressly, fifteen instances are selected, including

Algorithm 2 The Adapitive Information Feedback Mechanism
Require: ∈, Q 0 , N , n, X t Ensure: X l+1 1: do 2: Calculate Q t p for the population X t of generation t 3: Perform Subproblem-guided crossover and mutation operators with X t to generate the offspringX t+1 4:   [29], MOEA/D [40], SPEA2 [27], and MOCELL [41] as comparison algorithms. To comprehensively evaluate the performance of I-NSGAIII, I-NSGAIII was compared with CPLEX and comparison algorithms. The relative parameters in the comparison algorithms are identical to those in I-NSGAIII. For the stopping condition of these five algorithms, the maximum number of function evaluations (NFE) is set to 50,000. The stopping condition of CPLEX is the maximum time limit of 10800s. All experiments are performed on a PC with an Intel Core-i5 3.30 GHz processor with 4 GB memory, and algorithms are coded in Java with IntelliJ IDEA 2019. The detailed data used in this article can be downloaded at https://pan.baidu.com/s/ 1dV5-KlrQcr-M7b0we6zdzw.

B. PERFORMANCE METRICS
To evaluate the algorithms' performance, we adopt three important multi-objective performance metrics, Inverted Generational Distance (IGD) [42], Hyper Volume (HV) [43] and C-metric [44]. Specifically, C-metric considers the non-dominated sets to assess the dominance of the two algorithms; IGD and HV are used to assess the convergence and distribution of the solutions. We use the Wilcoxon rank-sum test with a 95% confidence level, the nonparametric Kruskal-Wallis test with a 95% confidence level and the Friedman ranking test to determine the statistical significance of the computation results [45], [46]. C(A, B) indicates the percentage of the solutions in the set B dominated by the solutions in the set A. If C(A, B) > C(B, A), then the algorithm to obtain A is regarded better.
• IGD is calculated as IGD(PF * , P) = x∈PF * min y∈P ∥x−y∥ |PF * | , using ∥x − y∥ to calculate the Euclidean distance between solutions x and y. Notability, the true PF of the real instances studied in this article is unknown. Therefore, we aggregated the set of non-dominated solutions obtained by all algorithms during 21 independent runs and set it as the reference PF * for each instance. A smaller IGD value points to the better algorithm performance.
• HV is calculated as HV(P) = ∪ x∈P vol(x), vol(x) representing the volume of the hypercube that is enclosed by VOLUME 11, 2023

Algorithm 3 Subproblem Guided Mutation and Crossover Operators
Require: HP, P t , p pm , ppc, p tm , p tc Ensure: U t+1 1: for i = 0 to |PS| // PS is the size for P t do 2: Randomly select an individual X j from HP 3: // Perform the Subproblem guided crossover and mutation in the production section 4: if random [0, 1] < p pc then 5: Perform Subproblem guidedcrossover between production segments of X i and X j 6: end if 7: if random [0, 1] < p pm then 8: Perform Subproblem guided -mutation with production segments of X j 9: end if 10: // Perform the crossover and mutation in the transportation section 11: if random [0, 1] < p tc then 12: Perform Cross Exchange between transportion segments of X i , and X j 13: end if 14: if random [0, 1] < P tm then 15: Randomly select mutation operator R 16: switch R do 17: case oropt 18: Perform or -optmutation with transportion segments of X i 19: case 2opt 20: Perform 2optmutation with transportion segments of X i 21: end if 22: end for individual x and the reference point, which is set as (1.0, 1.0, 1.0) in this article. The better performance of the algorithm, the greater value of HV.

C. PARAMETERS TUNNING
I-NSGAIII has 5 key parameters, population size (Ps), production crossover probability (Pc), transportation crossover probability (Tc), production mutation probability (Pm), and transportation mutation probability (Tm). For obtaining the appropriate parameter combination, the standard Design of Experiment (DOE) [47] method was performed on the representative medium instance with the scale of 120 × 3 × 8. Moreover, the mean value of HV is applied as the response value (RV). We set these parameters with different levels, as Specifically, the orthogonal array with the corresponding RV is listed in Table 3. We further analyse the extreme difference in RV values between the four levels for each parameter. The larger extreme difference of a parameter implies that the algorithm is more sensitive to the parameter and the parameter is more important. It can be observed that I-NSGAIII is most sensitive to Ps and second to Pc. The trends of RV with each parameter are shown in Figure 3. Referring to the results of DOE, the key parameters of I-NSGAIII are selected as follows: Pc = 0.9, Pm = 0.15, Tc = 0.9, Tm = 0.20, and Ps = 250. The rest parameter Q o = 0.1 and ε = 0.014 are obtained by the parameter sensitivity analysis in the next subsection.

D. PARAMETER SENSITIVITY ANALYSES
Two parameters in the adaptive information feedback mechanism can be found, including the affinity threshold Q o and feedback probability ε To verify the most suitable combination of Q o and ε, we used the Friedman ranking test to investigate the effect of parameter variation on the algorithm.
We set the affinity threshold and random feedback probability as Q o = [0.1,0.12,0.14,0.16,0.18] and ε = [0.006, 0.008,0.01,0.012,0.014], respectively. For each parameter combination, I-NSGAIII is run 21 times on all instances independently. Taking the HV metric as an example, Figure 4 presents the statistical results of the average Friedman ranking. Figure 4(a) shows that when ε is fixed, the algorithm performance deteriorates as Q o increases. With Q o = 0.1 and Q o = 0.12, the algorithm has a better performance compared to other situations. By fixing Q o ∈ {0.1,0.12}, Figure 4   This is because the algorithm needs pay more attention to exploration in the early stage of evolution. Meanwhile, smaller Q o and larger ε values make both locating local optimisation and taking the random feedback method to improve the diversity of the offspring population easier. This enables the algorithm to avoid local optimisation.

A. EFFECTIVENESS OF DYNAMIC PROGRAMMING AND INFORMATION FEEDBACK MECHANISM
To verify the effectiveness of proposed strategies, in all instances, this comparison was performed between I-NSGAIII, I-NSGAIII without the subproblem-guided dynamic programming (I-NSGAIII-ND) and I-NSGAIII without the adaptive feedback mechanism (I-NSGAIII-NA). Table 3 summarises the mean and standard deviation of the IGD and HV metrics for I-NSGAIII, I-NSGAIII-ND and I-NSGAIII-NA. Moreover, the Wilcoxon's rank-sum test was performed between I-NSGAIII, I-NSGAIII-ND and I-NSGAIII-NA at 95% confidence level. Marks '+', '−' and '∼' indicate that the performance of the comparison algorithm is significantly better, worse and similar to I-NSGAIII, respectively. Dark grey and light grey shading mark optimal and suboptimal performance, respectively. Table 4 shows that, for most instances, I-NSGAIII has the best performance, and I-NSGAIII-NA obtains suboptimal  performance. The Wilcoxon rank-sum test results clearly show that I-NSGAIII significantly outperforms I-NSGAIII-ND and I-NSGAIII-NA, especially as the scale increases. For suboptimal performance, I-NSGAIII-NA obtains more times than I-NSGAIII-ND; hence, I-NSGAIII-NA is more competitive than I-NSGAIII-ND. The above results fully demonstrate that both subproblem-guided dynamic programming and the adaptive feedback mechanism positively impact the proposed algorithm. Subproblem-guided dynamic programming can more significantly improve algorithm performance.

B. COMPARISON WITH EXACT SOLVER
In this subsection, we performed a comparison of the proposed algorithm with CPLEX in small-scale instances. To solve the multi-objective optimisation problem by CPLEX, we adopt the object weighting method and select 30 uniformly distributed weight vectors in the objective space. We then obtain the non-dominated solution set. As the PF of each real-world instance is unknown, we combine nondominated solutions obtained with each execution of all algorithms to form the approximate PF to calculate performance '   metrics. Table 5 summarises the Wilcoxon's rank-sum test results of IGD. Table 5 indicates that I-NSGAIII outperforms other algorithms for IGD metric. Moreover, as the scale of the problem grows, the quality of CPLEX's solution worsens. The reason is that the number of variables and constraints increases significantly with the size of the problem. Hence, CPLEX would then add more inequality constraints in the model, making it harder to solve. Table 6 reports three indicators for evaluating the superiority of the solution obtained by I-NSGAIII and CPLEX, including the relative gap %P, %N and the average running time (T ). %P indicates the gap between the lower and upper bounds obtained by CPLEX, and %N indicates the gap between the optimal solution obtained by I-NSGAIII and the upper bound in CPLEX.
I-NSGAIII is observed to always obtain a solution that is not worse than CPLEX. For example, in the case of 24 × 2 × 3 instance with socially sustainable optimisation as the highest priority goal, the average and maximum %P are 1.83% and 2.47%, respectively, indicating that the solution is excellent [48]. Table 7 visualises the scheduling result. Thus, the feasibility of our established MOMINLP model is verified. In the remaining three cases, I-NSGAIII can achieve a near-optimal solution within the given timeframe. As the scale increases, the advantages of I-NSGAIII become more obvious than CPLEX. VOLUME 11, 2023

C. COMPARISON WITH EXISTING ALGORITHMS
In this subsection, we performed the comparisons of I-NSGAIII with the other four compared algorithms in all instances. To intuitively show the performance of algorithms, Figure 5 presents approximate Pareto fronts obtained by these five algorithms based on two different scale instances. For both instances, the solutions obtained by I-NSGAIII are better than the other compared algorithms (no matter the quality or quantity of the solutions). Hence, I-NSGAIII can obtain more reasonable and competitive Pareto fronts.
To further evaluate the effectiveness of I-NSGAIII, the Wilcoxon's rank-sum tests of IGD and HV metrics were conducted. Table 8 and Figure 6(a) show that I-NSGAIII achieves the optimal median of IGD in all instances. Table 9 and Figure 6(b) show that I-NSGAIII achieved the best HV values except for the ninth instance. In the ninth instance, the HV performance of I-NSGAIII is similar to MOEA/D but has a smaller variance. Therefore, the convergence and diversity of I-NSGAIII are better. Table 10 summarises the C-metric performance of comparison algorithms, where I-NSGAIII, NSGAIII, MOEA/D, SPEA2 and MOCELL are denoted by A, B, C, D and E, respectively. We performed the nonparametric Kruskal-Wallis test with 95% confidence level with Cmetric; Table 10 lists p values. Table 10 shows that I-NSGAIII can obtain more non-dominant solutions than other algorithms in most cases. Thus, I-NSGAIII ensures a remarkable convergence while maintaining a justifiable diversity at its approximated Pareto front. Therefore, the effectiveness of I-NSGAIII is verified. These results confirm the ability of I-NSGAIII to solve the problems encountered in the real industry.

D. COMPARISON OF SUPPLY CHAINS' PERFORMANCE
Without loss of generality, the PI-H, PI-S and TS were tested on medium-scale instances. The Pareto solutions are shown to demonstrate their performance of sustainability. Specifically, the operating characteristics of the three SCs are analysed based on key operating cost indicators. Besides, the sustainability of these three SCs is discussed in different structures of market demands and customer scales. The improvement among PI-S, PI-H and TS was estimated by cost-saving [(C − C PI )/C PI ] × 100%. Figure 7 shows their approximate Pareto front (green dots, yellow dots, and red dots represent PI-H, PI-S, and TS, respectively).   Figure 7 demonstrates that the sustainability of PI-H in the social dimension is superior to PI-S. To further analyse the PI-H's sustainability advantage on whether PI-H can enhance sustainability performance without the additional cost of PI-S, the ideal point for PI-S is set to compare with PI-H. The economic, environmental and social object of the ideal point is set as the optimal value in Case A, Case B, and Case C, respectively. The optimal values for different optimisation cases are obtained by setting the priority of the corresponding objective function higher than the other two objectives. In this case, the ideal point of PI-S is (976330, 205.2214, 7504.264) marked in black. Moreover, Figure 7 intuitively shows PI-H has multiple solutions, outperforming the ideal point.
To further clarify PI-H's advantages, the cost performance of each operating factor of the different SCs in three optimisation cases was investigated. To show the changing trend clearly, the operating factors are set as the fixed production cost (FPC); the variable production cost (VPC); the transportation cost (STC); the loading cost (SLC); the inventory holding cost (SHC); the production emission cost (SPE); the transportation emission cost (STE); the production social contribution (SPS); the transportation social contribution (STS). The payoff table in Table 11 clearly shows that PI-H achieves the best performance in every optimisation case. Figures 8-10 demonstrate the operating cost composition of each SC in three optimisation cases specifically.
Under the circumstance of A, the total cost of PI-H and PI-S are lower than that in TS 24.48% and 14.22%, respectively. Obviously, STC and the total production cost are also lower than those in TS. In PI-SCs, the M2M-PT greatly reduces the number of gun switching, and the logical constraints from production are relaxed, which provides more options  for transportation. Besides, the total emission costs of both PI-SCs are also lower than TS. However, the GHG emissions from finished-vehicles holding and loading operations in PI-SCs are higher than those in TS. Furthermore, the emission cost from holding and loading operations, and the total cost in PI-H are lower than those in PI-S. PI-H executes the hyperconnected semi-finished transportation system to realise the connection of the material flow between different workshops in different plants and further integrate plants with different locations into a uniform order pool digitally. PI-H reduces the sum production period, the number of switching in PT, the fluctuation in AS, and transportation cost compared with PI-S. However, the VPC of PI-H is higher than that in PI-S, because the unit cost of the inter-plant transport is higher than the unit cost of painting and assembly. Meanwhile, the social influence of PI-H is also significantly better than PI-S and TS. Moreover, PI-H deepens the production's social advantage through the semi-finished transportation, which balances the uneven distribution of order series, colours, etc., and reduces the discrete coefficient of the production system. Compared with PI-S and TS, PI-H performs better in all three aspects of sustainability.
By analysing solutions in environmental-oriented and social-oriented cases (Figures 8 and 9), we can easily show that differences between these three different SCs vary slightly. Whether cases are environment-or social-oriented, the costs of PI-SCs remain lower than those in TS. Notably, PI-H improves more environmental benefits at a lower cost than PI-S. Conversely, variable production cost for PI-H increases faster than PI-S. PI-H reduces GHG from colour switching by performing a semi-finished transport, bringing greater VPC compared with PI-S. As PI-S can only adjust scheduling according to orders' colour within each plant, it would cause a delay in the offline time of some orders. Subsequently, to ensure timely arrival of orders, we decreased the vehicle utilisation of the overall transportation system. Consequently, STC of PI-S is much higher than that of PI-H. Notably, PI-S still lags far behind PI-H in the social benefits aspect. This is because PI-S cannot adjust the demand fluctuations caused by MTO, resulting in a larger discrete coefficient of production than PI-H. Generally, PI-H can significantly reduce the sustainable pressure of auto enterprises at a lower cost.

E. MTO'S PARAMETER SENSITIVITY ANALYSIS
The market structure and order scale are the two essential features of MTO. This section examines whether these two features influence PI-H's advantages. We performed analyses on nine combinations of three market structures and three instance sizes. Table 12 presents the experimental results. Table 12 indicates that TS's performance is always worse than PI-SCs. Therefore, the following analysis mainly   focuses on PI-H and PI-S. It is clear that the economic and social advantages of PI-H over PI-S become more obvious in all customer scales as the complexity of the market structure increases. However, as market structure complexity increases, the environmental benefit of PI-H over PI-S decreases from 37.837%, 41.088%, 47.8196% to 37.794%, 39.633%, 42.219% and 23.239%, 26.5403%, 41,446% respectively. In the same structure, different customer scales cannot affect the social performance of PI-H and PI-S, while economic advantage of PI-H is in a negative correlation with the customer scales. The advantages are decreasing from 11.967%, 9.256%, 12.387% to 9.489%, 9.1049%, 11.464% and 9.152%, 5.820%, 11.105%. This might be because scale economy gradually improves truck utilisation. Conversely, in the same structure, the environmental advantage of PI-H is in positive correlation with customer scales, and the advantage increases from 37.79%, 37.84%, 23,239% to 39.63%, 41.09%, 26.54% and 42.22%, 47.82%, 41.45%, respectively. Increased orders in the same demand classification impose greater scheduling pressure on the assembly constraints in PI-S compared to PI-H. Notably, PI-H always has superiority in all scenarios. Thus, PI-H can address the MTO market with higher sustainable performance.

VII. MANAGERIAL INSIGHTS
In the above discussion, the sustainability performance of PI-H has been fully demonstrated to practitioners. Beginning from the practical management of the OEM, we extract three important management implications: • For the OEM, adopting PI-H is recommended if they want to realise MTO more sustainably. PI-H breaking the maximum capacity limit of a single plant brings greater single-cycle flexible scheduling capabilities to address random demands, which have an uneven distribution of order series, colours and delivery time. The hyperconnected supply chain framework is also applicable to other OEMs.
• The OEM should adopt PI-H, especially with a more complex market structure. The above discussions have clearly shown that PI-H will obtain greater economic and social benefits from flexible production capacity as market structure complexity increases. On environmental sustainability, if trucks with new energy can be used for transportation, PI-H will perform better.
• Combining the modular concept with hyperconnected semi-finished transportation will have a considerably positive impact on the overall performance of the automobile supply chain. Existing discussions have found that the transportation cost and efficiency of semi-finished transportation have a significant impact on the supply chain's long-term performance.   The PI-enabled, hyperconnected modular system can distribute a wider variety of semi-finished products at a lower cost.

VIII. CONCLUSION
This study investigated the sustainable production-distribution joint scheduling problem in the PI-enabled hyperconnected order-to-delivery system to promote sustainability, which is crucial for the automobile industry. Hence, workshops in different plants were hyperconnected to produce various series of vehicles. Subsequently, the finished-vehicles were delivered to customers through the three-level transportation system by multiple transportation modes. Re-integration was then performed in PI-hub. We formulated the multipleobjective model based on practical features of PI-HPDS with digital-driven MMS, trading off economic, environmental and social dimensions for SSCM. In this article, we developed I-NSGAIII to solve the sustainable PI-HPDS. Firstly, we tuned the parameters by DOE test and sensitivity analysis. Secondly, we analysed the effectiveness of dynamic programming and information feedback mechanism in the proposed method. Next, we compared I-NSGAIII with CPLEX in the small-scale instances. We then evaluated the performance of I-NSGAIII by the IGD, HV and C-metric in all instances, comparing with the other four comparison algorithms. Moreover, we discussed the sustainable advantages of PI-H under different optimisation preference scenarios. Finally, we performed the sensitivity analysis to determine whether the market structure and customer demand-the two essential features of MTOinfluence PI-H's advantages. Through computational experiments, we obtained the following results: • I-NSGAIII ensures a remarkable convergence and reasonable distribution on its approximated Pareto front, providing better results than comparison algorithms.
• In practice, the proposed algorithm can provide satisfactory solutions for decision-makers to help enterprises weaken their energy consumption and thus generate higher profits.
• Results demonstrate that PI-H consistently has superior sustainability performance under MTO. PI-H will provide greater competitiveness in a more complex market structure.
• If automobile companies want to pursue significant sustainable advantages under the MTO, they should establish a PI-enabled hyperconnection production and distribution system.
Although I-NSGAIII outperforms other algorithms in comparison, it can be further improved using learning strategies. In the future, the algorithm performance should train a classification model to select promising solutions. The future direction of this study is to extend the PI-enabled hyperconnected modular application in the closed-loop supply chain and investigate its resiliency in addressing stochastic disruptions brought about by natural disasters.