Approximate model and algorithms for precast supply chain scheduling problem with time-dependent transportation times

This paper focuses on the precast supply chain scheduling problem with time-dependent transportation time to minimise the total weighted tardiness (PSCSP_TDT |TWT). In the problem, an order sequence and several job sequences are to be determined simultaneously. At first, through in-depth analysis of problem structure and real data from a precast manufacturer, we approximate the problem into a three-stage order scheduling problem by combining the seven production stages into one differentiation stage, and then explore some useful properties of the schedules for the approximate problem. Subsequently, to solve the small instances for the PSCSP_TDT |TWT, we propose an approximate model-based hybrid dynamic programming and heuristic (AMHDPH) and obtain a lower bound as a by-product of the algorithm. For dealing with medium-or large instances, with considering the complexity of the problem, we propose four approximate model-based hybrid iterated greedy (AMHIG) algorithms by integration of constructive heuristics, structural properties of solutions, an iterated greedy, and a correction heuristic. Comprehensive computational results show that the AMHDPH generates tight lower bounds for small instances and solves the most of small instances to optimality within 60 seconds. Whereas the best AMHIG generates feasible solutions with an average optimality gap below 5 percent for around 70 percent instances.


Introduction
Over the past decades, the demands for precast construction are on the fast increase especially in public housing and transportation infrastructure projects. Compared with conventional construction, precast construction has the main advantages of higher industrialisation level, lower building-energy consumption and less complex work, with the improved quality of precast products in a specialised off-site facility. With time-based competition and rapid technology advancements, scheduling has been recognised by precast manufacturers as an especially important tool that can significantly lower the tardiness penalty cost, reduce inventory and improve customer satisfaction. At the same time, as precast manufacturers are increasing relying on their suppliers, e.g. mould manufacturer, or business partners, e.g. transportation companies, precast scheduling in a factory has to be extended to coordinate members across chain of suppliers.
Advanced information technology, e.g. 5G, has changed the landscape of modern supply chain management, as it facilitates communication among members CONTACT Fuli Xiong xiongfuli@xauat.edu.cn Supplemental data for this article can be accessed here. https://doi.org/10. 1080/00207543.2022.2057254 such as suppliers, manufacturers and distributors (Dolgui and Ivanov 2022). As a result, supply chain scheduling, which is crucial in supply chain management, can be realised easier than before. Supply chain scheduling arises in different production and service systems, where limited production or service capacity and order-delivery requirements necessitate the use of minimisation of total weighted tardiness to meet distinct requirements of customers and maximising total profit. In a precast supply chain environment, which includes mould manufacturing, precast components (PCs) production, and transportation (Wang and Hu 2017;Wang et al. 2019), order delays often occur because of bottleneck process, tight due date, and traffic congestion. Thus integrated scheduling of mould manufacturing, PCs production, and transportation is crucial for precast manufacturers. However, precast supply chain scheduling problem is extremely difficult because of inherent complexity and uncertainty, e.g. due to transportation congestion. Scheduling problems are generally N P-hard in view of the combinatorial nature. In precast component production factory, there are seven stages including serial stages and parallel stages. The PC production process can be seen as a variant of flow shop, which is generally N P-hard. Transportation congestion is often met in precast supply chain systems, different departure time will generate different transportation time, which further increases the difficulty of the precast supply chain scheduling problem. How to effectively schedule and coordinate activities with and across precast supply chains is becoming a major challenge. This paper focuses on the precast supply chain scheduling problem with emphasis on time-dependent transportation and minimisation of total weighted tardiness (TWT). For simplicity, the problem is denoted by PSCSP_TDT |TWT. In the problem, each order must go through three phases sequentially: mould manufacturing, PCs production, and order delivery. First, the moulds for producing PCs are prepared in the mould manufacturing phase. Subsequently, orders from customers are fabricated in PCs production phase which includes seven stages. Finally, the orders from the same customer are grouped and delivered to the customer. Thus the problem can be seen as a nine stage flow shop scheduling problem with sequential and parallel stages. The problem can be seen as a variant of the classical flow shop scheduling problem, which has been proved to be strongly NP-hardness in most of cases (Emmons and Vairaktarakis 2012;Pinedo 2016;Baker and Trietsch 2019). In the classical flow shop scheduling problem, only a job sequence or an order sequence need to be optimised. However, for the considered problem, it requires one to simultaneously determine order sequence, job sequence in each order, and optimal departure time of each order for transportation. This makes the considered problem more complex than the classical one. The objective for the problem is to minimise the TWT value. It should be noticed that, optimising the departure time of each order is equivalent to minimise the sum of waiting time and transportation time for each order. Thus one can balance the inventory costs and reduction of environmental emissions generated by transportation (Kong et al. 2018). The contributions of this study are as follows.
• We propose a novel precast supply chain order scheduling problem with time-dependent transportation times to minimise the total weighted tardiness (PSCSP_TDT|TWT). • Based on the in-depth analysis of problem structure and real data from a PCs manufacturer, we approximate the PSCSP_TDT|TWT into a three-stage order scheduling problem (APSCSP_TDT|TWT), in which the complexity of the PSCSP_TDT|TWT is greatly reduced. Then we explore some useful structural properties of the optimal schedules for the problem. • We propose an approximate-model based dynamic programming heuristic (AMHDPH) which can generate tight lower bounds for small instances and solve the most of small instances to optimality within very short CPU times. • For medium-or large instances, by integrating structural properties of the schedules, constructive heuristics, different local search methods, a mechanism of destruction and construction, and a correction heuristic, we propose four approximate-model based hybrid iterated greedy (AMHIG) algorithms. To evaluate the performances of AMHIGs, a lower bound method based on Lagrangian relaxation is developed.
To the best of our knowledge, this work is the first attempt to optimise order sequence and job sequences simultaneously in PCs supply chain environments with time-dependent transportation times. It is also the first work to formulate an approximate model and algorithms to simply the solution complexity of the PSCSP_TDT|TWT. This study can direct PC manufacturers to minimise the total delay penalties and improve customers satisfaction effectively. The proposed methods can be extended to schedule other supply chain systems with a bottleneck stage and time-independent transportation times.
The remaining sections of this paper are organised as follows. In Section 2, we review the related literature. Problem description and formulation are presented in Section 3. In Section 4, an approximate model for the PSCSP_TDT|TWT is proposed based on the real data analysis from a PC manufacturer, and then some structural properties and properties are presented. Solution methodology based on the approximate model is developed for the PSCSP_TDT|TWT in Section 5. In Section 6, the experimental results are reported and analysed. Finally, the contributions of this work are summarised and several future considerations are outlined in Section 7. For easy understanding, we list all the abbreviations in Nomenclature of Supplemental Data.

Literature review
Over the past decades, a great deal of research on the precast production scheduling problem (PPSP), the supply chain scheduling problem (SCSP), and the precast supply chain scheduling problem (PSCSP) have been published extensively worldwide. In this section, we will review the PPSP at first and then review the SCSP and the PSCSP. On the aspect of the SCSP, we focus on briefly reviewing the integrated production and transportation scheduling problem (IPTSP), which is highly related to the PSCSP_TDT |TWT.
Several studies have been focused on the PPSP to meet on-time delivery in the past 20 years (Chan and Hu 2002;Leu and Ting Hwang 2002;Ko and Wang 2011;Yang, Ma, and Wu 2016). Precast flow shop scheduling can be seen as a variant of flow shop scheduling. However, the production of prefabricated components is quite different from traditional flow-shop scheduling. First, the supply of PCs is much more complex because activities such as jobs can be processed simultaneously in some operations, e.g. steam curing and storing. Second, the physical characteristics of large, bulky, and heavy of prefabricated components, which features multi-product, small lot and customised production, must be taken into account in all production activities of PCs. The precast production scheduling problem is considered as N P-hard (Yang, Ma, and Wu 2016), which is difficult to be solved with the traditional exact approaches, e.g. branch and cut for its large-sized instances. Therefore, many studies have been conducted on metaheuristics for the precast production scheduling problems. Hu (2001, 2002) formulated a precast flow shop scheduling model to minimise total earliness and tardiness penalties. In their model, many complex production situations, such as normal and abnormal working hours, interruption and noninterruption processes, and so on, are considered. Leu and Ting Hwang (2002) developed a scheduling model for the mixed production of PCs of several construction projects, and accordingly proposed a scheduling method based on genetic algorithm (GA). Benjaoran, Dawood, and Hobbs (2005) discussed a bespoke precast multi-objective flow-shop scheduling problem to minimise machine idle time, tardiness penalty, and makespan simultaneously. Ko and Wang (2010) formulated a precast flow shop scheduling model by considering the buffer size between production stations, and then they proposed a multiple objective GA to minimise makespan and tardiness penalties simultaneously. de Athayde Prata, Pitombeira-Neto, and de Moraes Sales (2015) modelled the precast concrete beams production problem to minimise the production waste via the compact formulation of the one-dimensional multi-period cutting stock problem. Yang, Ma, and Wu (2016) analysed several resource constraints including maintenance room size, buffer size, number of moulds and number of pallets, and formulate a flow-shop scheduling model for multiple precast production lines, and then proposed a GA to solve it. The approach was validated preliminarily by comparing with traditional scheduling approaches. Wang, Hu, and Gong (2018a) developed a two-hierarchy simulation-GA hybrid model for precast production to fill the gap in simulation system design and methodology for precast production, and increase the applicability of precast production scheduling methods in real construction projects. Wang, Hu, and Gong (2018b) also systematically integrated human factors in the form of worker competence into precast production scheduling. They proposed a two-hierarchy scheduling optimisation method considering worker competence to synchronise the production scheduling with workforce development and to optimise the scheduling from both the shortterm and long-term views. Wang and Hu (2018) formulated a two-level rescheduling model for precast production with multiple production lines and develop a GA to minimise the rescheduling costs. Ma et al. (2018) proposed a GA to optimise shop rescheduling of multiple production lines for the production of PCs in flow shop. To integrate the production planning problem of heterogeneous prestressed precast beams and the cutting problem of steel bars, an integer linear programming model, similar to the extended formulation of the one-dimensional cutting stock problem, was proposed by Araujo, Bonates, and Prata (2021), and Araújo et al. (2021). In their models, demand balance constraints and some specific restrictions on the use of patterns in each mould are considered. Signorini, Alexandre de Araujo, and Mara Melega (2021) considered the integrated production planning problem of hollow-core slabs and proposed two mathematical models for the arising problem with regarding the multiple manufacturing moulds. Heuristic methods are further proposed for dealing with these models to minimise production and inventory costs. Xiong et al. (2021) considered a distributed precast flow shop scheduling problem to minimise the total earliness and tardiness. To solve the problem, they formulated an MILP model for small size instances and proposed a hybrid iterated greedy and tabu search for large size instances.
Compared with conventional production or transportation scheduling, integrated production and transportation scheduling can significantly improve the overall performance of the supply chain system. The majority research on previous IPTSP was performed in the last 20 years. Hall and Potts (2003) provided the first research on supply chain scheduling from the perspective of general supply chain. They have shown that costs can be significantly reduced through co-operation between the two stages. Chang and Lee (2004) studied an IPTSP with single machine, single customer and single vehicle transportation. Zhong, Dósa, and Tan (2007) addressed two IPTSPs in single-machine and parallelmachine environments in which each job demands different amounts of storage space during transportation. Chen (2010) provided a detailed literature survey on previous IPTSP studies in 2010. Ullrich (2013 considered an IPTSP problem in a parallel machine environment with job-dependent processing time and delivery time window. Viergutz and Knust (2014) dealed with an IPTSP in an MTO manufacturing environment with considering a single production facility and a short product lifespan. The problem is to find a selection of customers to be supplied such to maximise the total satisfied demand. Karaoğlan and Kesen (2017) studied a variant of the IPTSP that single product with limited shelf life is produced at single facility. To solve the problem, they developed a branch and cut method in which lower bound and upper bound are updated by using several valid inequalities and simulated annealing method, respectively. Sağlam and Banerjee (2018) integrated the single machine lot scheduling problem with outbound shipment costs and formulate a mixed integer programming model. Lacomme et al. (2018) discussed an IPTSP by considering multiple vehicles, specific capacity constraints and the short lifespan of products to optimise supply chains. To deal with the problem, a greedy randomised adaptive search procedure with an evolutionary local search (ELS) is proposed. Liu et al. (2020) proposed a variable neighbourhood search to address an integrated single machine scheduling and vehicle routing problem. He, Li, and Ram Kumar (2021) developed a branch and price method to solve the IPTSP with integration of the 3D printing and JIT delivery. The objective is to minimise the weighted sum of total transportation time and delivery times of customers. To reduce the whole cost of the supply chain in the make-to-order business environment, Long, Pardalos, and Li (2021) studied a multi-objective integrated production scheduling and vehicle routing problem (IPSVRP) with inventory holding, delivery, and tardiness costs. To obtain a set of diverse non-dominated solutions in the Pareto-optimal front of the problem, they designed a level-based multi-objective particle swarm optimiser based on the derived structural properties. Liu, Guo, and Zhang (2021) investigated an integrated production and outbound distribution scheduling problem with flexible vehicle departure time in a time-sensitive make-to-order supply chain. They proposed a hybrid multi-level optimisation framework by decomposing the problem into three subproblems, including parallel machines scheduling, vehicle assignment and distribution scheduling. Feng et al. (2021) formulated an MILP for a novel crowdsource-enabled IPTSP, and they develop a GA and a lower bound based on problem properties to solve the problem. The majority of existing IPTSP studies have restricted their focus on simple machine configurations and transportation environments, such as a single plant consisting of a single machine (Hall and Potts 2003;Chang and Lee 2004;Low et al. 2014;Li et al. 2016;Lacomme et al. 2018;Long, Pardalos, and Li 2021;He, Li, and Ram Kumar 2021), identical parallel machines (Chang and Lee 2004), or parallel machines (Chen and Vairaktarakis 2005;Guo et al. 2017;Joo and Kim 2017;Liu, Guo, and Zhang 2021) and a single transportation mode with homogeneous vehicles (Viergutz and Knust 2014). On more complex production environments involving multiple operations, Ramezanian, Mohammadi, and Cheraghalikhani (2017) examined a flow shop production environment and compared two delivery methods (direct shipping and routing) to minimise the sum of the production cost and delivery cost. Wang et al. (2020) considered a problem in a hybrid flow shop system with a single capacitated vehicle under the objective of minimising the delivery time to the last customer in the last run. In more recently, Yağmur and Kesen (2021) presented the multi-trip heterogeneous vehicle routing problem under a permutation flow shop production environment and proposed a simulated annealing and a memetic algorithm. Yağmur and Kesen (2022) also studied an IPTSP in which jobs are finished in a job shop environment at first and then delivered by a heterogeneous fleet of vehicles. To solve it, they developed an augmented ε-Constraint method for small-size instances, and two metaheuristics based on local search and non-dominated sorting GA for practical sized instances. However, research on IPTSPs with one or more complicated realistic features, e.g. time-dependent transportation time and an order with multiple jobs, has received relatively little attention to date.
On the aspect of the PSCSP, Li et al. (2010) formulated a precast production planning model with considering installation schedule on-site and production constraints in the factory, which include mould and inventory to minimise the total production cost. They proposed a GA and a branch and bound method for the problem. Maximum of the utilisation of moulds has also been considered in some models. For example, Khalili and Chua (2014) included the manufacturing cost, operation cost, and replacement cost of moulds in their precast scheduling objective, and formulated a mixed integer linear programming model to minimise the production cost. de Athayde Prata, Pitombeira-Neto, and de Moraes Sales (2015) considered the constraints of the available capacity of casting beam moulds and formulated an integer linear programming to minimise production losses on orders. Wang and Hu (2017) integrated the mould manufacturing, PCs storing, and transportation processes to modify the traditional production scheduling model from the perspective of the whole PCs supply chain. Two case studies were conducted to test the validity of the proposed scheduling model based on the GA. Kong et al. (2018) discussed an integrated precast delivery and assembly scheduling problem with Just-in-Time strategy and time-dependent transportation time. They proposed a polynomial algorithm for optimal batch delivery to minimise the sum of earliness/tardiness, resource waste, and emission penalties with time constraints. However, the PCs production stage and mould manufacturing stage are not considered in their work, which might have a great effect on operational efficiency of precast supply chain system. Zhang and Yu (2020) studied the dynamic transportation problem based on the Just-in-Time strategy and proposed a particle swarm optimization to deal with the problem. However, the components production stage and time-dependent transportation times are also not considered in their work. From above we can also see that, most of the existing methods for the PPSP or the PSCSP are metaheuristics, e.g. GA, due to the complexity of the problem. These metaheuristics are computationally efficient. However, it is difficult to evaluate the quality of solutions. Thus, in this paper, we focus on developing efficient matheuristic, metaheuristics, and lower bound methods for the PSCSP_TDT |TWT.
The PSCSP_TDT |TWT is often met in real precast production environments. For example, in a PC manufacturing enterprise in Shanghai, China, once a set of orders arrives at the beginning of the production period, the manufacturer has to coordinate with mould manufacturing companies and logistic companies. To meet due dates and improve service level, job sequence in each order and order sequence should be optimised jointly. Although the PSCSP_TDT |TWT is crucial for increasing the revenue and improving customers satisfaction, to our best knowledge, it has seldom been addressed thus far. In this paper, we focus on the integrated order scheduling and job scheduling in a precast supply chain environment with considering multi-stage production process and time-dependent transportation (TDT). To meet due dates from different customers, our objective is to minimise the total weighted tardiness cost. To deal with this problem, some structural properties are explored based on in-depth analysis on the production and transportation data from a real precast supply chain system. Subsequently, we formulate an approximate model for the PSCSP_TDT |TWT. In the model, the seven production stages are simplified as one differentiation stage which greatly reduce the complexity of the original problem. Based on the approximate model, we propose a dynamic programming heuristic for small instances and four hybrid iterated greedy algorithms for medium-or large instances. We also develop a Lagrangian relaxation bound for evaluating the performances of the AMHIGs.

Problem description and formulation
Prefabricated components (PCs) usually include prefabricated stairs, prefabricated columns, prefabricated balconies, prefabricated force beams, etc. In each production stage, different PCs may have different processing time. To deliver orders in time, it is necessary to make a good schedule for a precast supply chain system. Thus, in this paper, we consider the precast supply chain scheduling problem with time-dependent transportation time and total weighted tardiness criteria (PSCSP_TDT |TWT). The problem can be described as follows. In a real precast supply chain environment, I customer orders from the set I are to be scheduled in nine stages including mould manufacturing (M 1 ), mould assembly (M 2 ), steel embedded (M 3 ), concrete casting (M 4 ), steam curing (M 5 ), mould stripping (M 6 ), repair & finishing (M 7 ), storage (M 8 ), and transportation (M 9 ) (Wang and Hu 2017). We use the pair (i, j) to refer to job j in order i. We let I denote the number of orders, N the number of jobs, and N i the number of jobs belongs to order i. Thus N 1 + N 2 + · · · + N I = N. Each order i consists of a set J i of jobs (i, 1), . . . , (i, j), . . . , (i, N i ) from the same customer. It should be noticed that the job sequence in each order needs to be scheduled. Then we have the set J = J 1 ∪ · · · ∪ J i ∪ · · · ∪ J I of all jobs. The nine stages are classified into two situations including: (1) Sequential and (2) Parallel. In sequential stages, at most one job can be processed at a time. Whereas, in the parallel stages (M 5 , M 8 , and M 9 ), the numbers of machines are large enough to process all available jobs simultaneously at any time slot. For easy management, all jobs of the same order are processed in a group without interleaving by jobs belongs to other orders. The processing time of job (i, j) at stage m ∈ M \ {9} is denoted by p m, (i,j) . Only if all jobs in an order are completed at stage 8 the order can be delivered to its customer. Transportation starts at an exact integer hour and transportation time is a step function of the delivery time. Associated with each order is a weight w i and a due date d i .
The objective is to determine job sequences in each order and an order sequence to minimise the total weighted tardiness. The following variables are defined for a given schedule π: C J m,(i,j) (π ): the completion time of job (i, j) at stage m; C O m,i (π ) = max j∈{1,...,N i } C J m,(i,j) (π ), the completion time of order i at stage m. When there are no ambiguity, we simplify C J m,(i,j) (π ) and C O m,i (π ) as C J m, (i,j) and C O m,i , respectively. Let T i denote the tardiness of order O i , and it can be calculated by max{0, C O 9,i − d i }. In this paper, i∈I w i T i is simplified as TWT, and optimal value is denoted by TWT * . Several assumptions are done as follows.
• All machines are available at time 0 and all orders are released at time 0 • All the information about all the orders are known in advance • The buffer size between two stages is infinite • Setup times and loading times are neglected • An order cannot be loaded and transported until all jobs in it are completed at the production stage • All customers locate at the same district. Thus they have the same transportation time which is only varied with the departure time • Traffic jam is regarded as the same in every day For easy understanding, the pertinent notations for the PSCSP_TDT|TWT are illustrated as follows: The order scheduled in the i th position for a permutation schedule The completion time of job (i, j) at stage m, The completion time of order i which is the time when all jobs of order i is completed at A solution for the PSCSP_TDT|TWT, where π = (π 1 , . . . , π i , . . . , π I ). π i is the job sequence in order [i ].
) TWT(π ) Total weighted tardiness of π for the PSCSP_ TDT|TWT S J m, [i ,j ] The start time of the job in position j of the i th order of a schedule S O m, [i ] The start time of the order in position i of a schedule C J m, [i ,j ] The completion time of the job in the j th position of the i th order of a schedule C O m, [i ] The completion time of the order in position i of a schedule Aimed at calculating the tardiness of an order, we need to determine its completion time. Given a permutation schedule (π 1 , . . . , π i , . . . , π I ) for the PSCSP_TDT |TWT, and [i , j ] to denote the order scheduled in the i th position, and the job scheduled in the j th position of order [i ], respectively. Let p m, [i ,j ] denote the processing time of job [i , j ] at stage m. Each job completes as early as possible at all stages. Thus S J 1,[1,1] = 0. Similar as permutation flow shops (Pinedo 2016 [i ,j ] ) can be computed easily through a set of recursive equations: Equation (1) gives the machine initialisation constraints for the job in the first position of the first order (job [1, 1]), and Equations (2)-(4) represent the jobs' initialisation constraints at stage 1. Both Equations (5) and (6) give the job availability constraints at all stages except stages 5, 8, and 9. Specially, Equation (5) also links the job in the first position of an order and the job in the last position of the immediately preceding order at the same stage. Whereas Equation (6) links two adjacent jobs in an order at the same stage. Only the machine availability constraints are given in Equation (7) because there are enough number of parallel machines at stages 5 and 8 for processing all jobs simultaneously. Equation (8)  ) should be determined. In this study, the transportation time is varied with the departure time. This assumption arises from the real crowded urban transportation environments in which the transportation times are largely affected by many factors, e.g.traffic density (Dabia et al. 2013;Kong et al. 2018). To meet precast transportation situations as real as possible, this study considers the time-dependent transportation time under different traffic flows. According to traffic levels (light, normal, heavy, and major traffic jams), the transportation time between single manufacturer and single customer is predicted. Considering the fixed manufacturers and customers, according to the statistical results of departure time in urban traffic flow, the transportation time (the processing time of the last stage) in the precast supply chain system is formulated as Equation (11) which is a step function to reflect the traffic congestion degree at different departure times.
where K is the total number of time intervals. Let transportation time varies with the start time. During the first 4 hours or the last 6 hours of a day, light traffic jams will be encountered. In this case, 2.5 hours on average will be costed for transportation. Major traffic jams will happen during 7:00 AM to 9:00 AM, which means a long transportation time (5 hours on average). If the orders are started to be transported during 5:00 AM to 6:00 AM or 9:00 AM to 16:00 PM, the traffic jams will be in normal level and it will take 3.2 hours on average for transportation. In other time slots, heavy traffic jams will be encountered and it will take 4.5 hours for transportation.
For an order [i ], once the start time of transportation (S O 9,[i ] ) is given, the transportation time (tr(S O 9,[i ] )) can be determined by Equation (11). Thus it is a key problem to choose the optimal start time S O 9,[i ] to minimise the sum of the waiting time and the transportation time for a given order sequence. Note that transportation stage is a parallel stage in which orders can be transported in parallel. The optimal departure time for an order is not necessarily the earliest departure time because of the time-dependent transportation times. For easily understanding, a simple example is given as follows. Assume that transportation time is shown in Figure 1. Suppose that the completion time of order [i ] at stage 8 is 6:30 AM. If order [i ] is departed at 7:00 AM, then the transportation time will be 5 hours. However, if the order is departed 1 hour later (8:00 AM), the transportation time will be 3.2 hours. Thus 5.0 − (1.0 + 3.2) = 0.8 hour will be saved. Let tr max − tr min . Based on the fact that we cannot gain more than tr max hours of travel time even by choosing the moment when the traffic density is minimal, there is no use of waiting more than tr max after the moment C O 8, [i ] when the transportation may start. We define S 9, [i ] where S O 9,[i ] is the start time when order [i ] is departed from the factory, and it is nonnegative integer hour.
With Equations (1)-(12), for a given schedule π , the completion time, i.e. delivery time of an order (C O 9,[i ] ) can be obtained. Thus the total weighted tardiness (TWT) can be further calculated by

Approximate problem and some properties
In this section, the real data from a PC manufacturer will be analysed at first, and then the PSCSP_TDT |TWT is approximated by a three-stage supply chain scheduling problem named as APSCSP_TDT |TWT. Finally, some important properties and propositions are explored for developing algorithms.

Data analysis and approximate PSCSP_TDT |TWT
In this section, we approximate the PSCSP_TDT |TWT as a three-stage scheduling problem based on observing and analysing the real data from a precast manufacturer. Table 1 lists the real production and transportation data of 10 types of PCs from a PCs manufacturer in Shanghai, China (Wang and Hu 2017). Provided that we have a set J of jobs which belong to these types of PCs. From Table 1, we observe that, min (u,v)∈J p 1,(u,v) = 4.0, and all jobs (i, j) ∈ J satisfy the conditions min (u,v)∈J p 1,(u,v) ≥ m∈{2,3,4} p m,(i,j) except the jobs with type 9. All jobs satisfying p 5,(i,j) ≥ 8(p 6,(i,j) + p 7,(i,j) )/3 implies that the processing time of stage 5 greatly dominates the sum of processing times of stage 6 and stage 7. Based on the above observations and analysis, we obtain the following properties.
Proposition 4.1: Given a permutation schedule π for the PSCSP_TDT |TWT, 3,4,5} p m,(i,j) ; and with Condition 2 and Equations (1)-(8), we have C J 8,(i,j) = C J 5,(i,j) + m∈{6,7,8} p m, (i,j) , ∀(i, j). Thus if both conditions (1) and (2) are satisfied, then we have Based on the above analysis and Propositions 4.1, we approximately simply the nine-stage supply chain system Table 1. Production data of processing time from a real precast plant in Shanghai, China (Wang and Hu 2017). into a three-stage supply chain system. In the approximate system, the first stage is mould manufacturing, which is actually a bottleneck stage. The second stage is a differentiation stage which is simplified by combining from the second to eighth stages of the original supply chain model. The last stage is a distribution stage which collects all jobs from the same order to deliver. Let p A k, (i,j) and C J,A k,(i,j) be the processing time and completion time of job (i, j) at stage k ∈ {1, 2, 3} for the approximate system, respectively. For easily understanding, given a simple example with three orders O 1 , O 2 and O 3 as depicted in Figure 2. There are three orders, each of which includes three or four jobs, to be processed. Each job must go through three stages. The first stage is the same as that in the PSCSP_TDT |TWT, all jobs are processed in sequential with a processing time p A 1,(i,j) = p 1,(i,j) . In the second stage, each job is processed in their dedicated machine with a processing time p A 2,(i,j) , where p A 2,(i,j) := 8 m=2 p m, (i,j) . Thus the second stage is in parallel. The last stage is transportation stage which is the same as that in PSCSP_TDT |TWT. Therefore, the PSCSP_TDT |TWT is approximately transformed into a three-stage precast supply chain scheduling problem. For simplicity, we denoted the approximate problem by APSCSP_TDT |TWT. Notice that, if all the conditions in Proposition 4.1 are satisfied for any schedule, the approximate problem will be equivalent to the considered original one. Let TWT A denote the total weighted tardiness for the approximate problem, then the approximate problem is to find the optimal job sequence in each order and the optimal order sequence to minimise the total weighted tardiness. At the transportation stage, the approximate problem has the same step transportation time function as Equation (11).
Given a permutation schedule π = (π 1 , . . . , π i , . . . , π I ) for the APSCSP_TDT |TWT, let p A k, [i ,j ] Equations (14)- (16) give the jobs' initialisation constraints at stage 1. Equation (17) where S O,A 3,[i ] is the start time when order [i ] is departed from the factory, and it is nonnegative integer hour. Thus the total weighted tardiness for the APSCSP_TDT |TWT (TWT A ) can be further calculated by Notice that, if the differentiation stage and the transportation stage are not considered, the APSCSP_TDT |TWT will be reduced to a single machine scheduling problem 1|| w i T i , which is strongly N P-hard (Pinedo 2016). Thus APSCSP_TDT |TWT is also strongly N P-hard.

Some structural properties of the schedules
In this section, some important properties, which can be applied to develop the algorithms, are presented. for the APSCSP_TDT |TWT, let S O,A 1,[i ] denote the start time of order i, which is the start time of the job in the first position in order i. Then we have the following propositions.

Proof: Given a permutation job sequence
Notice that for any sequence π we obtain above result. Set the job in last is an optimal solution for the APSCSP_TDT |TWT with a fixed order sequence σ .
Proof: Because order sequence is fixed and both the differentiation stage and transportation stage are in parallel, any operations in an order, e.g.swap operation, and insert operation, will not affect the completion time of other orders. We use TWT A to denote the total weighted tardiness for the APSCSP_TDT |TWT. * [i ] are the minimum of completion time and the minimum of tardiness penalty of order Thus Propositions 4.2 and 4.3 illustrate that, given a fixed order sequence APSCSP_TDT |TWT, we can easily obtain an optimal job sequence in each order by scheduling the job, which has the minimum processing time at stage 2, into the last position of the order. Thus, for the APSCSP_TDT |TWT, we can focus on dealing with the order scheduling problem and do not need to pay much time to solve the job scheduling problem in each order. This will substantially reduce the computational effort.

Proposition 4.4:
Given an optimal permutation schedule π * for the APSCSP_TDT |TWT, let the optimal TWT value be TWT A (π * ), then TWT A (π * ) is a lower bound for the corresponding PSCSP_TDT |TWT.
Proof: Suppose that the optimal solution for the PSCSP_ TDT |TWT is π * 0 and the optimal TWT value is and T i are the tardinesses of order i for the APSCSP_TDT |TWT and PSCSP_TDT |TWT, respectively. Therefore, we have Because π * 0 is also a feasible solution for the APSCSP_ TDT |TWT, we have With Equations (24) and (25), we have TWT A (π * ) ≤ TWT A (π * 0 ) ≤ TWT(π * 0 ). Thus TWT A (π * ) is a lower bound for the corresponding PSCSP_TDT |TWT.
Proposition 4.4 illustrates that, if we can obtain the optimal value for an APSCSP_TDT |TWT, then we also derive a lower bound for the corresponding PSCSP_TDT |TWT. Thus the optimal value for the APSCSP_TDT |TWT can be applied to evaluate the performance of a schedule for the PSCSP_TDT |TWT.

Solution methodology
From Proposition 4.4, we know that the optimal TWT value for the APSCSP_TDT |TWT is a lower bound for the PSCSP_TDT |TWT. Dynamic programming (DP) is a general exact optimisation method for solving scheduling problems that can be decomposed into subproblems, each involving a subset of the sequence decisions, in such a way that the optimality principle holds (Bellman 1966;Bellman and Dreyfus 2015;Blazewicz et al. 2019;Baker and Trietsch 2019). Therefore, in this study, with the structural properties of the optimal solutions, we design an approximate model-based hybrid dynamic programming and heuristic in two-fold motivations. The first one is that, according to Proposition 4.4, we can obtain an effective lower bound for small or middle size instances for the PSCSP_TDT |TWT by applying the DP to APSCSP_TDT |TWT. Whereas the second one is to obtain a better sequence by DP and then improve the sequence to achieve an optimal or near optimal solution for the PSCSP_TDT |TWT. However, it is difficult to obtain an optimal solution for the APSCSP_TDT |TWT within a reasonable CPU time by using exact methods, e.g. DP, especially for large-scale instances because of curse of dimensionality (Bertsekas 2017). Therefore, by integrating the properties of the optimal solution of the problem, dispatching rule and constructive heuristics, different local search methods, mechanism of destruction and reconstruction, and a correction heuristic, we propose four approximate model-based hybrid iterated greedy (AMHIG) algorithms which were developed based on the APSCSP_TDT |TWT. For medium-or large instances, to evaluate the performances of four AMHIGs, we formulate the APSCSP_TDT |TWT as a timeindexed integer programming model and obtain a lower bound via Lagrangian relaxation by relaxing the capacity constraints.

Approximate model-based hybrid dynamic programming and heuristic (AMHDPH) for the PSCSP_TDT |TWT
In this section, we proposed an approximate modelbased hybrid dynamic programming and heuristic (AMHDPH). The main idea is that, at first, obtain an optimal solution π * and its total weighted tardiness for the APSCSP_TDT |TWT (TWT A (π * )) by using a dynamic programming method. Then we calculate the total weighted tardiness for the PSCSP_TDT |TWT (TWT(π * )). If TWT(π * ) = TWT A (π * ), then π * is the optimal solution for the PSCSP_TDT |TWT; else, we run a correction heuristic (CH) to further improve the incumbent. Finally, we output the best solution π * , TWT(π * ) and the gap between TWT(π * ) and TWT A (π * ). Note that TWT A (π * ) is a lower bound for the PSCSP_TDT |TWT according to Proposition 4.4.
To apply dynamic programming for the APSCSP_TDT |TWT, let I sub denote some subset of the orders, and let C A 3,i denote the completion time of order i, which can be calculated by using Equations (15)-(22). For convenience, we use (I sub \ {i}) to denote the set I sub with the order i removed. Suppose that an order sequence has been constructed in which orders in set I set precede all other orders. Let V(I sub ) denote the minimum cost for the subproblem consisting of the orders in set I sub . Given the order i comes last, the value of V(I sub ) is the sum of two terms, the cost incurred by order i and the optimal cost incurred by the remaining orders. This latter term which can be written as V(I sub \ {i}) is the optimal cost obtained by the subproblem involving only the orders in set I sub \ {i}. If we compare all possible orders i that could come last in set I sub and select the best one, we shall find the minimum cost for the set I sub . That is where and ∅ denotes the empty set, and is the weighted tardiness of order i. Finally, because the cost function V is defined on subsets of orders, the minimum total cost can be written as V(I) which is computed as follows: At each stage, the function V(I sub ) measures the total cost contributed by the orders in set I sub , when set I sub occurs at the beginning of the schedules and is sequenced optimally. The recursion relation (26) shows that aimed at calculating the value of V for any particular subset of size k, we first have to know the value of V for subset of size k−1. Thus the DP procedure for the APSCSP_TDT |TWT begins with the value of V for a subset of size zero, from Equation (27). Then, the value of V for all subsets of size 1 can be calculated by using Equations (27) and (28), then the value of V for all subsets of size 2, and so on. In this manner, the procedure considers ever larger sets I, ultimately using Equation (29) to determine which order should be scheduled last. The optimal value of TWT A is V(I). If we keep track of where minima in Equation (26) occurs at each stage, then, after finding V(I), the optimal order sequence π * for the APSCSP_TDT |TWT can be reconstructed. After obtaining an optimal solution π * for the APSCSP_TDT |TWT, we will compare the values of TWT(π * ) and TWT A (π * ). If the values are equal, the optimal solution for the PSCSP_TDT |TWT is achieved. Otherwise, we use a heuristic to see whether the incumbent can be improved. For easy understanding, consider a simple example with 3 orders and 12 jobs, as shown in Figure 3. Note that at Stage 1 in Figure 3, the first two jobs belong to Order 1, the 3rd to 6th jobs belong to Order 2, and the 7th to 12th jobs belong to Order 3. From Figure 3(a), we observe that, in the Gantt chart of the schedule obtained after DP stage of AMHDPH, a job belongs to Order 3 is delayed by its preceding job at the same stage (Stage 6 or 7) as pointed in Figure 3(a), which may increase the TWT value. If two jobs of the same order are interchanged as shown in Figure 3(b), the job delay will be eliminated and the TWT value will be reduced. This case illustrates that the error of APSCSP_TDT |TWT to PSCSP_TDT |TWT might be caused by delayed operations in Stage 6 or 7. This highlights us to design a correction heuristic (CH) by rescheduling the jobs in the first N i − 1 positions in each order. The CH is a very simple and effective heuristic (as shown in the following experimental results in Section 6), and its main idea is that, to reduce delaying time of the job in last position of an order, for each order, we reschedule the jobs in the first N i − 1 positions by decreasing values of p 6,(i,j) + p 7,(i,j) . The CH procedure is summarised in Algorithm 1. Note that (π * ) = 0 implies that the solution π * cannot be improved for a fixed order sequence. Whereas (π * c ) − (π * ) < 0 implies that the solution π * is improved by the CH.
The computational complexity of the AMHDPH can be analysed as follows. The value of V(I sub ) has to be determined for all subsets that contains k orders. There are I!/k!(I − k)! subsets. Thus the total number of evaluations that have to be done are I k=1 I!/k!(I − k)! = O(2 I ), which means that the computational complexity of the AMHDPH increases exponentially with the increase in the number of orders (I). Therefore, to deal with large-size instances, we proposed an effective iterative greedy-based framework in the following section.

Approximate model-based hybrid iterated greedy method (AMHIG) for the PSCSP_TDT |TWT
Iterated greedy (IG) is a local search-based method that iterates through applications of construction heuristics using the repeated execution of two main phases, the partial destruction of a complete candidate solution and a subsequent reconstruction of a complete candidate solution. IG has many success applications to complex scheduling problems Stützle 2007, 2008;Minella, Ruiz, and Ciavotta 2011; and Naderi 2019; Fernandez-Viagas and Framinan 2019; Riahi, Chiong, and Zhang 2020). In this paper, we proposed an approximate model-based hybrid iterated greedy algorithm (AMHIG) by integrating the structural properties of schedules, constructive heuristics and an iterated greedy algorithm. The procedure of AMHIG is shown in Algorithm 2. Note that in Algorithm 2, π * is the best solution obtained after IG stage, and its TWT values for the APSCSP_TDT|TWT and PSCSP_TDT |TWT are TWT A (π * ) and TWT(π * ), respectively. π * c is the solution obtained after CH stage in AMHIG, and its TWT value for the PSCSP_TDT |TWT is TWT(π * c ). (π * ) and (π * c ) are defined by TWT(π * ) − TWT A (π * ), and TWT(π * c ) − TWT A (π * ), respectively. The main idea of AMHIG is that, at first, the optimal job sequence in each order for the APSCSP_TDT |TWT is determined by using Proposition 4.2, as shown in lines 2-5 in Algorithm 2. Subsequently, a constructive heuristic named HWGI, which will be introduced in Section 5.2.1, is executed to obtain a better initial solution for the IG in line 6. Then, we run an iterated greedy algorithm (lines 8-26 ) to obtain a better order sequence for the APSCSP_TDT |TWT and compute the cost of the sequence for the PSCSP_TDT |TWT. Finally, as shown in lines 27, we try to further improve solution quality by calling a CH which is the same as that in AMHDPH. Note that in AMHIGs, once the optimal job sequence in each order is determined, it will remain unchanged until the correction heuristic stage starts. This will save much computation time. Note that p A 2,(i,j) is the processing time of job (i, j) at stage 2 (Differentiation stage) for the APSCSP_TDT |TWT. In the following sections, some key components of our AMHIGs, such as solution representation, evaluation and initialisation, constructive heuristic, different local search methods, mechanism of destruction and construction, and correction heuristic, will be presented.

Solution representation, evaluation and initialisation
In our AMHIG, a solution is represented as permutation sequence of I orders because for the approximate Algorithm 2: AMHIG procedure for the PSCSP_TDT |TWT Input: Instance data, algorithm parameters: maximum runtime t max Output: π * , TWT A (π * ), TWT(π * ), π * c , TWT(π * c ), (π * ) and (π * c ). 1: Set t 0 ← cputime() 2: for all i = 1 to I do 3: Randomly generate a permutation job sequence for order i 4: In order i, find the job with the smallest p A 2, (i,j) value and interchange it with the job in the last position, and then obtain a job sequence π i 5: end for 6: π ← HWGI 7: π ← LocalSearch(π ) 8: π * ← π 9: Set t ← cputime() − t 0 10: while t ≤ t max do 11: π ← π 12: Set t ← cputime() − t 0 ; 26: end while 27: π * c ← CH(π * ) 28: Calculate TWT A (π * ), TWT(π * ) and TWT(π * c ) 29: three-stage scheduling model, the optimal job sequence in each order for the APSCSP_TDT |TWT can be easily obtained by using Property 4.2. Highlighted by Kanet and Li (2004), we generalise the weighted modified due date (WMDD) rule for the 1|w j T j problem to the APSCSP_TDT |TWT. That is to generate an order sequence in a nondecreasing order of the quantity max{d i − t, C O,A 3,i − t}/w i , where t is the start time of π ← (π , O k ) 8: end for 9: Let O f be the order in O 10: π ← (π , O f ) 11: for all i = 1 to I − 1 do 12: π ← π 13: Obtain a new sequence π by interchanging the orders in position i and i + 1; 14: if TWT A (π ) < TWT A (π 0 ) then 15: 17: end if 18: end for 19: Let π be the best permutation between (π 0 1 ,π 0 2 ) and (π 0 2 ,π 0 1 ) 20: for all i = 3 to I do 21: Insert order π 0 i in the position of π yielding the lowest TWT A value 22: end for 23: Calculate TWT A (π ) 24: Let π greedy = ∅ 25: for all i = N to 1 do 26: From order set choose the order π best yielding the lowest TWT A value when it finishes last and add it to the ith position of π greedy 27: = \ {π best } 28: end for 29: Calculate TWT A (π greedy ); 30: if TWT A (π greedy )) < TWT A (π ) then 31: π ← π greedy 32: TWT A (π ) ← TWT A (π greedy ) 33: end if 34: return π, TWT A (π ) order i, and C O,A 3,i can be obtained by Equations (15)-(22). Dispatching by this rule involves sorting the orders and choosing the order with the lowest WMDD as the next order. The order sequence is dynamic because the dispatching criterion depends on the start time t.
We develop an effective heuristic named HWGI by integration of a WMDD rule, an NEH-like heuristic (Nawaz, Emory Enscore Jr, and Ham 1983) and a greedy heuristic, as shown in Algorithm 3. At first, the procedure start a WMDD to generate an initial solution (lines 1-10), then improve the solution by using a correction algorithm in which two adjacent orders are interchanged from the first position to the last position (lines 11-18). Subsequently, we use an NEH-like insertion method to further improve the solution (lines 19-23). Above process can be seen as a WMDD_Correction_Insertion procedure, denoted as WMDD_C_I. Finally, we applied the greedy heuristic to construct an another order sequence from the last position to the first position (lines 24-29). We compare the TWT values of the two sequences obtained by WMDD_C_I and Greedy, and then choose the better one. A metaheuristic needs at least one start point which may be randomly generated. Because a random solution is expected to be very low quality, it is usually common to initialise metaheursitics with good known constructive heuristics to improve search efficiency (Framinan, Leisten, and Ruiz García 2014). Thus, in this paper, we use the solution obtained by the HWGI as the initial solution of the AMHIG. Note that, the job sequence in each order is generated according to Proposition 4.2, as shown in lines 2-5 in Algorithm 2.

Local search
Local search methods are generally used to improve solution quality in IG Stützle 2007, 2008;Ruiz, Pan, and Naderi 2019;Riahi, Chiong, and Zhang 2020). Swap operator is to select two positions in a sequence and interchange the two jobs in the two positions. This operator is one of most often used operators in local search-based method (Pinedo 2016;Martí, Pardalos, and Resende 2018;Lourenço, Martin, and Stützle 2019). Insert operator is to extract one order from a sequence and insert it in another position. It is also one of the most often used operators in local search-based method Stützle 2007, 2008;Ruiz, Pan, and Naderi 2019;Li et al. 2019). In this paper, based on swap and insert operators, three different local search methods (LS 1 , LS 2 and LS 3 ) are developed for AMHIG to improve search efficiency.
LS 1 exchanges orders between two positions in an order sequence. For each position k = 1, 2, . . . , I − 1, the order in position k is swapped with its following jobs to find the sequence yielding the lowest TWT. If the objective value is improved, update the current solution and then set k: = k + 1. Repeat the same operations as previous steps.
LS 2 , which is extended from Ruiz and Stützle (2007), search for better solutions by extracting an order from one position and inserting it into another position. In LS 2 , an order is randomly extracted from the sequence and find the position yielding the lowest TWT A . If the best TWT A is better than the starting TWT A , the order is relocated and the search starts again from the beginning. Otherwise, the order is reinserted back into its original position and the search continues. The procedure iterates until all orders from the sequence have been tested.
LS 3 is a variable neighbourhood search (VNS) which incorporates parts of the previous LS 1 and LS 2 . The main idea of LS 3 is to employ LS 1 and LS 2 alternatively until stopping criterion is met. The procedure starts from LS 1 . If the starting solution is improved, the search will be continued in LS 1 . Otherwise, LS 2 will be run. If the solution is improved by LS 2 , the local search goes back to LS 1 . Otherwise, the procedure will be terminated.
Note that the pseudo codes of LS 1 , LS 2 and LS 3 can be found in Algorithms S1, S2 and S3 of Supplemental Data.

Destruction and reconstruction (DR) in AMHIG
The main principle of iterated greedy is to iterate over construction methods by first generating a complete candidate solution and then cycling through a main loop that consists of two main steps. In the first step, which is called the destruction step, some solution components are removed from the current complete candidate solution π to result in some intermediate partial candidate solution π D . In the next step, which is called reconstruction step, starting from π D , a constructive heuristic is used to generate a new complete solution π 0 . Then an acceptance test is conducted on the two solutions, π or π 0 , and one of them are chosen for accepted. The next destruction step applies. Given a permutation solution π for the APSCSP_TDT|TWT, in this paper, we extend the standard destruction and reconstruction (SDR) steps in Stützle (2007, 2008) to our AMHIG as follows. (1) Destruction phase: Randomly choose d orders from π . Delete them from π and add them to π R in the order in which they are chosen, where π R is the sequence of the d removed orders. (2) Reconstruction phase: Sequentially reinsert the orders of π R into π D until a complete solution is obtained, where π D is the partial sequence of π after removing the sampled orders. During the reconstruction procedure, the best (partial) solution, π , obtained by inserting the orders of π R into all possible positions of the current subsequence, is recorded.
For tardiness criterion, Riahi, Chiong, and Zhang (2020) proposed a tardiness-guided destruction and reconstruction (TG_DR) method in their IG algorithm to deal with the no-idle permutation flow shop scheduling problem. Computational results verified the effectiveness of the TG_DR method. Highlighted by Riahi, Chiong, and Zhang (2020), in the AMHIG_VNS, in addition to the SDR method presented as above, we also apply a weighted tardiness-guided destruction and reconstruction (WTG_DR) method. The basic idea of the WTG_DR is that, in destruction stage, instead of random selection, orders with the highest weighed tardiness values are selected. First, orders are sorted by decreasing values of their weighted tardiness values and saved in list θ. Next, the first d/2 orders of list θ are selected and inserted into π d . The remaining orders are then chosen at random to guarantee a different solution at each iteration. In the reconstruction stage, both the standard insert move and swap move are used. With the swap move, two orders are chosen in either a random or greedy manner from the current sequence and then, their positions are exchanged. The insert and swap moves work as follows. First, a number between 0 and 1 is generated randomly. If the number is less than or equal to 0.5, then the insert move is used. Otherwise, the swap move is applied. For more details of TG_DR, please refer to Riahi, Chiong, and Zhang (2020).
Based on different local search methods and DR methods, we proposed four variants of AMHIG, named AMHIG_Swap, AMHIG_Insert, AMHIG_VNS, and AMHIG_VNS_TG. LS 1 , LS 2 , and LS 3 are used in the local search steps of AMHIG_Swap, AMHIG_Insert and AMHIG_VNS, respectively. And all the AMHIGs except AMHIG_VNS_TG use SDR method. For AMHIG_VNS_ TG, LS 3 is used as local search and TG_DR method is applied in DR steps.

Acceptance criterion, stopping criterion and correction heuristic (CH)
The acceptance criterion has a strong influence on the diversification and intensification behaviour of an IG algorithm. An appropriate choice may be crucial to the performance of an IG algorithm. In this paper, we adopt the Metropolis criterion which is a popular choice when accepting worse solutions Stützle 2007, 2008;Ruiz, Pan, and Naderi 2019). That is, if a new solution is better or equal the incumbent, it is accepted; otherwise it is accepted with a probability exp{(TWT A (π ) − TWT A (π ))/T}, where T is called temperature. It can be calculated as follows: where T para is a parameter for controlling temperature T, and p 9,(i,j) := tr min . The stopping criterion is set as maximum allowable CPU runtime. At the last stage of AMHIG, we also use a correction heuristic (CH), which is the same as Algorithm 1 in AMHDPH, to improve solution quality for the PSCSP_TDT |TWT.

Computational experiments
In this section, computational experiments are designed and conducted to verify the effectiveness of the proposed AMHDPH, and the four AMHIG algorithms. In Section 6.1, computational experiments are designed according to the characteristics of the problem, and then experimental results are presented and analysed in Section 6.2. Finally, managerial insights and limitations are presented in Section 6.3.

Experimental design
In this paper, all the proposed algorithms have been  Table 1. For each job in each instance, we randomly choose one of the 10 types of prefabricated components listed in Table 1 as the product type of the job. Then the processing time of the job can be determined from the table. However, although the processing times of all test instances are from the same set of real-world data, they may be very different because of different demands of customers. For example, we generate two instances (Ins1 and Ins2), each of which has one order and five jobs. In Ins1, there are two jobs with type 2, two jobs with type 3, one job with type 10. Whereas, in Ins2, there are two jobs with type 3, one job with type 2, one job with type 8, one job with type 9. Four traffic levels are considered (light, normal, heavy, and major traffic jams). The transportation time tr(b i ) from the PCs factory to the construction site is set as Equation (31) which illustrates the distribution of transportation time at different departure time.
where D i = b i /24 and b i is the departure time of order i when it is departed from the factory. The per unit time tardiness penalty of each order is set as w i = j∈{1,...,N i } β (i,j) /10, where β (i,j) denotes per unit time tardiness penalty of job (i, j) which is generated from a uniform distribution on the interval of [2,10]. The due date of each example is generated from a uniform distribution on the interval of [P(1 − τ − R/2), P(1 − τ + R/2)]. Parameters τ and R take values inside the set 0.4, 0.6 and 0.2, 0.6, respectively. P is a good makespan obtained by the proposed AMHIGs in which the objective is replaced by makespan. The parameter τ denotes the tardiness factor, which controls average due-date values and parameter R controls the relative range of due dates. Generating due dates from such kind of interval follows the methods described in Potts and Van Wassenhove (1982) and Vallada, Ruiz, and Minella (2008). For each combination of I × N × τ × R, we generate five instances. Each instance is run 10 times for the proposed AMHIGs. The time limit for the AMHDPH is set as 7200s. Aimed at allowing more computation time as the number of orders increases, the termination criterion of metaheuristics for many scheduling problems are set to a CPU time which is dependent to number of orders (Pan and Ruiz 2012;Ruiz, Pan, and Naderi 2019). Therefore, in this study, the termination criterion used to test all the proposed AMHIGs for the same instance is set to a CPU time limited to 0.6 × I seconds.
It is crucial to tune the parameters because the performance of a metaheuristic is usually sensitive to them. In this paper, each AMHIG algorithm has two parameters: d and T. Taguchi method is an effective design of experiment approach (Montgomery 2013) and has been widely used to parameter tuning of metaheuristics with multiple parameters (Ruiz, Pan, and Naderi 2019). Therefore, in this study, we applied Taguchi method to tune the parameters of four AMHIG algorithms. Because there are only two parameters for each algorithm, we use full factor tests for parameter tuning. The process of parameter tuning is omitted in this paper due to space limitation, and it can be found in Supplemental Data. The parameters for the AMHIG_VNS, AMHIG_VNS_TG, AMHIG_Swap and and AMHIG_Insert were determined as follows.
To evaluate the performances of the proposed algorithms, the optimality gap between the TWT value and the LB of an algorithm for an instance is defined as follows: GAP alg (l) = (TWT alg (l) − LB(l)) × 100/(TWT alg (l) + 10 −10 ), where TWT alg (l) is the TWT value obtained for a given algorithm alg and instance l for the PSCSP_TDT |TWTs. For small size instances with I ≤ 20, the LB is obtained by the DP procedure. Whereas, for the middle or large size instances with I > 20, the LB is obtained by using LR method because DP cannot yield the optimal solution for the PSCSP_TDT |TWT within 7200 s (Polyak 1969;Chen and Luh 2003;Bertsekas 2016). The details for the lower bound obtained by the LR method can be found in Supplemental Data. The constant 10 −10 is added just to avoid dividing by zero. To facilitate statistical analysis, we also define the average relative percentage deviation (ARPD) as follows: (TWT(l) − TWT best (l)) * 100/(TWT best + 10 −10 ), where L is the total number of instances and TWT best (l) is the best TWT value obtained so far for instance l.

Computational results and analysis on the AMHDPH
Computational results of the LB DP , TWT(π * ) and TWT(π * c ) appear in Table 2, where LB DP denotes the lower bound generated by the DP stage in the AMHDPH. Suppose that π * is the optimal sequence after DP stage. Thus we have LB DP = TWT A (π * ). π * c denotes the solution generated by using the correction heuristic on π * . To classify the instances in the data set, we used a three digit notation, e.g. I_20_80_8. The first two numbers in this notation indicate how many orders and jobs are used in a given instance, respectively. Last digit of this notation indicates the index of the instance for each combination. GAP 1 denotes the gap between TWT(π * ) and LB DP , and GAP 2 denotes the gap between TWT(π * c ) and LB DP . As shown in Table 2, for each combination of I and N, we test all algorithms by using 25 instances with different τ and R. For each combination of I × N, τ = 0.4 and R = 0.2 for the first five instances, τ = 0.4 and R = 0.6 for the second quarter instances, τ = 0.6 and R = 0.2 for the third quarter instances, and τ = 0.6 and R = 0.6 for the last five instances.
As seen from Table 2, the lower bound generated by AMHDPH (LB DP ) is reasonable tight. The optimality gap between TWT(π * ) and LB DP (GAP 1 ) is well below 0.5 percent in 39 out of 40 instances, and 32 instances are solved to optimality, i.e. 80% of all small instances.
Whereas the gap between TWT(π * c ) and LB DP (GAP 2 ) is well below 0.5 percent for all 40 instances, and the optimal values are achieved for 36 out of 40 instances. It also can been seen from Table 2, CH improves the solution qualities for 6 out of 8 instances in which the gaps between TWT(π * ) and LB DP are not zeros. Among the 6 instances, optimal solutions are achieved by CH for 4 instances. This verifies the effectiveness of CH in AMHDPH. Note that although optimal solutions are not proved for 4 out of 40 instances, the optimality gap is very small and might be caused by inherent errors generated by the approximate model. Because the solution obtained after the DP stage in the AMHDPH is optimal for the APSCSP_TDT |TWT, all the results as shown in the tables strongly support the high performance of the approximate model up to 20 orders.
To further verify the effectiveness of the AMHDPH, we formulate two mixed integer linear programming (MILP) models for the APSCSP_TDT|TWT and the PSCSP_TDT|TWT, respectively (see Supplemental Data). For brevity, we denote them as MILP_APSCSP and MILP_PSCSP. The two MILP models are solved in Gurobi 9.1.2 coded in Julia+JuMP (Lubin and Dunning 2015;Dunning, Huchette, and Lubin 2017), and run with a time limit of 3600s. The computational results on the 40 instances with 10 orders and 20 orders are listed in Table 3. As shown in Table 3, the optimal results of the MILP_APSCSP are obtained by Gurobi in 17 of 20 small-size instances with 10 orders, which are marked by ' * '. Whereas the MILP_PSCSP can only find optima in three instances. For each solved instance, it costs at least 96.14s to solve the MILP models via Gurobi. Whereas AMHDPH only needs at most 0.05 s for instances with 10 orders. For all instances with 20 orders, no MILP model can be solved optimality via Gurobi within 1 hour. However, the AMHDPH finds optima in all instances with 20 orders within at most 60s. It should be noticed that, although some instances, e.g. I_20_80_1, have obtain optimal solutions actually, the optimality gaps generated by Gurobi (GAP Gurobi ) are still very large (100% for I_20_80_1).
The CPU time will increase dramatically as order number increases. AMHDPH cannot obtain optimal solutions within 7200 s at order number I ≥ 30. Thus, to deal with medium-or large instances, we develop effective metaheuristics by integration of problem-specific knowledge, constructive heuristics, local search methods and mechanism of destruction and construction and CH. To evaluate their performances, we proposed a Lagrangian relaxation bound based on the separable structure of the problem. The computational results and analysis are displayed in the following sections. Table 2. Computational results of LB DP , TWT(π * ), TWT(π * c ), GAP 1 and GAP 2 for instances with I = 10 and 20.

Comparison of constructive heuristics for the APSCSP_TDT |TWT
The nine constructive heuristics were tested on 140 instances for different combinations of I × N × τ × R. Figure 4 summarises computational results of EDD, MST, SPT, SWPT, WMDD, WMDD_C, WMDD_C_I, Greedy and the proposed HWGI for the PSCSP_TDT|TWT. Note that EDD, SPT, SWPT, and MST are to schedule orders by increasing values of d i , (Pinedo 2016;Baker and Trietsch 2019). WMDD_C executes adjacent pairwise interchange from the first position to the last position of the sequence obtained by the WMDD rule and then achieves an improved sequence. WMDD_C_I is to execute an NEHlike heuristic on the sequence generated by WMDD_C. P i is defined by the completion time of the last job in order i if the first job in order i start at time 0. The average GAP and CPU time for different heuristics are depicted in Figure 4. From Figure 4, we observe that, on average, the proposed HWGI generally outperforms the other eight heuristics for all the instances. Whereas, among the nine heuristics, HWGI cost the most CPU time. However, the maximum CPU time required for HWGI is less than 0.5 s, which is tolerable. Therefore, in this paper, HWGI is chosen to generate a better initial solution for each AMHIG algorithm.

Performance comparison of four AMHIG algorithms
In this section, we evaluate the performances of AMHIG_ VNS, AMHIG_VNS_TG, AMHIG_Swap, and AMHIG_ Insert against each other for the different size instances. Due to space limitation, the computational results of 140 instances for different combinations of I and N are presented in Supplemental Data (Tables S6-S9). The results include an average gap of 10 runs (GAP avg ), a minimum gap among 10 runs (GAP min ), and an average relative percentage deviation of 10 runs (ARPD). Table S6 (see Supplemental Data) shows that, for small instances up to 20 orders, 36 out of 40 instances can be solved to optimality by all AMHIG algorithms. The maximum GAP min of each AMHIG is below 0.34  percent for all 40 small instances. It can also be seen from   AMHIGs achieve the same GAP min , GAP avg and ARPD for all instances except for I_30_100_19, I_40_150_5 and I_40_150_18. For I_30_100_19 and I_40_150_18, AMHIG_TG performs slightly worse than the other three algorithms. Whereas, for I_40_150_5, AMHIG_Swap performs slightly worse than the other three algorithms. Tables S8, S9 (see Supplemental Data) and 4 show that, for large instances with I ≥ 60, AMHIG_VNS yields the best or equal GAP avg and ARPD for 51 out of 60 instances, AMHIG_Swap for 21 out of 60 instances, AMHIG_Insert for 13 out of 60 instances, and AMHIG_VNS_TG for 29 out of 60 instances, respectively. AMHIG_VNS obtains the best overall ARPD (0.01), AMHIG_VNS_TG ranks the second (0.02), and AMHIG_Swap and AMHIG_Insert yield the worst overall ARPD (0.04). It seems that AMHIG_VNS performs the best among all the AMHIGs, AMHIG_VNS_TG rank the second, and AMHIG_Swap and AMHIG_Insert performs the worst ( Figure 5).
To determine whether the observed differences are indeed statistically significant, we carried out a series of Wilcoxon signed-rank tests at the significance level α = 0.05, based on the results in Tables S8 and S9 (see Supplemental Data). As we have 60 large instances, the sample size is 60 for each pairwise comparison of algorithms. Let m D denote the median of the difference between the ARPDs generated by two different algorithms, the null hypothesis is defined by H 0 : m D = 0 which means that there is no difference between the ARPDs of two algorithms compared. On the other hand, the alternative hypothesis is defined by H 1 : m D = 0 which means that there is a difference between the ARPDs generated by two algorithms compared. The Wilcoxon signed-rank test results for the pairwise comparisons of algorithms are presented in Table 5. As shown in Table 5, all the pairwise differences are statistically significant at the α = 0.05 level, except AMHIG_Swap vs AMHIG_VNS_TG. For this pair, the difference between the algorithms is not meaningful at the α = 0.05 level, implying AMHIG_Swap and AMHIG_VNS_TG have similar performances. In summary, AMHIG_VNS significantly and statistically performs better than other three AMHIGs. AMHIG_Swap and AMHIG_VNS_TG are statistically equivalent and rank the second. Whereas AMHIG_Insert statistically performs the worst. On the perspective of application, if the user puts more emphasis on the quality of the solution than on CPU time consumption, the AMHDPH is recommended for solving small size instances of the PSCSP_TDT|TWT. Otherwise, the AMHIG_VNS is a good choice when computing time is very limited ( Table 6).

Effect of the HWGI in the AMHIGs
In this section, to reveal the effect of the HWGI at the initial stage in the AMHIGs, we compare the AMHIGs with the AMHIG * s (AMHIGs without HWGI). Note that the same parameters are used in the AMHIG * s as those of the AMHIGs. From Figure 6, we observe that, the overall ARPD values obtained by the AMHIG_VNS and the AMHIG_VNS_TG are obviously smaller than those obtained by the AMHIG_VNS * and the AMHIG_VNS_TG * , respectively. This might be due to the better initial solution generated by heuristics in the AMHIGs. Whereas the AMHIG_Swap and the AMHIG_Insert achieve slightly smaller ARPD values than the AMHIG_Swap * and the AMHIG_Insert * , respectively. It seems that the AMHIG with single neighbourhood structure in local search steps are not much affected by their initial solutions.
To further test the performances of AMHIGs and AMHIG * s, a series of Wilcoxon signed-rank tests at the significance level α = 0.05 are carried out, using the same hypothesis testing method in Section 6.2.3. The sample size is also 60 for each pairwise comparison of algorithms. The Wilcoxon signed-rank test results for the pairwise comparisons of algorithms are listed in Table 7. Table 7 shows that, for each AMHIG and AMHIG * except AMHIG_Insert and AMHIG_Insert * , the pairwise difference is statistically significant at the α = 0.05, which implies that the performance of all AMHIGs except AMHIG_Insert is statistically improved by the heuristic, AMHIG_Insert is not much affected by its initial solution. From Table 7, it also can been seen that the pairwise differences of AMHIG_Swap vs AMHIG_VNS_TG * , AMHIG_Swap vs AMHIG_VNS * , AMHIG_Insert vs AMHIG_Swap * , AMHIG_VNS_TG vs AMHIG_VNS * , AMHIG_Swap * vs AMHIG_Insert * , and ANHIG_Swap * vs AMHIG_VNS_TG * are statistically equivalent. Pairwise differences of other algorithms are statistically significant. In summary, AMHIG_VNS statistically outperforms other algorithms. Whereas the performances of AMHIG_VNS_TG and AMHIG_VNS * are very close and rank the second.

Effect of the correction heuristic (CH) in the AMHIG_VNS
Based on the computational results in Section 6.2.4, in this section, we focus on the effect of the correction heuristic in the AMHIG_VNS. To test whether the best solution π * obtained by the AMHIG_VNS can be further improved if the order sequence is fixed, we develop a hybrid brute force and VNS method (HBV), in which brute force (full enumeration) or VNS are applied to the job sequence of each order according to the number of jobs in the order. Note that the brute force method is to list all possible sequences and choose the sequence with lowest TWT value for the PSCSP_TDT |TWT. The VNS Table 7. Wilcoxon signed-rank test results for the proposed AMHIGs with or without heuristics in initial stage. two models is very small, the probability to obtain the same optimal solution by the two models will greatly increase. For easily understanding, we suppose that there are two solutions π 1 and π 2 . Let π 1 ≺ π 2 denote solution π 1 is better than solution π 2 . The corresponding values evaluated by the original model are TWT(π 1 ) = 128.00 and TWT(π 1 ) = 139.00. Because TWT(π 1 ) < TWT(π 2 ), we have π 1 ≺ π 2 . Whereas if the corresponding values evaluated by the approximate model are TWT A (π 1 ) = 128.13 and TWT A (π 1 ) = 138.82, we also have π 1 ≺ π 2 . This implies that, although there are a small difference on objective values generated by the two models, the dominance relationship between the two solutions remains the same. In this case, the optimal or near-optimal solution for the original considered model can be estimated precisely. Thus, based on this idea, we can search for the optimal or near-optimal solution of the approximate model at first, and then evaluate it by the original model. From a practical point of view, the approximate model and algorithms are proposed based on deepen-analysis on real world data, which is very helpful to solve the original complex scheduling problem efficiently.
The approximate model and algorithms are proposed based on real world data and Proposition 4.1-4.4. In this paper, we assume that all parameters of the problem are known in advance. In real world production environments, problem parameters, e.g. processing time, may be not precise known in advance. However, the dominance relationships in Proposition 4.1 could still hold with a high probability in some cases. For easily understanding, we consider a simple example in which the processing time of a job at each stage is uniformly distributed during the interval [p m,j −p m,j , p m,j +p m,j ], where p m,j is the mean value of the processing time of a job j at stage m and p m,j restricts the bounds of p m,j . By considering the worst case that p m,j = p m,j +p m,j , ∀m ∈ {2, 3, 4}, we define the probability Pr{A} := Pr{p 1,j ≥ m∈{2,3,4} p m,j |p 1,j ≥ m∈{2,3,4} p m,j , p k,j = p k,j +p k,j , ∀k ∈ {2, 3, 4}}. Because p 1,j is uniformly distributed during the interval [p 1,j − p 1,j , p 1,j +p 1,j ], we have Pr{A} = (p 1,j +p 1,j − max { m∈{2,3,4} p m,j , p 1,j −p 1,j })/(p 1,j +p 1,j − (p 1,j −p 1,j )). If m∈{2,3,4} p m,j ≤ p 1,j −p 1,j , then Pr{A} = 1, which means the dominance relationship between processing times at the first stage and its following three stages remains the same; Otherwise, we have Pr{A} = 1 − θ 1,j , where θ 1,j := ( m∈{2,3,4} p m,j + m∈{1,2,3,4}p m,j − p 1,j )/(2p 1,j ). In the latter case, if the value of θ 1,j is very small, e.g. 0.1, the dominance relationships will hold with a high probability (0.9). It's true that if the value of θ 1,j is close to 1, the dominance relationships are unlike to hold. Although the dominance relationships may fail to hold, the proposed approximate methods can be easily integrated into two algorithm frameworks to provide more robust performance guarantee. The first framework is to generate a good initial solution by the approximate methods at first, and then update job sequences and order sequence alternatively by local search-based methods. The second one includes three phases: (1) remove the jobs that do not satisfy the dominance relationships in Proposition 4.1 from their orders; (2) schedule the remaining jobs via the approximate methods; (3) reinsert the removed jobs into the position with the lowest total weighted tardiness.

Performances comparison of the AMHIG_VNS and the hybrid GA
To further verify the effectiveness of the proposed AMHIG_VNS, we tailor the GA-based methods, which are often used in existing works on precast scheduling problems (Chan and Hu 2001;Wang and Hu 2017;Wang, Hu, and Gong 2018a;Ma et al. 2018;Araujo, Bonates, and Prata 2021) to adapt the PSCSP_TDT|TWT. For brevity, the tailored GA-based method is named as HGA. In the HGA, a hybrid GA and VNS is applied to deal with order sequence problem, and a local search method based on swap and insertion operations is used to determine the job sequence in each order. The roulette selection method is used to select individuals from the population. The two-point crossover method and the two-point swap method are adopted for crossover and mutation, respectively. Similar as the AMHIGs, Taguchi method is applied to tune the parameters of the HGA. The parameters for the HGA are set as: PS = 20, p c = 0.9, and p m = 0.1. The process of parameter tuning can be found in Supplemental Data. The run time for each instance is also set to be 0.6 × Is. The differences between the AMHIG_VNS and the HGA for different instance scale groups are provided in Figure 7. Note that X_Y denotes the instance group with X orders and Y jobs. As shown in Figure 7, the AMHIG_VNS outperforms the HGA in term of ARPD values for all scale instance groups. Another observation is that the average performance of the AMHIG_VNS is very robust around 0.4% across various instance groups.

Managerial insights and limitations
The studied problem is motivated by the real scenario of precast industry. Thus, scheduling mould manufacturing, precast production and transportation operations simultaneously is a direct example of the study. However, it can also be realised in many other sectors in which resources or components preparation, production and distribution operations are need to be cooperated. As an example, car manufacturing firms are fabricating components, assembling them into cars, and delivering cars to meet customers' tight due dates. Electronic firms with make to order policy might consider the time-dependent transportation times and use the model to compromise the waiting times and the transportation times of products. Another example is the production and distribution of chemotherapy drugs for cancer treatment. An important factor determining the quality of service of chemotherapy treatment is the time the patient must wait to receive his or her injection of the chemotherapy drug. Chemotherapy production and delivery can be modelled as a production scheduling problem with multiple stages combined with a transportation scheduling problem.
The potential benefits brought by the usage of the model and algorithms can be summarised as follows: • Integrated decisions will help decision-makers to prepare resources, production, and deliver orders in time, resulting in an increase in service quality level. • Results of the study enables decision-makers to be able to compromise production, inventory, and transportation times of orders, which is crucial in real-life problems. • The time to make a decision on the complex supply chain scheduling problem can be greatly shortened while the solution quality is guaranteed.
There are also some limitations of the study.
• In real life, production and distribution operations involve many types of uncertainties such as machine breakdowns, stochastic processing times or transportation times, emergency coming orders, and random customer demand. Our model cannot handle these uncertainties. • Our model is single factory-oriented. Multiple factories in distributed production environments are not considered.

Conclusions and future research
This paper studies a novel precast supply chain order scheduling problem with time-dependent transportations times to minimise the total weighted tardiness of all orders. At first, inspired from real case study and data extracted from a precast manufacturer, we approximate the original problem into a three-stage scheduling problem based on the structural properties. Subsequently, we explore some properties of the schedules for the approximate problem. To solve the PSCSP_TDT|TWT, we propose an AMHDPH for dealing with small size instances. Whereas for medium-or large size instances, we propose an effective AMHIG framework by integrating problem properties, constructive heuristics, an IG and a CH. Based on the framework, we develop four AMHIG algorithms and provided a Lagrangian relaxation bound for performances evaluation especially for medium or large instances. To the best of our knowledge, this study is the first attempt to deal with the PSCSP_TDT |TWT.
Computational results show that the AMHDPH can solve most of small instances to optimality within 60s. Whereas AMHIG can deal with the problems with an average gap below 5 percent in around 70 out of 100 medium or large instances. Statistically tests show that the AMHIG_VNS outperforms other three AMHIG algorithms. The effectiveness of CH at the last stage in AMHIG_VNS is also verified by testing 140 instances. Finally, we test an HBV and RIG on the best solution obtained by the AMHIG_VNS and no solution is improved, which further verify the effectiveness of the approximate model and the proposed algorithms.
It should be noted that, in this paper, we focus on deterministic scheduling models and algorithms for a precast supply chain scheduling problem (PSCSP). At the steady state in transportation stage, the transportation time can be estimated precisely. Thus, in this case, our model can be used in a PSCSP. Our work provides the basis for the PSCSP with more complicated situations such as uncertain transportation times. Thus, to deal with the PSCSP with uncertain transportation time, a potential direction for future studies is to integrate robust optimisation methods, e.g. scenario-based methods (Bertsimas and Sim 2004) into our model and algorithms.
The results are helpful for the precast manufacturers to make planning and scheduling decisions to get more profit and improve customer's satisfaction. Our methods can also be extended to solve similar order scheduling problems especially with a bottleneck stage and time-dependent transportation times. Potential directions for future studies include incorporating additional real-world constraints, such as limited buffer size and emergency orders, machine breakdown, existing in precast supply chain systems into our model; multiple production lines and distributed multiple precast factories will also be discussed; it is also worth developing tighter lower bound methods for evaluating the AMHIG methods especially for large instances; based on previous tighter lower bound methods, we will develop more efficient matheuristics by integration of metaheuristic and exact method, e.g. branch and cut.  Xi'an, China, in 2002, 2006. He is currently a professor with the School of Information and Control Engineering, Xi'an University of Architecture and Technology, Xi'an, China. His research interests include combinatorial optimisation, machine learning and pattern recognition. Linlin Li received the B.E. degree in Software Engineering from Xi'an University of Architecture and Technology, Xi'an, China. She is currently working toward the Master degree with the School of Information and Control Engineering, Xi'an University of Architecture and Technology, Xi'an, China. Her research interests include production planning, scheduling, and supply chain management.