An Innovative Formulation Tightening Approach for Job-Shop Scheduling

Job shops are an important production environment for low-volume high-variety manufacturing. Its scheduling has recently been formulated as an integer linear programming (ILP) problem to take advantages of popular mixed-integer linear programming (MILP) methods, e.g., branch-and-cut. When considering a large number of parts, MILP methods may experience difficulties. To address this, a critical but much overlooked issue is formulation tightening. The idea is that if problem constraints can be transformed to directly delineate the problem convex hull in the data preprocessing stage, then a solution can be obtained by using linear programming (LP) methods without combinatorial difficulties. The tightening process, however, is fundamentally challenging because of the existence of integer variables. In this article, an innovative and systematic approach is established for the first time to tighten the formulations of individual parts, each with multiple operations, in the data preprocessing stage. It is a major advancement of our previous work on problems with binary and continuous variables to integer variables. The idea is to first link integer variables to binary variables by innovatively combining constraints so that the integer variables are uniquely determined by the binary variables. With binary and continuous variables only, it is proved that the vertices of the convex hull can be obtained based on vertices of the LP problem after relaxing binary requirements. These vertices are then converted to tightened constraints for general use. This approach significantly improves our previous results on tightening individual operations. Numerical results demonstrate significant benefits on solution quality and computational efficiency. This approach also applies to other complex ILP and MILP problems with similar characteristics and fundamentally changes the way how such problems are formulated and solved. Note to Practitioners—Scheduling is an important but difficult problem in planning and operation of job shops. The problem has been recently formulated in an integer linear programming (ILP) form to take advantage of popular mixed-integer linear programming methods. Given an ILP problem, there must exist a linear programming (LP) formulation so that all of its vertices are also the vertices to the ILP problem. If such an LP problem can be found in the data preprocessing stage, then the corresponding ILP problem is tight and can be solved by using an LP method without difficulties. In this article, an innovative and systematic approach is established to tighten the formulations of individual parts, each with one or multiple operations. It is a major advancement of our previous work on problems with binary and continuous variables by novel exploitation of the relationship between integer and binary variables in job-shop scheduling. The resulting tightened constraints are characterized by part parameters and the length of the scheduling horizon and can be easily adjusted for other data sets. Results demonstrate significant benefits on solution quality and computational efficiency. This approach also applies to other complex ILP and MILP problems with similar characteristics and fundamentally changes the way how such problems are formulated and solved.

ILP problem is tight and can be solved by using an LP method without difficulties. In this article, an innovative and systematic approach is established to tighten the formulations of individual parts, each with one or multiple operations. It is a major advancement of our previous work on problems with binary and continuous variables by novel exploitation of the relationship between integer and binary variables in job-shop scheduling. The resulting tightened constraints are characterized by part parameters and the length of the scheduling horizon and can be easily adjusted for other data sets. Results demonstrate significant benefits on solution quality and computational efficiency. This approach also applies to other complex ILP and MILP problems with similar characteristics and fundamentally changes the way how such problems are formulated and solved.

I. INTRODUCTION
J OB shops are an important production environment for low-volume high-variety manufacturing. In a job shop, machines are usually categorized into different types based on their functions. With these machines, multiple parts with different due dates are processed, and each part needs a sequence of operations to be completed [1]. To meet on-time deliveries, scheduling of parts is critical. The problem is to minimize the required objective, e.g., the total weighted tardiness and the total cycle time, by assigning parts to machines while satisfying part processing time requirements, and operation precedence and machine capacity constraints. It is one of the hardest scheduling problems. Only a few special cases are polynomially solvable, such as two machines or two parts, and slightly generalized versions of these problems are NP-hard [2]. For practical-sized job-shop problems, the optimal solution is difficult to get and the goal is usually to obtain near-optimal solutions with quantifiable quality.
As reviewed in Section II, some nonlinear job-shop scheduling formulations were established and efficiently exploited by decomposition and coordination methods. To take advantage of popular mixed-integer linear programming (MILP) methods, e.g., branch-and-cut, the problem is recently formulated in an integer linear programming (ILP) form and is not convex 1 with the existence of integer variables. Branch-and-cut first solves the linear programming (LP) problem without integrality requirements. If the solution is feasible to the original MILP problem, then it is optimal. If not, valid cuts are performed around the optimal solution of the LP problem on the fly to get 1 A set is convex if the line segment between any two points in it lies in it. A function defined on a convex set is convex if the line segment in-between any two points on the graph of the function lies above the graph between the two points.
solutions to the MILP problem. If such solutions are obtained, the problem is directly solved. If not, the method relies on time-consuming branching operations. When considering a large number of parts, the state-ofthe-art and practice MILP methods may experience convergence and quality difficulties. To obtain near-optimal job-shop schedules fast, a critical but much overlooked issue is formulation transformation. The idea is to transform problem constraints to directly delineate the convex hull (the smallest convex set that contains all feasible solutions [3]) in the data preprocessing stage. If this can be done (i.e., the formulation is "tight"), then a solution can be obtained by using an LP method without combinatorial difficulties. Theoretically, the tighter the formulation, the less the time to obtain solutions with the same quality. The tightening process, however, is fundamentally challenging without a systematic approach because of the existence of integer variables (e.g., beginning time).
In the literature, a few tightened single-part formulations were reported and shown computationally efficient for the overall problems. In our previous work [4], [5], a systematic formulation tightening approach was developed for mixed-binary linear programming (MBLP) for the first time, and the method was realized based on power system unit comment problems. The idea is to derive vertices of the convex hull of a unit (a generator) without binary requirements. From them, vertices of the original convex hull are innovatively obtained. These vertices are converted to tightened constraints, which are then parameterized based on unit characteristics for general use in unit commitment problems. The method was then extended to job-shop scheduling problems with integer variables in [6]. Our idea is to obtain the vertices of the convex hull of a part without integrality requirements and then approximate noninteger values in vertices as nearest feasible integers. For each operation, two sets of tightened constraints related to processing time were obtained. They were then parameterized by analyzing their patterns for general use in job-shop scheduling problems. Precedence constraints were modeled but not tightened. The results were then extended to consider energy costs as a component of the objective function in energy-efficient job-shop scheduling problems in [7]. Results in [4]- [7] demonstrate computational efficiency and solution quality benefits of formulation tightening.
In this article, the job-shop scheduling problem is first formulated in an integer programming form in Section III following [6], and the objective is to minimize the total weighted tardiness and the total cycle time. Since tardiness is a nonlinear function of part completion times, it is linearized by introducing new binary and continuous variables and the corresponding constraints to make effective use of MILP methods. Then, the problem becomes an MILP problem. Since the complexity of tightening increases with problem sizes, the focus is on individual parts.
To tighten the formulations of parts with multiple operations in the data preprocessing stage, an innovative and systematic approach is established for the first time in Section IV. This is a major advancement of our previous work as it is generalized from MBLP to MILP problems with special structures. The idea is to first link integer variables (e.g., beginning time) to binary variables (e.g., part statues) by innovatively combining constraints so that the integer variables are uniquely determined by the binary variables. With binary and continuous variables only, it is proved that the vertices of the convex hull can be obtained based on vertices of the LP problem after relaxing binary requirements [4], [5]. Thus, there is no approximation in the tightening process as those in [6] and [7]. These vertices are then converted to tightened constraints. The numbers of resulting tightened constraints and variables involved, and the constraint coefficients depend on part parameters. Since all parts must be processed within the scheduling horizon, the above also depends on the length of the horizon. For general use purposes, these tightened constraints are characterized by analyzing constraint structures and the relationship between coefficients and part parameters as well as the scheduling horizon.
Since it is difficult to directly tighten the formulation with multiple operations, the idea is to first tighten the formulation of a single operation to explore relations among part status and beginning/completion time. The resulting processing time and beginning/completion time-related tightened constraints can be applied to every operation of parts with multiple operations. Then, the same method is used to tighten the formulation of two successive operations to explore their interactions. The resulting precedence-related tightened constraints can be used for every two consecutive operations of parts with multiple operations. The process can be repeated for the formulation with three and more operations. The tightening process only needs to be performed once, and the resulting tightened constraints can be easily adjusted for other data sets after parameterization and can be directly applied in the data preprocessing stage, tremendously reducing online computational requirements. This approach significantly improves our previous results on tightening individual operations [6].
Three examples are considered in Section V. The first is to tighten the formulations for single parts to illustrate the tightening idea and present insights. Robustness of formulation tightening is shown in the second example. The last example is to demonstrate the performance of tightened single-part formulations when solving the overall job-shop scheduling problems. Results demonstrate significant benefits on solution quality and computational efficiency.
Beyond MILP job-shop scheduling problems under consideration, this approach also applies to the other complex ILP and MILP problems with similar characteristics between integer and binary variables. It fundamentally changes the way how such problems are formulated and solved. This approach goes naturally with decomposition and coordination approaches, a subject worthy of further exploration.

II. LITERATURE REVIEW
Existing job-shop formulations and solution methodologies are reviewed in Section II-A. Tightened constraints are reviewed in Section II-B.

A. Problem Formulations and Solution Methodologies
With large numbers of decision variables and constraints in job-shop scheduling, developing efficient formulations is complex [8]. "Separable" and nonlinear formulations were established and efficiently exploited by decomposition and coordination-based Lagrangian relaxation methods in [9]- [14]. ILP models were also developed in [15]- [25]. Considering sequence-dependent setups, an ILP model was established in [15]. With additional variables, job successors and predecessors were modeled. In our previous work on high-volume and low-variety manufacturing [25], an ILP model was developed. However, large-scale problems cannot be effectively solved by using those models.
To solve job-shop scheduling problems, branch-and-cut has been widely used. The method solves the LP problem without integral requirements by using an LP method first. If the solution has integer values for all integer decision variables, it is optimal with respect to the original problem. If not, the method tries to obtain the convex hull by adding valid cuts to cut off regions outside the convex hull without cutting off feasible solutions. If successful, the problem is directly solved. If not, time-consuming branching operations are performed, resulting in very slow convergence. In [15] and [17]- [25], the problems were solved by branch-and-cut implemented in commercial software CPLEX or Gurobi. When considering a large number of parts, branch-and-cut may experience convergence and quality difficulties.

B. Tightened Constraints
Obtaining a tight formulation is fundamentally difficult and NP-hard without clear ways. In the literature, few tightening studies exist on general problems. For traveling salesman problems, a tightened formulation was obtained based on subtour elimination [26]. For knapsack problems, tight formulations were obtained through the use of "structural" disjunctive cuts based on the problem structure [27].
For manufacturing scheduling, a few tightened constraints were presented for single parts without explaining how they were generated. For traditional job-shop scheduling, a few valid cuts were developed by analyzing problem structures in [17]. The major idea is to find a ceiling for inventory shortage and the longest working procedure sequence till completion for parts. Testing results based on randomly generated data for 325 instances with 3-5 machines and 4-6 parts demonstrate the computational efficiency of the cuts. For flow-shop scheduling, subtour elimination constraints and lower/upper bound mixed-integer inequalities were developed by analyzing formulation structures in [15], and some of them are facet-defining cuts. Testing results based on randomly generated data for problems with 2-6 machines and 7-10 parts show that the computational time is much reduced with these tightened constraints. For both studies, testing results demonstrate the computational efficiency of these tightened constraints.
In our previous work [4], [5], a systematic formulation tightening approach was developed for MBLP problems for the first time, and the method was realized based on power system unit comment problems. It is based on novel integration of "constraint-and-vertex conversion," "vertex elimination," and "parameterization" processes to tighten single-unit formulations. Our idea is to derive vertices of the convex hull of a unit (a generator) without binary requirements. From them, vertices of the original convex hull are innovatively obtained by eliminating vertices with fractional values for binary variables with rigorous proof. These vertices are converted to tightened constraints, which are then parameterized based on unit characteristics for general use in unit commitment problems. The method was then extended to job-shop scheduling problems with integer variables in [6]. Our idea is to obtain the vertices of the convex hull of a part without integrality requirements and then approximate noninteger values in vertices as nearest feasible integers. For each operation, two sets of tightened constraints related to processing time were obtained. They were then parameterized by analyzing their patterns for general use in job-shop scheduling problems. Precedence constraints were modeled but not tightened. The results were then extended to consider energy costs as a component of the objective function in energy-efficient job-shop scheduling problems in [7]. Results in [4]- [7] demonstrate computational efficiency and solution quality benefits of formulation tightening.

III. JOB-SHOP SCHEDULING FORMULATION
Consider a job shop with multiple machines categorized into M types based on their functionalities. By using these machines, I parts with different due dates need to be processed, and the part index is i . Part i requires J i operations, and the operation index is j . It is assumed that the scheduling horizon is long enough so that all parts can be processed. The horizon is discretized into K time slots and let k denote the time index. Assuming that system-level machine capacity constraints are relaxed, a single-part scheduling problem is formulated based on our previous work [6] in Section III-A (part index i is omitted for brevity). Machine capacity constraints and the objective function are briefly described in Section III-B.

A. Single-Part Formulation
For a part with J operations, the main decision variables are beginning time b j and completion time c j for each operation j , both integer variables. To capture the status of operation j at time k, i.e., active (processed) or not, binary variables δ jmk with operation j , machine group m, and time k indices are considered Here, the machine group concept is considered instead of individual machines since machines are usually categorized into different types based on their functions in a job shop. With machine groups, the decision space is much reduced, and the problem complexity is thus much reduced. The choice of machines within the same group will be based on heuristics after optimization and is not captured in the formulation.
Since the processing time p jm is machine-type-dependent, a new set of binary decision variables is defined to identify the assignment of operation j to machine type m as follows: x jm = 1, if operation j is assigned to machine type m 0, otherwise.
Part-level constraints are processing time requirements, operation precedence constraints, and machine-type assignment constraints. Modeling of linearized tardiness is also included.
1) Processing Time Requirements: Because of the "nonpreemption," a contiguous time period with length of p jm is needed to process operation j , i.e., where M j denotes the set of machine types that can process operation j of part i (part index i is omitted). Since processing time p jm is generally machine-type-dependent, the actual processing time depends on the assignment of machine types.
Since δ jmk represents the status of the part, δ jmk must be 1 within [b j , c j ] if machine type m is assigned to process operation j and 0 otherwise, i.e., Logical (2) is linearized as follows: where N is a big number. It can be seen that (3)-(5) guarantee that δ jmk = 1 iff b j ≤ k ≤ c j and x jm = 1, and δ jmk = 0 otherwise.

2) Operation Precedence Constraints:
It is assumed that the operation sequence of a part is fixed, and operation j + 1 cannot start until j is finished, i.e., Also, the part cannot start the process of the first operation until it arrives at time a, i.e., 3) Machine-Type Assignment Constraints: Operation j can only be assigned to one machine type, i.e., Also, if a machine type cannot process operation j of part i , then the assignment variable should be 0, i.e., 4) Linearized Tardiness: Tardiness T is formulated as follows: where d is the due date. To represent this, a piecewise linear function is used, as shown in Fig. 1.
As shown in the figure, the lower and upper bounds of c J −d are j min m∈M j p jm − d and K − d, and the corresponding tardiness is 0 and K −d. The three break points of this function on the x-axis are , and the corresponding values at the y-axis are 0, 0, and K − d. This piecewise linear function is linearized by using special ordered set techniques [28]. Three continuous variables w 1 , w 2 , and w 3 (0 ≤ w 1 , w 2 , w 3 ≤ 1) are considered to represent the weights of the three points. In addition, three binary variables α 1 , α 2 , and α 3 are used to set up upper bounds for these weights. The constraints are as follows: For simplicity, instead of j min m∈M j p jm − d and K − d, two break points −K and 2K are used for all parts (d could be negative). Here and later in the article, the max function is kept for compactness of notation.

B. Machine Capacity Constraints and Objective Function
For completeness, machine capacity constraints and the objective function are briefly described in this section.

1) Machine Capacity Constraints:
For each machine type m, the total number of active parts cannot exceed machine capacity M m at any time slot, i.e., In the above, (i , j ) denotes operation j of part i , and O m denotes the set of (i , j ) that can be processed by machine type m.
2) Objective Function: The objective function is to minimize the weighted sum of total tardiness and total cycle time, i.e., where ω is the weight for total tardiness and ω T i is for part i . The job-shop scheduling problem with (1), (3)-(9), and (11)-(18) established above is an MILP problem. Most of the  decision variables are binary (e.g., δ and x). There are also a few integers (e.g., b and c) and continuous variables (i.e., w). If every operation of each part can only be processed on one machine type, then there is no need to consider machine-type assignment variable x, and machine-type index m for part status variable δ can be deleted. The machine capacity constraints and objective function are linear but irrelevant for tightening.

IV. FORMULATION TIGHTENING
Building upon our previous work [4]- [7], an innovative and systematic method is established to tighten the above single-part formulation in Section IV-A. A numerical example is also presented to illustrate the tightening idea. Tightness is proved in Section IV-B.

A. Formulation Tightening
In our previous work on power system unit comment problems, a systematic approach is developed to tighten MBLP problems [4], [5]. To illustrate the idea, consider a simple binary linear programming (BLP) problem in Fig. 2 with two binary variables x 1 and x 2 , and x 1 + x 2 ≥ 0.5. After relaxing binary requirements, the vertices [blue dots in Fig. 2 Fig. 2 [5]. These vertices are then converted to tightened constraints for general use. The idea to tighten MBLP problems is the same.
For ease of presentation, the following terms are defined. Definition 1: For an MBLP problem, if the integrality requirements are relaxed, the resulting convex hull is defined as the "integer-relaxed convex hull." In terms of the simple example described above, the integer-relaxed convex hulls are defined by blue lines in Fig. 2 Definition 2: For an integer-relaxed convex hull of an MBLP problem, a vertex consists of integral and real components. If all integral components have integer values, then it is called an "integral vertex." Otherwise, it is called "fractional vertex." In terms of the simple example above, integral and fractional vertices are denoted by solid and open blue dots, respectively, in Fig. 2(b).
The above definitions can apply to an MILP problem. To apply the MBLP tightening idea to tighten the MILP problem under consideration, the relationship that integer variables (e.g., beginning time b) are uniquely determined by binary variables (e.g., part status δ and machine-type assignment variable x) is innovatively established as to be proved in Section IV-B. Therefore, the MBLP principle of eliminating fractional vertices with respect to δ and x described above can be applied.
Since it is difficult to directly tighten the formulation with multiple operations, the idea is to first tighten the formulation of a single operation to explore relations among part status and beginning/completion time. The resulting processing time and beginning/completion time-related tightened constraints can be applied to every operation of parts with multiple operations. Then, the same method is used to tighten the formulation of two successive operations to explore their interactions. The resulting precedence-related tightened constraints can be used for every two consecutive operations of parts with multiple operations. The process can be repeated for the formulation with three and more operations.
1) One Operation: Given part parameters (due date d, processing time p, and arrival time a) and the length of the scheduling horizon (K ) in numerical values, tightened constraints are established by an innovative and systematic method through four steps, as shown in Fig. 3.
Step 1 (Constraint-to-Vertex Conversion): After relaxing integrality requirements, the constraints are converted to the vertices of the integer-relaxed convex hull. The conversion is done by algebraic manipulation of part parameters and the scheduling horizon length with algorithms [29] well established in existing software Porta [30]. The constraints are input, and the output is vertices in numerical values.
Step 2 (Vertex Elimination): If all vertices obtained in Step 1 are integral, the formulation is tight. If not, fractional vertices are projected onto the original convex hull. For this particular problem, all integral vertices of the integer-relaxed convex hull are the same as the vertices of the original convex hull and vice versa, as will be proved in Section IV-B. Thus, vertex projection can be done by eliminating factional vertices.
Step 3 (Vertex-to-Constraint Conversion): In this step, the vertices obtained in Step 2 are converted back to tightened constraints by using Porta as a reverse process of that in Step 1. The resulting formulation with those constraints should be tight.
Step 4 (Parameterization): Constraints obtained above have coefficients in numerical values. To make them reusable for other parts, the idea is to convert numerical coefficients to part parameters (e.g., processing time) and the total number of time slots in the scheduling horizon. This parameterization is done by analyzing constraints and the relationship between numerical coefficients and part parameters and the scheduling horizon length. It is verified by checking the physical meanings of the resulting constraints with coefficients in part parameters and the scheduling horizon length under all possible combinations of binary variables. The resulting tightened constraints can be easily adjusted for problems with other data sets.
For a single part, the number of tightened constraints, the number of variables involved, and constraint coefficients depend on part parameters and the length of the scheduling horizon. For example, consider the first operation of a part with p = 3 and K = 5, and assume that this operation can only be processed on one machine type. Because of "nonpreemption," a contiguous time period with length of 3 is needed to process this operation. If the first time block is taken, then the contiguous time period cannot go beyond time block 3; otherwise, the process is disjunctive. Therefore, δ 1 + δ 4 ≤ 1 and δ 2 + δ 5 ≤ 1. Since the assumption is that the scheduling horizon is long enough so that the operation can be processed, δ 1 + δ 4 = 1 and δ 2 + δ 5 = 1. For the same part with K = 6, there is one more similar constraint. Note that after parameterization, the resulting tightened constraints can be used for individual operations of parts with multiple operations.
Numerical Example: To illustrate the above approach, a numerical example is presented. Consider the first operation of a prat with p = 3 and K = 8, and assume that this operation can only be processed on one machine type. Decision variables include part status δ k , beginning time b, and completion time c. Constraints are processing time requirements (1) and (3)- (5). Without integrality requirements, the constraints to Porta are shown in Fig. 4.
By constraint-to-vertex conversion, 1234 vertices are obtained and the last ten are shown in Fig. 5. Six integral vertices remain after eliminating factional vertices, as shown in Fig. 6. By vertex-to-constraint conversion, tight constraints are generated by Porta in Fig. 7.   Fig. 7 are converted to a set of processing time-related tightened constraints as follows: Because of "nonpreemption," a contiguous time period with length of 3 is needed to process this operation. If the first time block is taken, then the contiguous time period cannot go beyond time block 3, and otherwise, the process is disjunctive. Therefore, δ 1 + δ 4 + δ 7 ≤ 1. Since it is assumed that the scheduling horizon is long enough so that the part can be processed, δ 1 + δ 4 + δ 7 = 1 as shown in (19a). Similarly, one δ from time slots 2, 5, and 8 must be 1 as shown in (19b), and one δ from time slots 3 and 6 must be 1 as shown in (19c). Given (19) and binary requirements of δ, inequalities (11) and (12) in Fig. 7 are redundant. The processing time-related tightened constraint set has been reported in our previous work [6].
The above set of tightened constraints can be generalized for all operations with different processing time as follows: 2) Two Operations: Now, consider two operations. Given part parameters (due date d, processing time p 1 and p 2 , and arrival time a) and the scheduling horizon length in numerical values, tightened constraints are established as follows.
For the first and second operations, they have their own constraints, such as processing time requirements. There is also an operation precedence constraint that couples the two operations together. Denote the operation-level constraints for the first and second operations as C 1 and C 2 , respectively, and the coupling constraint as C 1−2 . Apply the tightened constraints obtained by tightening the single operation formulation to C 1 and C 2 , and obtain TC 1 and TC 2 , respectively. With the constraint set {TC 1 , TC 2 , C 1−2 }, tighten the two-operation formulation through the four steps presented above and obtain tightened constraints across two operations as TC 1−2 . Note that after parameterization, TC 1−2 can be used for every two consecutive operations of parts with multiple operations.
Similar to the tightened constraints for every operation, the tightened constraints across two operations also depend on part parameters and the length of the scheduling horizon. For example, consider a part with p 1 = 3 and p 2 = 1, and K = 5. Because the part must be processed in the scheduling horizon, the latest completion time of operation one is 4 as operation two needs one time slot after it, and thus, δ 1,5 = 0. Similarly, the earliest beginning time of operation two is 4 as operation one needs three time slots before it, and thus, δ 2,1 = δ 2,2 = δ 2,3 = 0.
3) Multiple Operations: With tightened constraints for individual operations and every two consecutive operations, the tightening process is repeated for parts with more operations. Since the number of vertices increases exponentially in constraint-and-vertex conversion and so does the number of constraints, it is difficult to obtain a tight formulation. Our goal is thus to obtain "near-tight" formulations by partially tightening.

B. Tightness Proof
Tightness proof is established in Theorem 1. Proof: The proof will be conducted in two steps. The major step is to show that the values of integer decision variables can be uniquely determined by the values of binary variables. The remaining step is to prove that integral vertices of the integer-relaxed convex hull Conv(P MILP−IR ) are the vertices of the original convex hull Conv(P MILP ) based on the theorems developed for MBLP problems in our previous work [5].

1) Step 1-Integer Variables Can Be Uniquely Determined by Binary Variables:
An integral vertex of the integer-relaxed convex hull is feasible to the original MILP problem. Now, consider one operation. Since "nonpreemptive" processing time requirements modeled in (3)-(5) are satisfied, a contiguous time period of length m∈M x m p m (operation index j is omitted for brevity, and it is assigned to one machine type in set M that can process this operation) should be assigned to process this operation. For any time k 0 such that 1 ≤ k 0 ≤ K -m∈M x m p m + 1, without loss of generality, it is assumed that this operation is processed during the time interval [k 0 , k 0 + m∈M x m p m −1]. Because of (5), m∈M x m p m time slots will be occupied, and because of (3) and (4), these time slots will be contiguous. Since the operation is processed during time slot k ∈ [k 0 , k 0 + m∈M x m p m − 1], then m∈M δ mk equals 1 during these time slots and 0 otherwise as required by the processing time requirements (5). Therefore, the terms with big number N in (3) and (4) disappear. By replacing c as b + m∈M x m p m − 1 in (3) and combining (3) and (4), the following inequality is obtained: To analyze the relationship between b and k 0 , (21) is evaluated at the two extreme points of k in [k 0 , k 0 + m∈M x m p m − 1]. First, set k as k 0 , and (21) Then, set k as k 0 + m∈M x m p m − 1, and (21) becomes Equation (23) can be rewritten as Combining (22) and (24), obtain k 0 ≤ b ≤ k 0 , which implies that b = k 0 . Then, c can be described as k 0 + m∈M x m p m −1. 2) Step 2-Tightness of MILP Problems: For an MBLP problem, it has been proved that the integral vertices of its integer-relaxed convex hull are all vertices of the original convex hull in our previous work [5] and vice versa. Since the MILP problem under consideration can be treated as an MBLP problem in the tightening process, integral vertices of its integer-relaxed convex hull Conv(P MILP−IR ) are the vertices of original convex hull Conv(P MILP ) and vice versa.
End: Based on Theorem 1, vertex projection can be simply done by eliminating factional vertices in Step 2 to tighten the single operation formulation. For parts with multiple operations, since the relations between b, δ, and x within individual operations still hold, the values of b and c can be uniquely determined by the values of δ and x. Thus, the formulation is still tight by applying the same idea as that for the single operation.
Generalization: Beyond MILP job-shop scheduling problems under consideration, this approach also applies to other complex ILP and MILP problems with similar characteristics between integer and binary variables.

V. NUMERICAL RESULTS
The above tightening method is implemented by using Porta [30]. The job-shop scheduling problems are solved by using IBM ILOG CPLEX Optimization Studio V 12.8.0.0 [31] on a PC with 2.40-GHz Intel Xeon E-2286M CPU and 32-GB RAM. Three examples are presented. The first is to tighten the formulations of single parts to illustrate the idea and present the insights. Robustness of formulation tightening is shown in the second example. The last example is to demonstrate performance of tightened single-part formulations when solving the overall problems.
A. Example 1: Single Part 1) One Operation: Consider the scheduling problem for the first operation of a prat with p = 3 and K = 8 used in Section IV-A. Constraints (2), (3), and (5) in Fig. 7 have been explored in Section IV-A, and inequalities (11) and (12) are redundant given the binary requirements of δ. Therefore, the focus is on exploring inequalities (6)-(10) and equalities (1) and (4) in Fig. 7.
The above set of tightened constraints can be generalized for all operations with different processing times as follows: Equation (26) implies that if δ jmτ are all 1, either δ jmk or δ jm,k+ p j m must be 1.
When examining inequalities (6)-(10) in Fig. 7 as a group, they can be put together in the matrix form as follows: ⎛ Note that δ 1 -δ 3 do not show up in (27). The reason is that if δ 4 -δ 8 are properly regulated to satisfy the "nonpreemptive"  processing time requirements by (27), and then, δ 1 -δ 3 are expected to satisfy the requirements because of the processing time-related tightened constraint (19). In addition, although the first row of (27), denoted as , is redundant given the binary requirements of δ, it helps put the constraints together in a square matrix form. Physical meanings of (27-2)-  are analyzed one by one in the following.
Physical meanings of (27-2) under all combinations of binary variables involved are shown in Table I.
It can be seen that  guarantees that if δ 8 is 1, δ 7 must be 1 as implied by the second row of Table I. This is reasonable because the last possible three time slots to process the operation are 6-8, given that the processing time is 3. If the eighth time slot is taken, then the seventh must be taken too, and otherwise, the operation cannot be completed within the scheduling horizon. The physical meaning of (27-3) is similar, and if δ 7 is 1, δ 6 must be 1, given that the processing time is 3.
Physical meanings of  under all combinations of binary variables involved are shown in Table II.
It can be seen that  guarantees that δ 6 cannot be 1 when δ 5 and δ 8 are both 0 as implied by the third row of Table II. In other words, if δ 6 is 1, one of δ 5 and δ 8 has to be 1. This is reasonable because if neither the fifth time slot nor the eighth is taken, then the sixth cannot be taken based on the processing time requirement. Note that Case 8 is not feasible because δ 5 and δ 8 cannot be 1 at the same time as guaranteed by (19b).
Physical meanings of  under all combinations of binary variables involved are shown in Table VI in the Appendix, as well as the analysis. The physical meanings of (27) can be intuitively shown in Fig. 8.
In Fig. 8(a), if δ 8 is 1, then δ 7 must be 1. In Fig. 8(b), if δ 7 is 1, then δ 6 must be 1. In Fig. 8(c), if δ 6 is 1, then one of δ 5 and δ 8 has to be 1. In Fig. 8(d), if δ 5 is 1, then one of δ 4 and δ 7 has to be 1. Equation (27) with the processing time-related tightened constraint (19) guarantees a contiguous time period with a length of 3 to process the operation. For example, if δ 6 is 1, there are two possibilities: 1) δ 8 is 1 or 2) δ 5 is 1. The first situation is simple since δ 7 will be 1 when δ 8 is 1, and thus, time slots 6-8 are used to process the operation. For the second situation, δ 5 is 1, and there are two possibilities again: 1) δ 4 is 1 or 2) δ 7 is 1. For the first situation, time slots 4-6 are used to process the operation. For the second situation, time slots 5-7 are used. Among all situations described above, a contiguous time period with length of 3 is guaranteed to process the operation. A similar analysis can be performed for other δ.
With further analysis on the meanings of (27), it can be extended to a set of tightened constraints (28) for problems with longer time periods as follows: ⎛ The entries y nh of the above matrix are presented as follows.
The meaning of (28) is similar to (27). For example, when δ k is 1, the meaning is intuitively shown in Fig. 9.
For example, if δ k is 1, there are two possibilities: 1) δ k−1 is 1 or 2) δ k+2 is 1. For Case 1), there are two possibilities again: 1) δ k−2 is 1 or 2) δ k+1 is 1. For each of them, a contiguous time period with length of 3 is guaranteed to process the operation. For Case 2), δ k+2 is 1, and there are two possibilities again: 1) δ k+1 is 1 or 2) δ k+4 is 1. For Case (2.1), a contiguous time period with length of 3 is guaranteed to process the operation. For Case (2.2), according to (28), one of δ k+3 and δ k+6 has to be 1. This is contradictory with (20), where δ k+3 and δ k+6 have to be zero as δ k is 1. Among all the situations described above, a contiguous time period with length of 3 is guaranteed. A similar analysis can be performed for other δ.
Since the processing time is p and the operation must be completed within the scheduling horizon, the largest beginning time is K − p + 1 with δ 6 -δ 8 as 1 as implied in (34). When the starting of nonzero δ moves earlier, b gets smaller. The earlier δ, the larger the negative impacts on b. The meaning of (35) is similar. The largest completion time is K with δ 6 -δ 8 as 1. When the starting of nonzero δ moves earlier, c gets smaller. It can be verified that these constraints are meaningful under all possible part statuses of δ 1 -δ 8 .
The above two tightened constraints can be generalized for all operations with different processing times as follows: Equations (20), (26), (36), and (37) directly constrain δ 1 -δ 8 , b, and c within one operation. For the single-operation part problem with p = 3 and K = 8 under consideration, with (19) and (25), the total number of vertices decreases from 1234 to 250 in a major way. With (34) and (35), it is further reduced to 42. After replacing (25) by (27), the total number of vertices is 6, and all of the vertices are integral vertices, implying that the formulation is tight. For the single operation scheduling problem with the processing time as 1-4, tight formulations are obtained.
2) Two Operations: Now, add another operation with a processing time of 3 to the problem in Section V-A1, with additional operation precedence constraint (6). For each operation, decision variables include b, c, and one set of δ k . With the standard formulation established in Section III, after relaxing integrality requirements, 63 872 vertices are obtained by constraint-to-vertex conversion. After applying (20) and (26) obtained in the one-operation example to both operations, a total number of 23 206 vertices remains. With (36) and (37) applied to both operations, 333 vertices remain. After eliminating factional vertices, there are six integral vertices. After vertex-to-constraint conversion, the resulting tight constraints are obtained, as shown in Fig. 10.
After analyzing the meanings of equalities (2)-(7) in Fig. 10 under possible part statuses, two sets of operation precedence-related tightened constraints are obtained as follows: Since the two operations need to be completed in the scheduling horizon, the largest completion time for operation two is K , with the beginning time of Kp 2 + 1. Therefore, operation one must be completed by that time, and δ 1,k must be 0 for period of k ∈ [K − p 2 + 1, K ] as implied in (38). The meaning of (39) is similar. The smallest beginning time of operation one is 1, with completion time of p 1 . Therefore, operation two cannot start before p 1 , and δ 2,k must be 0 for period of k ∈ [1, p 1 ].
Equations (38) and (39) directly constrain δ 1,k and δ 2,k across operations. With them, the total number of vertices decreases to 14 from 333 in a major way.
The above two tightened constraints can be generalized for all operations with different processing times as follows: 3) One Operation With Linearized Tardiness: Now, add the constraints associated with tardiness to the problem in Section V-A1, assuming due date d is 2. Constraints under consideration are processing time requirements (1) and (3)-(5) and tardiness constraints (11)- (16). Decision variables include b, c, a set of δ k , a set of continuous variables w, a set of binary variables α, and tardiness T . With the standard formulation established in Section III, after relaxing integrality requirements, 6170 vertices are obtained by constraintto-vertex conversion. With (19), (25), (34), and (35) obtained in the one-operation example, a total number of 210 vertices remain. Eliminate factional vertices, and obtain eight integral vertices. After vertex-to-constraint conversion, the resulting tight constraints are shown in Fig. 11.
By combining equalities (2), (5), and (10) in Fig. 11, the following constraint is obtained: Since the processing time is 3 and the due date is 2, the smallest tardiness is 1 with δ 1 -δ 3 as 1, and the largest tardiness is 6 with δ 6 -δ 8 as 1 as implied in (42). When the starting of nonzero δ moves earlier, T gets smaller. The earlier δ, the smaller the impacts on T . Equation (42) directly constrains δ 1 -δ 8 and tardiness T . With it, the total number of vertices decreases to 84 from 210.
After analyzing the physical meanings, (42) is converted to the following constraint in a generic form (43), as shown at the bottom of the next page: The tightened constraints obtained in Sections V-A1-V-A3 tighten the formulation. However, they can hardly be obtained manually without going through the above tightening process. The formulation with them is much tighter (not tight yet) than the original one. Those tightened constraints can be extended to other parts with more operations and processing time other than 3, whose due date is positive [only for (43)].

B. Example 2: Medium-Sized Problems
This medium-sized example is to demonstrate the effectiveness and robustness of formulation tightening. The instance is  [9]. The largest number of operations of parts is 6. According to which parts/operations that machines can process, machines are categorized into 17 types, and each type has 1-6 machines with the same function. Assume that each operation of each pat can only be processed on one machine type, and machines are always available for simplicity. The number of time slots under consideration is 220 so that all the parts can be processed. There are three values for tardiness weights, 1, 10, and 100, and they are randomly assigned to parts with the percentage of 50%, 40%, and 10%, respectively. The weight for the total tardiness is 0.9.
Before and after adding tightened constraints, the overall job-shop scheduling problems are solved by using branch-andcut implemented in CPLEX. In CPLEX, optimization stops when computational time reaches the preset stop time or the relative Mixed-Integer Programming (MIP) gap (relative difference between the objectives of the optimal relaxed solution and current integer solution) falls below the preset gap. Here, the stop MIP gap is set as 0.01% and there is no time limit.
With different formulations, the results are presented in Table III According to Table III,  When replacing (26) by (28), (30), and (32) (the formulation becomes tighter), a similar solution is obtained in 6 s, with cutting time of 0.9 s. The reason is that when solving the problem by using branch-and-cut, cuts are performed around the optimal solution to the integer-relaxed problem [31], not on the entire feasible region as mentioned in Section I. Therefore, more tightened constraints may not guarantee better computational efficiency. There is a tradeoff between tightness and computational efficiency.
The problem is also solved with randomly assigned part tardiness weights considering formulations (a), (e), and (f) presented above. Cutting and branching time for problems with different sets of weights are compared in Fig. 12. Then, for each part, a random variable following U (−5, 5) is generated and added to the part due date, while tardiness weights are the same as the original problem. The results are shown in Fig. 13.
For every instance, solutions with the same quality are obtained with and without tightened constraints. By adding tightened constraints, the total cutting and branching time is significantly reduced, and the reduction is mainly from the reduction of branching time. Results demonstrate the effectiveness and robustness of our formulation tightening.

C. Example 3: Large-Sized Problems
This large-sized example is to demonstrate the performance of tightened single-part formulations. The instance is created based on all 127 parts and all machines in [9]. The number of time slots under consideration is 300 so that all the parts can be processed. The other problem setup is the same as that in  Example 2. Before and after adding tightened constraints, the job-shop scheduling problems are solved by using branch-andcut, and the results are shown in Table IV. Since the problem is more complicated than Example 2, two stopping criteria are considered: 1200-s CPU time or 0.01% MIP gap. According to Table IV, both the solution quality and computational efficiency are significantly improved by adding tightened constraints. With the standard formulation, no feasible solution can be found in 20 min. By adding new tightened constraints (20), (26), (36), (37), (40), (41), and (43), a feasible solution with an MIP gap of 0.01% is obtained in 12 s. The results show that tightening single parts also improves the solution quality and computational efficiency when solving the overall problems. Similar to Example 2, when replacing (26) by (28), (30), and (32) (the formulation becomes tighter), both the CPU and solving time increases. Now, consider that some of the 127 parts can be processed on multiple machine types. The other problem setup is the same as that in the above. Before and after adding tightened constraints, the job-shop scheduling problems are solved by using branch-and-cut, and the results are shown in Table V. According to Table V, both the solution quality and computational efficiency are significantly improved by adding   tightened constraints. With the standard formulation, no feasible solution can be found in 20 min. By adding new tightened constraints, a feasible solution with an MIP gap of 0.01% is obtained in 50 s.
The above results demonstrate the great potential of our formulation tightening method for complex ILP and MILP problems where the values of integer variables are uniquely determined by the values of binary variables.

VI. CONCLUSION
In this article, an innovative and systematic method is established for the first time to tighten the formulations of individual parts with multiple operations in the data preprocessing stage. It is a major advancement of our previous work on problems with binary and continuous variables to integer variables. The idea is to first link integer variables to binary variables by innovatively combining constraints so that the integer variables are uniquely determined by the binary variables. With binary variables only, the vertices of the convex hull can be obtained based on the vertices of the linear problem after relaxing binary requirements with proved tightness. These vertices are then converted back to tightened constraints with coefficients characterized by part parameters and the length of the scheduling horizon. The tightening process only needs to perform once, and the resulting tightened constraints can be easily adjusted for other data sets after parameterization and can be directly applied in the data preprocessing stage, tremendously reducing online computational requirements. This method significantly improves our previous results on tightening individual operations. Numerical results demonstrate significant benefits on solution quality and computational efficiency.
Beyond MILP job-shop scheduling problems under consideration, this approach also applies to other complex ILP and MILP problems with similar characteristics between integer and binary variables, such as job-shop scheduling problems with other features like sequence-dependent setups. For practical applications, the idea is to obtain "near-tight" formulations by partially tightening. The approach fundamentally changes the way how such problems are formulated and solved. In addition, this method goes naturally with part-based decomposition and coordination approaches, a subject worthy of further exploration.

A. Physical Meanings of Tight Constraints
Physical meanings of  under all combinations of binary variables involved are shown in Table VI.
It can be seen that  guarantees that δ 5 cannot be 1 when δ 4 and δ 7 are both 0 as implied by the fifth row of Table VI. In other words, if δ 5 is 1, one of δ 4 and δ 7 has to be 1. This is reasonable because if neither of the fourth time slot or the seventh is taken, the fifth cannot be taken based on the processing time requirement. Note that Case 10 in Table VI is not feasible because if δ 8 is 1, δ 7 must be 1 as guaranteed by . Cases 11,12,15, and 16 are not feasible as guaranteed by (19). For Case 3, since δ 6 cannot be 1 when δ 5 and δ 8 are both 0 as guaranteed by , δ 6 has to be 0. However, if δ 7 is 1, δ 6 must be 1 as guaranteed by . Therefore, Case 3 is not infeasible either.
ACKNOWLEDGMENT Any opinions, findings, and conclusions or recommendations expressed in this material are those of the author(s) and do not necessarily reflect the views of NSF or DoE.