# Integration of discrete-event dynamics and machining dynamics for machine
tool: modeling, analysis and algorithms

Mason Ma  Alisa Ren  Christopher Tyler  Jaydeep Karandikar  Michael Gomez
Tony Shi<EMAIL_ADDRESS>Tony Schmitz  Department of Industrial and Systems
Engineering, University of Tennessee Knoxville, Knoxville, TN 37996, USA
Manufacturing Science Division, Oak Ridge National Laboratory, Oak Ridge, TN
37830, USA  MSC Industrial Supply Company, Knoxville, TN 37932, USA
Department of Mechanical, Aerospace and Biomedical Engineering, University of
Tennessee Knoxville, Knoxville, TN 37996, USA

###### Abstract

Machining dynamics research lays a solid foundation for machining operations
by providing stable combinations of spindle speed and depth of cut.
Furthermore, machine learning has been applied to predict tool life as a
function of cutting speed. However, the existing research does not consider
the discrete-event dynamics in machine shop, i.e., the machine tool needs to
process a series of parts in queue under various practical production
requirements. This paper addresses the integration of discrete-event dynamics
and machining dynamics to achieve cost savings in machining. We first propose
a learning-based cost function for the studied integrated optimization problem
of machine tool. The proposed cost function utilizes the predicted tool life
under different stable cutting speeds for further optimizing speed selection
of machine tool to deal with the discrete-event dynamics in machine shop.
Then, according to the practical production requirements, we develop several
mathematical optimization models for the related integrated optimization
problems with the consideration of cost, makespan and due date. Numerical
results show the effectiveness of our proposed methods and also the potential
to be used in practice.

###### keywords:

Machine tool , machining dynamics , discrete-event dynamics , mathematical
optimization , operational excellence

##  1 Introduction

The daily operations of machine tools aim to produce products to fulfill
orders with various production requirements at the minimum cost. This involves
both the machining dynamics that govern the material cutting by the machine
tool, and the discrete-event dynamics that govern the production of a series
of parts queued at the machine tool. Studies in the two fields have
independently made important contributions to the operational excellence of
machine tools in production environments.

The research of machining dynamics has laid a solid foundation for the
operational excellence of machine tool by providing a set of stable cutting
parameters, including spindle speed. The related studies focus on the physical
constraints and modeling of machining processes [20, 21], the stable and
chatter machining parameters identification [1, 8], tool wear and tool life
prediction [9, 10], cutting force modeling and measurement [13, 6] and etc.
The machining of a single part is typically considered in these research
studies. However, few research efforts have considered the practical
production requirements for hundreds or even thousands of parts from multiple
accepted orders.

In the meantime, the discrete-event dynamics research for manufacturing
typically considers the production scheduling of a series of parts that are in
the queue of the machine tool. The operational excellence is usually achieved
at the machine shop level by considering various production requirements, such
as order scheduling on parallel machines [23, 18, 17], production scheduling
on batch machines [16, 14], buffers for machines [25, 12], due date
requirement [4] and etc. However, these studies do not consider the physical
constraints and practical operating of the machine tools, such as the
machining stability and tool life, which can significantly impact the accuracy
of processing times. Most existing studies assume that the processing times
are fixed and known. As such, the solutions from these studies are often
impractical.

Therefore, the integration of machining dynamics and discrete-event dynamics
is critical and promising. In the literature, a few studies have considered
the production scheduling of machine tool, such as machine tool feed rate
scheduling for processing parts with various geometry [3, 2], computer
numerical control (CNC) machine scheduling with controllable processing times
constrained by tool life [11, 7], tool replacement scheduling to balance cost
of machining conditions, product quality and tool life [22, 24]. However, to
our knowledge, no research to date addresses how to integrate machining
dynamics, including the machining stability, various cutting speeds and tool
life, into discrete-event dynamics for the operational excellence of machine
tool. We seek to fill this gap.

In this paper, we address the problem of integration of discrete-event
dynamics and machining dynamics, denoted by IDM, to increase the profit of
machine shop. The machine shop wants to reduce the total cost for producing
all fulfilled orders by selecting an appropriate cutting speed for each order.
We analyze the problem under two scenarios of the practical production
constraints with increasing complexities: 1) The integrated optimization of
cost and makespan; 2) The integrated optimization of cost and due date. While
only optimizing cost is the foundation, the above two scenarios increase the
fidelity of our studied problem to practical production. We show that we can
take advantage of high speed machining and proper production schedule to
reduce the overall cost.

We organize the rest of the paper as follows. In Section 2, we present the
detailed problem settings for both the discrete-event dynamics and the
machining dynamics. A learning-based cost function for machine tool is given
in Section 3. Different production scenarios are studied in Sections 4 and 5.
The mathematical optimization models, analysis and algorithms are provided. In
Section 6, we present numerical results for the proposed methods. Section 7
concludes this paper.

##  2 Problem Settings

In this section, we describe our problem settings of IDM in terms of orders
and parts, stable spindle speed, processing time and tool life, and machine
shop production and the cost.

###  2.1 Discrete-Event Dynamics: Order and Part Setting

We specify the production for $n$ accumulated orders, each consisting of a
number of parts. For example, a machine shop may receive more than 10 orders
from the master production plan for daily operations. We address how to
schedule the production of the machine tool by identifying the best speeds for
each of these orders and scheduling their sequences in the queue of the
machine tool to save total production cost.

Mathematically, we define the set of $n$ orders as
$\mathcal{J}=\\{1,2,\dots,n\\}$. Each order consists of a specific quantity of
parts, denoted by $l_{j}$, $j\in\mathcal{J}$. In this paper, we assume all the
parts of an order are the same type, but parts for different orders can be
different. The different part types can be resulted by different situations,
such as material type, part geometry, and etc. For the single type of part in
an order $j$, let the volume to be removed by machining for each part be
$v_{j}$, such that the total volume to be removed (i.e. workload of order $j$)
can be calculated by $v_{j}l_{j}$. All the orders are released at time zero. A
due date $d_{j}$ is associated for order $j$ before which all parts of the
order must be produced.

The challenges for such discrete-event dynamics research are that the
processing times of orders are various due to different cutting speeds for
different part types. To obtain the operational excellence for machine tool,
the critical issue is to obtain the accurate estimation of processing time for
orders.

###  2.2 Machining Dynamics: Spindle Speed and Tool Life Setting

Machining dynamics provides the foundation for accurate processing time
estimation via the machining stability and tool life estimation under a wide
variety of spindle speeds. Without loss of generality, we adopt the milling
settings throughout this paper. The stability lobe diagram (SLD) enables to
select the best stable spindle speeds at increased axial depth of cut, as
shown in Fig. 1. The machining process can either be stable, that are the
green circles below the stability boundary (blue curve), and chatter (red
cross). In practice, for a given tool-material combination, the SLD can be
obtained by conducting impact test. For more details, see the SLD and its
governing relationship in [15]. In this paper, we assume the stable spindle
speeds are selected from SLD generated by impact test.

Figure 1: Stability lobe diagram (SLD) to determine the stable and unstable
cutting (chatter) by choosing parameters spindle speed ($\Omega$) and axial
depth of cut ($b$).

Given a stable spindle speed $\Omega$ by SLD and a part with the volume to be
removed as $v$, the machining time then can be calculated as

$t=\frac{v}{MRR}=\frac{v}{abf_{t}N_{t}\Omega},$ (1)

where MRR is the material removal rate defined as the average volume of
material removed per unit time with a chosen stable spindle speed $\Omega$
from the SLD, $a$ and $b$ the radial and axial depth of cut, and $f_{t}$ is
the feed per tooth. In this paper, we calculate $t$ by considering various
$\Omega$ with fixed $a$, $b$ and $f_{t}$ as the same setting in [10].

Notably, the tool life is decided directly given by the machining parameters,
like the spindle speed $\Omega$. In general, this is given by the Taylor’s
tool life equation [19]:

$VT^{N}=C,$ (2)

where the cutting speed $V$ is the peripheral velocity by $V=\pi D\Omega$, $D$
is the diameter of the tool, and $T$ is the tool life. $N$ and $C$ are
constants. Whilst higher spindle speed can increase the machining efficiency
of jobs with less machining time $t$, it leads to reduction of tool life.
Thus, there exists a balance of tool life and the machining time which will
eventually impact the cost. To obtain the tool life model is usually by
machine learning methods. In this paper, we use a state-of-the-art machine
learning model in [10] to generate the tool life as input for our method, as
discussed in Section 3.

###  2.3 Production and Cost Settings

For production setting, we specify there is a single machine tool. There are
large quantity of tool edges of the same type that are capable of processing
all the $n$ orders. In practice, it is often infeasible to complete a
fractional number of a part for a tool edge as it requires to change the tool
edge in the middle of the cut. Typically, when the rest tool life is not long
enough for machining a single part, the tool edge needs to be changed with a
new one. The orders are produced sequentially in the queue of the machine
tool. No preemption and interruption of orders are allowed. Furthermore, the
worker would use the same spindle speed to machine all the parts in one order,
regardless of changing tool edges in between the machining process. This is
realistic and easy to implement by the on-site workers, because the quality of
all parts within an order need to be maintained uniformly by keeping the same
cutting conditions.

The objective is to minimize the overall cost introduced by machine tool
operating cost for machining parts plus the tool edges. In [10], a cost model
per part was provided as

$C_{p}=r_{m}t+\frac{(r_{m}t_{ch}+C_{te})t}{T},$ (3)

where we denote $r_{m}$, $t_{ch}$, $C_{te}$ and $C_{p}$ as the machine tool
operating cost in $ per unit time, the time required for changing a tool edge,
the cost in $ per tool edge and the total cost per part, respectively. This
model, however, considers the cost for a single part, rather than the overall
cost of orders for hundreds or even thousands of parts made by different
materials.

In this study, we extend the above model to represent the overall cost for
producing multiple parts in multiple orders as

$C_{o}=\sum_{j\in\mathcal{J}}\big{[}r_{m}p_{j}+(r_{m}t_{ch}+C_{te})q_{j}\big{]},$
(4)

where the first term $r_{m}p_{j}$ calculates the machine tool operating cost
for processing order $j$ and $p_{j}$ refers to the total machining time with a
chosen speed. The second term $(r_{m}t_{ch}+C_{te})q_{j}$ calculates the cost
induced by changing tool edges, and $q_{j}$ is the number of tool edges needed
for processing order $j$.

Eq. (4) is the generalized form of Eq. (3) to include all orders in the
practical machine shop production. In next section, we will present a new cost
function through the calculation of $p_{j}$ and $q_{j}$ under the context of
multiple orders in machine shop.

##  3 Learning-based Cost Function for Machine Tool

The determination of cost of the orders depends on the machining time and
number of tool edges used, which is directly decided by the prediction of tool
life. In this section, we present a new cost function based on the tool life
prediction in the sense of expectation with a state-of-the-art machine
learning model.

###  3.1 A Machine Learning Model for Tool Life Prediction

Karandikar et al. (2021) [10] developed a physics-guided logistic regression
model to predict tool wear probability with cutting speed $V$ and cutting time
$T$ as inputs. Through introducing into machine learning Taylor’s power law-
based tool life description (see Eq. (2)), the model has a powerful
performance in prediction accuracy. The logistic decision boundary of the
model, that is, when the probability of tool worn is 0.5, provides an accurate
estimation for tool life in the sense of expectation as

$\log(T)=-\frac{\theta_{1}}{\theta_{2}}\log(V)-\frac{\theta_{0}}{\theta_{2}},$
(5)

where the coefficients $\theta_{0}=-89.57$, $\theta_{1}=13.57$ and
$\theta_{2}=5.26$ are estimated through experimental data for a single-insert
endmill (Kennametal KICR073SD30333C) with $D=19.05$ mm diameter and a square
uncoated carbide insert (Kennametal 107,888,126 C9 JC). Equivalently, this
model can be transformed into the below Taylor’s tool life model after using
logarithmic transformation,

$T=\bigg{(}\frac{C}{V}\bigg{)}^{1/N}=\bigg{(}\frac{C}{\pi
D\Omega}\bigg{)}^{1/N},$ (6)

where $N=0.388$ and $C=735.2$ m/min using the above values of $\theta_{0}$,
$\theta_{1}$ and $\theta_{2}$. The logistic model in the logarithmic and
original space of $V$ and $T$ is shown in the two panels of Fig. 2. The gray
scale indicates the probability of tool worn with a cutting speed $V$
(horizontal axis) and cutting time $T$ (vertical axis). We will use Eq. (6)
throughout this paper for tool life prediction.

Figure 2: Probability of tool worn as a function of cutting speed from the
trained logistic model; the train data points are shown in light and the test
data points are shown in dark. The left panel shows the results in the
logarithmic space and the right panel in the original space. Figures are
reproduced from [10].

###  3.2 Calculation of Learning-based Cost Function

Based on the learning-based tool life model in [10], a cost function is
developed for multiple orders. As the parts for different orders can be made
by different materials, the SLD of the same tool edge for different orders
(parts) are different. As a result, the candidate stable spindle speeds for
orders are also different. Thus, we define the matrix of spindle speed as

${\bf\Omega}=\left[\begin{array}[]{ccccc}\Omega_{11}&\Omega_{12}&\cdots&\Omega_{1n}\\\
\Omega_{21}&\Omega_{22}&\cdots&\Omega_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
\Omega_{m1}&\Omega_{m2}&\cdots&\Omega_{mn}\\\ \end{array}\right],$ (7)

where the entry $\Omega_{ij}$ refers to the spindle speed for processing order
$j$ with spindle speed $i$, $i\in\mathcal{I}$, and
$\mathcal{I}=\\{1,2,\dots,m\\}$ is defined as the index set of available
stable spindle speed.

For a given spindle speed $\Omega_{ij}$ to process order $j$, the tool life
$T_{ij}$ is calculated by Eq. (6). Also, the total machining time $p_{ij}$ for
processing order $j$ with spindle speed $\Omega_{ij}$ calculated as

$p_{ij}=\frac{v_{j}l_{j}}{abf_{t}N_{t}\Omega_{ij}},$ (8)

where $v_{j}l_{j}$ indicates total volume of order $j$ by multiplying its part
quantity $l_{j}$ with the part volume $v_{j}$. When arranged $T_{ij}$ and
$p_{ij}$ in matrix form, we can obtain the below matrices for tool life and
machining time:

${\bf T}=\left[\begin{array}[]{ccccc}T_{11}&T_{12}&\cdots&T_{1n}\\\
T_{21}&T_{22}&\cdots&T_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
T_{m1}&T_{m2}&\cdots&T_{mn}\\\ \end{array}\right]\text{ and }{\bf
P}=\left[\begin{array}[]{ccccc}p_{11}&p_{12}&\cdots&p_{1n}\\\
p_{21}&p_{22}&\cdots&p_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
p_{m1}&p_{m2}&\cdots&p_{mn}\\\ \end{array}\right],$ (9)

Using $\bf T$, we can now compute the number of tool edges needed for each
order with different spindle speeds as matrix

${\bf Q}=\left[\begin{array}[]{ccccc}q_{11}&q_{12}&\cdots&q_{1n}\\\
q_{21}&q_{22}&\cdots&q_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
q_{m1}&q_{m2}&\cdots&q_{mn}\\\ \end{array}\right],$ (10)

where each entry is calculated by

$q_{ij}=\Biggl{\lceil}\frac{l_{j}}{\big{\lfloor}T_{ij}/t_{ij}\big{\rfloor}}\Biggl{\rceil},\quad\text{and}\quad
t_{ij}=\frac{v_{j}}{abf_{t}N_{t}\Omega_{ij}}.$ (11)

The term $t_{ij}$ is the machining time for a single part of order $j$ with
spindle speed $\Omega_{ij}$. Thus, with floor function, $\lfloor
T_{ij}/t_{ij}\big{\rfloor}$ refers to the number of parts of order $j$ that
can be finished by a single tool edge. Finally, the number of tool edges for
processing $l_{j}$ parts can be obtained by the ceiling function.

We can further obtain the summation of total machining time and tool edge
change time for each order under all its candidate spindle speeds. This
summation, denoted by $h_{ij}$, represents the total processing time occupied
by order $j$ if with speed $i$:

$h_{ij}=p_{ij}+t_{ch}q_{ij}.$ (12)

The associated total time matrix is denoted by

${\mathbf{H}}=\left[\begin{array}[]{ccccc}h_{11}&h_{12}&\cdots&h_{1n}\\\
h_{21}&h_{22}&\cdots&h_{2n}\\\ \vdots&\vdots&\ddots&\vdots\\\
h_{m1}&h_{m2}&\cdots&h_{mn}\\\ \end{array}\right].$ (13)

Finally, with matrices $\bf P$, $\bf Q$ and $\bf H$, we can obtain the cost
function for processing all orders as

$f_{cost}(x)=\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}(r_{m}h_{ij}+C_{te}q_{ij})x_{ij}.$
(14)

Note here we introduce an indicator $x_{ij}$ defined as

$x_{ij}:=\left\\{\begin{aligned} &1,\quad\text{ if speed $i$ is chosen for
processing order $j$},\\\ &0,\quad\text{ otherwise.}\end{aligned}\right.$ (15)

Eq. (14) integrates both the machining dynamics, that is, term
$r_{m}h_{ij}+C_{te}q_{ij}$ based on the cost function for parts of orders, and
the discrete-event dynamics, that is, discrete variable $x_{ij}$ in Eq. (15)
for choosing different speeds for different orders to achieve the reduction of
production cost.

In next section, we will optimize the speed selection through controlling
$x_{ij}$ under the discrete-event dynamics environment of machine shop such
that $f_{cost}(x)$ can be minimized for processing all orders. Before that, we
show an illustrative example.

###  3.3 Illustration of Cost Saving by the Proposed Integration

As an illustrative example of our proposed integration, we show two scenarios
in Fig. 3 for processing an order: with faster speed in (a) obtained by the
expected cost method in [10], and slower speed in (b) by solving our developed
M1 model in Section 4. The light grey period is for tool machining, the light
green period is for changing the tool edge. The red dashed sections refer to
the waste time that is not enough for machining a single part before its tool
life comes to end, while the one in the last tool edge is wasted because all
parts are done.

Figure 3: Illustrative example for choosing different speeds to process
multiple parts of an order by the expected cost method in [10] and solving our
developed model M1 in Section 4. By choosing slower speed, we process one more
part per tool edge in (b) in comparison with (a). The last three parts on the
4th tool edge of (a) are thus moved back to the 3rd tool edge. This leads to
saving the last tool edge and reducing the waste time for each tool edge in
(a), and importantly, reducing the cost from $244.18 to $223.14. See details
in Section 6.

By choosing slower speed in (b), we process one more part per tool edge in (b)
in comparison with (a), namely from 9 part per tool edge to 10. The last three
parts on the 4th tool edge of (a) are thus moved back to be processed with the
3rd tool edge, saving the last tool edge and reducing the waste time for each
tool edge in (a). Notably, the production cost is reduced from $244.18 to
$223.14. The makespan is increased from 80.82 min to 81.09 min as a sacrifice
of the reduced cost. Such trade-off between cost and makespan is meaningful in
practical production. Note Figure 3 only gives the case when smaller speed
will lead to cost reduction, it is more often and realistic to reduce cost by
increasing the speed for high speed machining, and this exactly shows our
motivation for this study.

##  4 Integrated Optimization of Cost and Makespan

We denote our problem when only cost is considered as IDM-C, which is studied
in Sections 4.1 \- 4.2 . In section 4.3, we study the IDM problem when
integrated optimization of cost and makespan is considered, denoted by IDM-CM.
Integer optimization (IO) models are presented for both problems. For clarity,
we first summarize all the notations that will be used for all the models
proposed in this paper throughout.

Sets:

$\mathcal{J}:$

set of orders, $\mathcal{J}=\\{1,2,\dots,n\\}$. Each order $j$ consists of a
quantity of $l_{j}$ parts that are the same part type with the volume to be
removed as $v_{j}$;

$\mathcal{I}:$

index set of stable candidate spindle speed for arbitrary order,
$\mathcal{I}=\\{1,2,\dots,m\\}$;

$\bf\Omega:$

${\bf\Omega}\in\mathbb{R}^{m\times n}$ is the matrix of spindle speed where
entry $\Omega_{ij}$ refers to the spindle speed $i$ for processing order $j$,
$\forall j\in\mathcal{J}$, $i\in\mathcal{I}$. A column $\Omega_{\bm{\cdot}j}$
refers to all candidate spindle speeds that are stable for processing order
$j$;

$\mathbf{T}:$

$\mathbf{T}\in\mathbb{R}^{m\times n}$ is the matrix of tool edge life if
selecting spindle speed $\Omega_{ij}$ for order $j$, $\forall
j\in\mathcal{J}$, $i\in\mathcal{I}$. Each entry $T_{ij}$ is calculated by Eq.
(6);

$\mathbf{P}:$

$\mathbf{P}\in\mathbb{R}^{m\times n}$ is the matrix of total machining time
for all orders where entry $p_{ij}$ refers to the total machining time if
selecting $\Omega_{ij}$ to process order $j$, as in Eq. (8);

$\mathbf{Q}:$

$\mathbf{Q}\in\mathbb{Z}_{+}^{m\times n}$ is the matrix of the number of tool
edges needed for order $j$ with entry $q_{ij}$ calculated by Eq. (11);

$\mathbf{H}:$

$\mathbf{H}\in\mathbb{R}^{m\times n}$ is the matrix of total processing time
where entry $h_{ij}$ consists of two parts: total machining time $p_{ij}$ and
the tool change time $t_{ch}q_{ij}$, as calculated by Eq. (12).

Parameters:

$r_{m}:$

machine tool operating cost ratio in $ per unit time;

$t_{ch}:$

time required for changing a tool edge;

$C_{te}:$

cost in $ per tool edge.

Decision variables:

$x_{ij}:$

$x_{ij}=1$ if spindle speed $\Omega_{ij}$ is selected for processing order
$j$, $\forall j\in\mathcal{J}$ and $i\in\mathcal{I}$. Otherwise, $x_{ij}=0$.

###  4.1 Model with Only Cost Consideration

When only cost is under consideration, the below IO model M1 aims to find the
optimal spindle speed for each order, such that the cost function can be
minimized.

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M1) s.t.
$\displaystyle\sum_{i\in\mathcal{I}}x_{ij}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (16) $\displaystyle
x_{ij}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j\in\mathcal{J}.$ (17)

The objective is to minimize the overall cost $f_{cost}(x)$ of producing $n$
orders. Constraints (16) ensure that each order must be processed by only one
spindle speed. Constraints (17) provide the binary variable restriction on
$x_{ij}$.

We show the constraint matrix of formulation M1 is totally unimodular. Thus,
the IO model can be solved as a linear programming (LP) model. To do this, we
arrange the variables in the sequence of
$x_{11},x_{21},\dots,x_{m1};x_{12},\dots x_{m2};\dots;x_{1n},\dots,x_{mn}$ as
columns of the constraint matrix. The constraint matrix
$\mathbf{A}\in\mathbb{Z}^{n\times mn}$ for constraints (16) is as below:

$\bf A=\bNiceMatrix[first-row]1^{th}&\dots m^{th}(m+1)^{th}\dots
2m^{th}\dots((n-1)m+1)^{th}\dots nm^{th}\\\ 1\dots 10\dots 0\dots 0\dots 0\\\
0\dots 01\dots 1\dots 0\dots 0\\\
\vdots\vdots\vdots\vdots\vdots\vdots\vdots\vdots\vdots\vdots\\\ 0\dots 00\dots
0\dots 1\dots 1\\\ .$ (18)

The above 0, 1 matrix $\bf A$ is totally unimodular because it satisfies the
following two conditions: (1) there are at most two nonzero elements in each
column of $\bf A$, which is number 1 and there is only one entry for each
column to be 1; (2) $\bf A$ admits an equitable row-bicoloring. In particular,
the rows of $\bf A$ can be partitioned into two sets, say “red set” as row 1
and “blue set” as rows 2 to n. It’s easy to check that the sum of the red rows
minus the sum of the blue rows is a vector with entries $\pm 1$.

As a result, the polyhedron for formulation M1 is integral, which has integral
extreme points. Thus, it always admits an integral solution and is equivalent
to solving the LP relaxation:

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M1-LP) s.t.
$\displaystyle\sum_{i\in\mathcal{I}}x_{ij}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (19) $\displaystyle
x_{ij}\in[0,1],~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j\in\mathcal{J},$ (20)

where we relax the binary variable constraints $x_{ij}\in\\{0,1\\}$ into the
interval $x_{ij}\in[0,1]$. Next, we develop an exact algorithm.

###  4.2 Decomposition-based Greedy Algorithm

Formulation M1 naturally adopts a decomposition for both the objective and the
constraints for each order. While every order must be assigned with a spindle
speed, there are no additional constraints that couples the assignment of
spindle speed for different orders. Therefore, the assignment of spindle speed
for each order is independent. That is, for each order $j$, we can select the
optimal speed from $m$ candidates that minimize its cost
$r_{m}h_{ij}+C_{te}q_{ij}$ in the objective. In general, we denote the
described procedure as Decomposition-based Greedy (DG) and summarize as the
below Algorithm 1.

Input : Data matrices $\bf\Omega$; order set $\mathcal{J}$ and associated
$l_{j}$ and $v_{j}$; the cost parameters $r_{m}$, $t_{ch}$ and $C_{te}$.

1 Initialize $f_{cost}(x)=0$; $x_{ij}=0,\forall
i\in\mathcal{I},j\in\mathcal{J}$;

2 Calculate matrices $\mathbf{T}$, $\mathbf{Q}$, $\mathbf{P}$ and $\bf H$ with
Eqs. (6) - (12);

3 for  _$j\in\mathcal{J}$_ do /* Decompose as $n$ subproblems */

4 for  _$i\in\mathcal{I}$_ do /* Find best speed for $j$ */

5 Compute $Obj_{ij}=r_{m}h_{ij}+C_{te}q_{ij}$;

6

7 end for

8 Find the best speed $i^{*}$ of job $j$ by
$i^{*}:=\operatorname*{arg\,min}_{i\in\mathcal{I}}Obj_{ij}$;

9 Update $x_{i^{*}j}=1$;

10 $f_{cost}(x)\leftarrow f_{cost}(x)+Obj_{i^{*}j}$

11 end for

Output Optimal cost $f_{cost}(x)$, speed $x_{ij},\forall i\in\mathcal{I}$,
$j\in\mathcal{J}$

Algorithm 1 The DG algorithm

For Algorithm 1, we have the below lemma to show the optimality of its
solution.

######  Lemma 4.1.

Given $n$ orders and $m$ candidate speeds, Algorithm 1 solves problem IDM-C
optimally in $O(18mn+n)$ time.

###### Proof.

To prove DG algorithm is optimal, we do it by contradiction. Suppose the
solution of DG is not optimal, then there must exist at least one order,
denoted by $\hat{j}$, that selects a spindle speed with smaller cost
$r_{m}h_{ij}+C_{te}q_{ij}$ than $\hat{j}$. This contradicts with steps 7-9 of
DG, which selects the spindle speed with the smallest cost. Thus, Algorithm 1
is an exact algorithm.

Next, we consider the solution time of Algorithm 1. For preparing
$\mathbf{T}$, $\mathbf{Q}$, $\mathbf{P}$ and $\bf H$, it takes 2 operations to
calculate $T_{ij}$, 6 operations to calculate $q_{ij}$, 3 operations to
calculate $p_{ij}$ and 2 operations to calculate $h_{ij}$. Together with the
coupled loops, this takes $14mn$ times. For each $j$, it takes $3m$ operations
to compute $Obj_{ij}$ for all $i\in\mathcal{I}$, $m$ operations to find the
speed $i^{*}$ with minimum $Obj_{ij}$, and 1 operation to update the overall
cost $f_{cost}(x)$, respectively. To sum up, we can get the total time
complexity of Algorithm 1 is $O(18mn+n)$, which completes the proof. ∎

###  4.3 Integration with Makespan

Makespan constraint refers to the hard constraint that all the orders must be
completed within a given time period. For example, this can describe the fixed
time length for the worker shift. To model problem IDM-CM, we need an
additional parameter:

$c:$

a given constant of the upper bound for the makespan.

With this additional parameter, we can model problem IDM-CM by modifying
formulation M1 with makespan constraint as below IO model, denoted by M2:

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M2) s.t.
$\displaystyle\sum_{i\in\mathcal{I}}x_{ij}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (21)
$\displaystyle\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}\leq c,$
(22) $\displaystyle
x_{ij}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j\in\mathcal{J}.$ (23)

Constraint (22) ensures the makespan of producing all the orders not exceeding
the given time length $c$. The makespan constraint adds additional complexity
for the problem. We remark the complexity of problem IDM-CM is still open.

##  5 Integrated Optimization of Cost and Due Date

We consider problem IDM when the integrated optimization of cost and due date
is to be achieved, denoted by IDM-CD. In this case, the sequence of orders
needs to be considered. We provide three categories of model based on the
choices of decision variables. That are, completion time variables (5.1),
linear order variables (5.2) and positional variables (5.3) of orders,
respectively. In the first category, we propose a mixed-integer optimization
(MIO) model. For each of the latter two, we present a quadratic integer
optimization model and an MIO model by applying linearization technique.

###  5.1 Modeling with Completion Time Variables

In the first model, we use completion time variables to indicate the order
sequence for problem IDM-CD. The needed additional variables are defined as
follows:

$y_{jk}:$

$y_{jk}=1$ if order $j$ is processed before order $k$, $\forall
j,k\in\mathcal{J}$ and $j\neq k$. Otherwise $y_{jk}=0$;

$C_{j}:$

completion time of order $j$, $\forall j\in\mathcal{J}$.

The problem IDM-CD can be formulated as below.

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M3-1) s.t.
$\displaystyle\sum_{i\in\mathcal{I}}x_{ij}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (24) $\displaystyle
C_{j}\geq\sum_{i\in\mathcal{I}}h_{ij}x_{ij},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (25) $\displaystyle
C_{j}+\sum_{i\in\mathcal{I}}h_{ik}x_{ik}\leq C_{k}+M(1-y_{jk}),~{}\forall
j<k\in\mathcal{J},$ (26) $\displaystyle
C_{k}+\sum_{i\in\mathcal{I}}h_{ij}x_{ij}\leq
C_{j}+My_{jk},~{}~{}~{}~{}~{}~{}~{}\forall j<k\in\mathcal{J},$ (27)
$\displaystyle C_{j}\in[0,d_{j}],~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (28) $\displaystyle
x_{ij},y_{jk}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j,k\in\mathcal{J}.$ (29)

Constraints (25) ensure that the completion time of each order is no later
than its processing time. Constraints (26) and (27) are disjunctive
constraints which ensure that either order $j$ is processed before order $k$,
i.e. $y_{jk}=1$ in constraints (26), or order $k$ is before order $j$, that is
$y_{jk}=0$ in constraints (27). $M$ is a sufficiently large number.
Constraints (28) enforce the completion time is non-negative and no larger
than the due date. In this formulation, the $M$ is a natural upper bound for
$C_{j}$ for arbitrary $j\in\mathcal{J}$. Thus, we can set
$M=\sum_{j\in\mathcal{J}}\max_{i\in\mathcal{I}}\\{h_{ij}\\}$. The total number
of variables for formulation M3-1 is $n^{2}+mn+n$, which includes $n^{2}+mn$
binary variables and $n^{2}+3n$ constraints.

######  Remark 1.

If $d_{j}=c$ for $j\in\mathcal{J}$, problem IDM-CD is reduced to problem IDM-
CM.

To check the above remark is valid, it’s obvious to see that when $d_{j}=c$
for $j\in\mathcal{J}$, then the makespan of any feasible solutions to problem
IDM-CD also satisfy the makespan limit $c$. Thus, problem IDM-CM is a special
case of problem IDM-CD. It is more general for us to study problem IDM-CD.

###  5.2 Modeling with Linear Ordering Variables

####  5.2.1 Proposed Integer Optimization Model

This formulation is based on the linear order variables $y_{jk}$ and $y_{kj}$
that either of them should be equal to 1 for arbitrary $j<k\in\mathcal{J}$.
The following formulation M3-2 is considered.

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M3-2) s.t.
$\displaystyle\sum_{i\in\mathcal{I}}x_{ij}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (30) $\displaystyle
y_{jk}+y_{kj}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall j,k\in\mathcal{J},j<k,$
(31) $\displaystyle y_{jk}+y_{kl}+y_{lj}\leq
2,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall j,k,l\in\mathcal{J},j\neq k\neq l,$
(32)
$\displaystyle\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}y_{jk}+\sum_{i\in\mathcal{I}}h_{ik}x_{ik}\leq
d_{k},~{}~{}~{}~{}~{}\forall k\in\mathcal{J},$ (33) $\displaystyle
x_{ij},y_{jk}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j,k\in\mathcal{J}.$ (34)

Constraints (31) are a set of conflict constraints, which ensure that either
order $j$ is processed before order $k$ or vice versa. Note here $j<k$ because
of symmetry. Constraints (32) indicate the transitivity constraints that
ensure a linear order of three arbitrary orders. Constraints (33) indicate the
starting time of order $k$ plus its processing time must be no larger than its
due date. The total number of variables of formulation M3-2 with $m$ speeds
and $n$ orders is $n^{2}+mn$, which are all binary variables and involve
$n^{3}-2.5n^{2}+3.5n$ constraints.

Note that M3-2 is a quadratically constrained IO formulation due to
constraints (33). Considering that the current optimization solvers like CPLEX
cannot efficiently deal with this kind of nonlinear problem, we adopt
linearization strategies for this formulation in the following section.

####  5.2.2 Reformulation with Linearization Technique

To enhance the solvability of formulation M3-2, we first apply the
linearization technique proposed in [5] to derive an MIO model, which is
equivalent to M3-2. For arbitrary $j\in\mathcal{J}$, $k\in\mathcal{J}$, let
$\alpha_{jk}=(\sum_{i\in\mathcal{I}}h_{ij}x_{ij})y_{jk}$ and substitute the
cross-product $\sum_{i\in\mathcal{I}}h_{ij}x_{ij}y_{jk}$ in M3-2 using the
following inequalities:

$\left\\{\begin{aligned}
&\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M(1-y_{jk})\leq\alpha_{jk}\leq\sum_{i\in\mathcal{I}}h_{ij}x_{ij}+M(1-y_{jk}),\\\
&-My_{jk}\leq\alpha_{jk}\leq My_{jk},\end{aligned}\right.$ (35)

where $M$ is a sufficiently large number. For arbitrary $j,k\in\mathcal{J}$,
if $y_{jk}=1$, then $\alpha_{jk}=\sum_{i\in\mathcal{I}}h_{ij}x_{ij}$, and the
inequalities become
$\sum_{i\in\mathcal{I}}h_{ij}x_{ij}\leq\alpha_{jk}\leq\sum_{i\in\mathcal{I}}h_{ij}x_{ij}$
and $-M\leq\alpha_{jk}\leq M$. Clearly, the equivalence holds. Otherwise, if
$y_{jk}=0$, the inequalities are
$\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M\leq\alpha_{jk}\leq\sum_{i\in\mathcal{I}}h_{ij}x_{ij}+M$
and $0\leq\alpha_{jk}\leq 0$, which indicates $\alpha_{jk}=0$. Thus, we get an
equivalent model to M3-2.

We further reduce the above inequality constraints (35) as

$\left\\{\begin{aligned}
&\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M(1-y_{jk})\leq\alpha_{jk},\\\
&0\leq\alpha_{jk}\leq My_{jk}.\end{aligned}\right.$ (36)

Thus, we obtain a new linearized formulation M3-2-L as below:

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M3-2-L) s.t.
$\displaystyle\text{Constraints }\eqref{cons:M3-2-1}-\eqref{cons:M3-2-3},$
(37)
$\displaystyle\sum_{j\in\mathcal{J}}\alpha_{jk}+\sum_{i\in\mathcal{I}}h_{ik}x_{ik}\leq
d_{k},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall k\in\mathcal{J},$ (38)
$\displaystyle\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M(1-y_{jk})\leq\alpha_{jk}\leq
My_{jk},~{}~{}~{}\forall j,k\in\mathcal{J},$ (39)
$\displaystyle\alpha_{jk}\geq 0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j,k\in\mathcal{J},$ (40) $\displaystyle
x_{ij},y_{jk}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j,k\in\mathcal{J}.$ (41)

This formulation is obtained by eliminating the right-hand-side of the first
inequality and the left-hand-side of the second inequality in (35). The value
of $M$ is set as the maximum total processing time of the orders, that is
$M=\max_{i,j}\\{h_{ij}\\}$. The total number of variables is $2n^{2}+mn$,
which involves $n^{2}+mn$ binary variables and $n^{3}+0.5n^{2}+3.5n$
constraints. For formulations M3-2-L and M3-2, we have the below Lemma 5.1.

######  Lemma 5.1.

For arbitrary optimal solution to M3-2-L, there exists an optimal solution to
M3-2 with the same objective value.

###### Proof.

For arbitrary optimal solution $(x^{*},y^{*},\alpha^{*})$ to M3-2-L, we can
show $(x^{*},y^{*})$ as an optimal solution to M3-2 with the same objective
value.

We first show $(x^{*},y^{*})$ is a feasible solution to M3-2. For arbitrary
$j,k\in\mathcal{J}$, if $y_{jk}^{*}=0$, then $\alpha_{jk}^{*}=0$ based on
constraints (39) and (40). Then we have
$\sum_{i\in\mathcal{I}}h_{ik}x_{ik}^{*}\leq d_{k}$ in constraints (38). This
indicates $(x^{*},y^{*})$ satisfies constraints (33) in formulation M3-2 when
$y_{jk}^{*}=0$. If $y_{jk}^{*}=1$, then by substituting into constraints (39)
first and then (38), we have

$\displaystyle\sum_{i\in\mathcal{I}}h_{ij}x_{ij}^{*}\leq\alpha_{jk}^{*}\leq
M,$ (42)
$\displaystyle\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}^{*}+\sum_{i\in\mathcal{I}}h_{ik}x_{ik}^{*}\leq\sum_{j\in\mathcal{J}}\alpha_{jk}^{*}+\sum_{i\in\mathcal{I}}h_{ik}x_{ik}^{*}\leq
d_{k}.$ (43)

Inequalities (43) indicate $(x^{*},y^{*})$ satisfies constraints (33) when
$y_{jk}^{*}=1$. To sum up, $(x^{*},y^{*})$ is a feasible solution to M3-2.

Next, we prove the objective value of $(x^{*},y^{*})$ is less than or equal to
that of any feasible solutions to M3-2. Let $\mathcal{S}_{2}$ be the solution
space of formulation M3-2, i.e.,

$\mathcal{S}_{2}:=\\{(x,y)~{}:~{}\text{Constraints \eqref{cons:M3-2-1} -
\eqref{cons:M3-2-5}}\\}.$ (44)

As M3-2-L is a relaxation of M3-2, the objective value of the optimal solution
$(x^{*},y^{*},\alpha^{*})$ will be less than or equal to that of any feasible
solutions to M3-2. That is

$f_{cost}(x^{*})=\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}(r_{m}h_{ij}+C_{te}q_{ij})x_{ij}^{*}\leq
f_{cost}(x),\;\forall(x,y)\in\mathcal{S}_{2}.$ (45)

Note $f_{cost}(x^{*})$ is also the objective value of feasible solution
$(x^{*},y^{*})$ to M3-2. This indicates the objective value of $(x^{*},y^{*})$
is no larger than any feasible solutions to M3-2.

Therefore, $(x^{*},y^{*})$ is an optimal solution to M3-2 with the same
objective value $f_{cost}(x^{*})$, which completes the proof. ∎

###  5.3 Modeling with Positional Variables

####  5.3.1 Proposed Integer Optimization Model

This model is based on the positional variables that describe the relationship
between arbitrary order with the position it is assigned in the schedule. The
additional decision variables are defined as follows.

$z_{jk}:$

$z_{jk}=1$ if order $j$ is assigned on position $k$, $\forall
j,k\in\mathcal{J}$. Otherwise $z_{jk}=0$;

$C_{k}:$

completion time of the order at position $k$, $\forall k\in\mathcal{J}$.

A new IO model can be formulated as below.

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M3-3) s.t.
$\displaystyle\sum_{i\in\mathcal{I}}x_{ij}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (46)
$\displaystyle\sum_{k\in\mathcal{J}}z_{jk}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j\in\mathcal{J},$ (47)
$\displaystyle\sum_{j\in\mathcal{J}}z_{jk}=1,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{J},$ (48) $\displaystyle 0\leq
C_{k}\leq\sum_{j\in\mathcal{J}}d_{j}z_{jk},~{}~{}~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{J},$ (49) $\displaystyle
C_{1}\geq\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}z_{j1},$ (50)
$\displaystyle C_{k}\geq
C_{k-1}+\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}z_{jk},~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{J}/\\{1\\},$ (51) $\displaystyle
x_{ij},z_{jk}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j,k\in\mathcal{J}.$ (52)

Constraints (47) and (48) ensure that a particular order can only be assigned
to exactly one position and each position can be assigned with exactly one
order. Constraints (49) ensure the completion time of order at position $k$ is
non-negative and no later than its due time. Constraints (50) and (51)
describe the completion time of the order at position $k$. The total number of
variables in this formulation is $n^{2}+mn+n$, which involves $n^{2}+mn$
binary variables and $6n$ constraints. Note that M3-3 is a quadratically
constrained IO formulation. We will also apply linearization technique to
reformulate M3-3.

####  5.3.2 Reformulation with Linearization Technique

For arbitrary $j,~{}k\in\mathcal{J}$, define a new variable
$\beta_{jk}=(\sum_{i\in\mathcal{I}}h_{ij}x_{ij})z_{jk}$. An equivalent MIO
model can be obtained by substituting
$\sum_{i\in\mathcal{I}}h_{ij}x_{ij}z_{jk}$ in M3-3 and with the below
additional constraints

$\left\\{\begin{aligned}
&\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M(1-z_{jk})\leq\beta_{jk}\leq\sum_{i\in\mathcal{I}}h_{ij}x_{ij}+M(1-z_{jk}),\\\
&-Mz_{jk}\leq\beta_{jk}\leq Mz_{jk}.\end{aligned}\right.$ (53)

Note the equivalence validation of the above inequalities is the same with
(35). For concise presentation we skip it.

We further reduce the above inequality constraints (53) as

$\left\\{\begin{aligned}
&\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M(1-z_{jk})\leq\beta_{jk},\\\
&0\leq\beta_{jk}\leq Mz_{jk}.\end{aligned}\right.$ (54)

Thus, a new linearized formulation M3-3-L can be obtained as:

$\displaystyle\quad\min\quad$ $\displaystyle f_{cost}(x)$ (M3-3-L) s.t.
$\displaystyle\text{Constraints }\eqref{cons:M3-3-1}-\eqref{cons:M3-3-4},$
(55) $\displaystyle C_{1}\geq\sum_{j\in\mathcal{J}}\beta_{j1},$ (56)
$\displaystyle C_{k}\geq
C_{k-1}+\sum_{j\in\mathcal{J}}\beta_{jk},~{}~{}~{}~{}~{}~{}\forall
k\in\mathcal{J}/\\{1\\},$ (57)
$\displaystyle\sum_{i\in\mathcal{I}}h_{ij}x_{ij}-M(1-z_{jk})\leq\beta_{jk}\leq
Mz_{jk},~{}~{}~{}~{}\forall j,k\in\mathcal{J},$ (58)
$\displaystyle\beta_{jk}\geq 0,~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
j,k\in\mathcal{J},$ (59) $\displaystyle
x_{ij},z_{jk}\in\\{0,1\\},~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}~{}\forall
i\in\mathcal{I},j,k\in\mathcal{J}.$ (60)

The value of $M$ is set as the maximum total processing time among all the
orders, that is $M=\max_{i,j}\\{h_{ij}\\}$. The total number of variables is
$2n^{2}+mn+n$, including $n^{2}+mn$ binary variables and $2n^{2}+6n$
constraints. For formulations M3-3 and M3-3-L, we have the below Lemma 5.2.

######  Lemma 5.2.

For arbitrary optimal solution to M3-3-L, there exists an optimal solution to
M3-3 with the same objective value.

###### Proof.

The proof idea is similar to Lemma 5.1. For arbitrary optimal solution
$(x^{*},z^{*},C^{*},\beta^{*})$ to M3-3-L, we will construct
$(x^{*},z^{*},C^{*})$ as an optimal solution to M3-3 with the same objective
value.

We first show $(x^{*},z^{*},C^{*})$ is a feasible solution to M3-3. It is
simple to check $(x^{*},z^{*},C^{*})$ is feasible when $z_{jk}^{*}=0$ by
substituting $z_{jk}^{*}=0$ into (58) and then (56) and (57). Otherwise, when
substituting $z_{jk}^{*}=1$ we have

$\displaystyle\sum_{i\in\mathcal{I}}h_{ij}x_{ij}^{*}\leq\beta_{jk}^{*}\leq
M,~{}~{}~{}~{}\forall j,k\in\mathcal{J},$ (61) $\displaystyle
C_{1}^{*}\geq\sum_{j\in\mathcal{J}}\beta_{j1}^{*}\geq\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}^{*},$
(62) $\displaystyle C_{k}^{*}\geq
C_{k-1}^{*}+\sum_{j\in\mathcal{J}}\sum_{i\in\mathcal{I}}h_{ij}x_{ij}^{*},~{}\forall
k\in\mathcal{J}/\\{1\\}$ (63)

Inequalities (62) and (63) indicate $(x^{*},z^{*},C^{*})$ is a feasible
solution when $z_{jk}^{*}=1$.

Next, similar to Lemma 5.1, we define the solution space of formulation M3-3
as $\mathcal{S}_{3}:=\\{(x,z,C):\text{Constraints \eqref{cons:M3-3-1} -
\eqref{cons:M3-3-7}}\\},$ and we can get that $f_{cost}(x^{*})\leq
f_{cost}(x),\forall(x,z,C)\in\mathcal{S}_{3}$ as M3-3-L is a relaxation of
M3-3. Thus, we conclude $(x^{*},z^{*},C^{*})$ is an optimal solution to M3-3.
∎

Based on Lemmas 5.1 and 5.2, we will directly solve models M3-2-L and M3-3-L
using MIO solver instead of M3-2 and M3-3, respectively. The comparison of the
three MIO formulations is discussed in next section.

##  6 Computational Experiments

Comprehensive computational experiments are conducted for different problem
settings. For problem IDM-C, we compare the performance of methods by solving
models M1 and M1-LP with modern optimization solver CPLEX, DG algorithm and
the expected cost method (ECM henceforth) from [10]. For problem IDM-CM, we
compare methods by solving M2 using CPLEX with ECM. The three MIO models for
problem IDM-CD, i.e., M3-1, M3-2-L and M3-3-L, are solved by CPLEX and
compared in terms of their computational efficiency.

All the methods are coded in Python. The optimization solver CPLEX 20.1.0 is
used for solving all the IO and MIO models with the Python interface docplex.
All the settings for CPLEX remain default except for time limit as 600
seconds. The experiments are deployed on a mobile workstation with Intel(R)
Xeon(R) W-10885M CPU @ 2.40GHz, 128 GB memory, 64 bit Windows 10 Pro operating
system for workstations.

###  6.1 Data Generation

_Machine tool._ The settings for machine tool including the tool geometry,
cutting parameters (except for spindle speed set) and tool life model are the
same with that in [10]. In particular, the tool is a single-insert endmill
with diameter $D=19.05$ mm and $N_{t}=1$. The radial depth of cut is $a=4.7$
mm, the axial depth of cut is $b=3$ mm and the feed per tooth is $f_{t}=0.06$
mm, respectively. The constants in the Taylor tool life equation are $N=0.388$
and $C=735.2$ m/min.

_Order and part._ The number of orders $n$ is from a set
$\\{2,4,6,8,10,12,14,16,18,20\\}$, with one single part type for each order.
The quantity $l_{j}$ of the part for one order is an integer number randomly
sampled from uniform distribution $U[10,200]$. The volume $v_{j}$ to be
removed per part is an integer number randomly sampled from uniform
distribution $U[5000,15000]$ with unit mm3. The makespan limitation $c$ under
each case of $n$ is decided by considering two values: the makespan of the
solution by solving problem IDM-C which only considers the cost, and the
minimum makespan when adopting the minimum $h_{ij}$ for each order $j$. We set
$c$ as the integer number by rounding the average of these two values.

_Due date._ We specify the due date by generating the data without regarding
for feasibility and then discarding instances that have no feasible schedule.
Specifically, we first obtain the expected cutting speed for each order using
the ECM in [10]. Henceforth, the makespan $\bar{C}_{max}$ and the mean value
of total processing time $\bar{h}$ among the $n$ orders can be obtained. The
initial due date $d_{j}$ is generated by randomly sampling from uniform
distribution $U[0.4\bar{C}_{max},0.5\bar{C}_{max}]$. A schedule is then
generated using earliest due date rule based on $d_{j}$, and the completion
time $C_{j}$ for each order is thus obtained. We then check the feasibility of
$d_{j}$, i.e., $d_{j}$ is accepted if $d_{j}\geq C_{j}$; otherwise
$d_{j}=C_{j}+r$, where $r$ is a randomly sampled integer number from uniform
distribution $U[0,\bar{h}]$. By this way of setting $d_{j}$, we ensure at
least the expected speeds are feasible solution to problem IDM-CD.

_Production and cost._ The candidate cutting speed $V$ is evenly spaced in
$[220,420]$ m/min with the interval 2 m/min. Thus, in total $m=100$ speeds are
used for one order. Note here without loss of generality, we use the same set
of cutting speed for all the $n$ orders. The spindle speed matrix $\bf\Omega$
is thus constructed, with the equivalent minimum spindle speed 3676.02 rpm and
the maximum spindle speed 6984.44 rpm. Lastly, the coefficients of the cost
model are set as $r_{m}=2$ $/min, $t_{ch}=2$ min and $C_{te}=\$12.5$ per tool
edge.

For each order number, 10 independent data instances using different random
seeds are generated with the above data generation procedure. Therefore, in
total 100 instances are used. All the reported results are the average values
obtained by the 10 independent instances for each order number.

###  6.2 Results for Problem IDM-C

Table 1: Numerical results of solving problem IDM-C by the expected cost
method (ECM) [10], solving models M1 and M1-LP, and the exact algorithm DG.

[t]  Order # ECM M1 M1-LP DG Obj: Cost Time (s) Obj: Cost CRR (%) Time (s)
Obj: Cost CRR (%) Time (s) Obj: Cost CRR (%) Time (s) 2 2082.55 2.81 1977.23
5.06 0.01 1977.23 5.06 0.01 1977.23 5.06 0.000 4 3301.91 5.58 3174.74 3.85
0.01 3174.74 3.85 0.02 3174.74 3.85 0.000 6 5418.85 8.28 5184.44 4.33 0.01
5184.44 4.33 0.03 5184.44 4.33 0.000 8 6670.49 11.01 6406.96 3.95 0.01 6406.96
3.95 0.03 6406.96 3.95 0.003 10 9688.88 13.87 9314.33 3.87 0.01 9314.33 3.87
0.03 9314.33 3.87 0.003 12 11872.06 16.66 11363.36 4.28 0.01 11363.36 4.28
0.05 11363.36 4.28 0.003 14 13597.67 19.27 13061.32 3.94 0.02 13061.32 3.94
0.05 13061.32 3.94 0.003 16 15417.88 21.77 14768.81 4.21 0.02 14768.81 4.21
0.06 14768.81 4.21 0.005 18 16849.63 24.57 16127.64 4.28 0.02 16127.64 4.28
0.06 16127.64 4.28 0.003 20 20654.23 27.29 19883.29 3.73 0.02 19883.29 3.73
0.08 19883.29 3.73 0.003

  * *

All values in all columns are the average results of 10 independent runs for
each order number.

  * *

All the 10 independent instances for all the order numbers are solved to
optimality with Gap 0.00% by M1, M1-LP and DG (shown by bold numbers).

Table 1 presents the results for solving problem IDM-C by the ECM in [10],
solving models M1 and M1-LP, and the DG algorithm 1. The columns show the
objectives (minimum cost in $ by each method) and the solution time in second.
In the column “CRR”, we calculate the cost reduction ratio (CRR) as

$CRR=\frac{Obj_{ECM}-Obj}{Obj_{ECM}}\times 100\%,$ (64)

where $Obj_{ECM}$ and $Obj$, respectively, refer to the objective value from
ECM and other methods under comparison.

Figure 4: The comparison of cost reduction ratio (CRR) for problem IDM-C. M1,
M1-LP and DG are overlapped as they all obtain the optimal solutions.

As shown in Table 1, solving M1 and M1-LP by CPLEX, and the DG algorithm all
find the optimal solution for all the 100 instances within 0.1 second, shown
as the bold numbers. Remarkably, the DG algorithm preforms as the fastest
among the three, within 0.005 second for all the 100 instances. In comparison
with the ECM, formulations M1 and M1-LP and DG algorithm show that, by
selecting different cutting speeds, the overall production cost can be reduced
by 3% - 5% on average. Figure 4 shows results of CRR in Table 1. Note the
objectives and the associated CRR for formulations M1, M1-LP and DG algorithm
are the same and thus are overlapped.

###  6.3 Results for Problem IDM-CM

Table 2: Numerical results of solving problem IDM-CM by expected cost method
(ECM) [10] and solving model M2.

[t]  Order # Makespan ECM M2 Limit: $c$ Obj: Cost Time (s) Makespan Obj: Cost
Num Gap (%) LP Relaxation Time (s) Makespan 2 633 2082.55 2.81 763.15 2052.42
10 0.00 2034.58 0.00 624.96 4 1027 3301.91 5.58 1209.71 3245.44 10 0.00
3233.46 0.05 1021.47 6 1678 5418.85 8.28 2000.68 5281.53 10 0.00 5271.44 0.08
1672.02 8 2070 6670.49 11.01 2449.62 6532.52 10 0.00 6526.56 0.08 2067.51 10
3007 9688.88 13.87 3580.07 9512.87 10 0.00 9502.75 0.11 3003.93 12 3669
11872.06 16.66 4362.28 11571.73 10 0.00 11566.20 0.08 3667.11 14 4198 13597.67
19.27 5013.83 13364.26 10 0.00 13354.87 0.05 4195.26 16 4775 15417.88 21.77
5687.06 15047.90 10 0.00 15042.45 0.08 4772.70 18 5202 16849.63 24.57 6221.69
16426.20 10 0.00 16421.71 0.05 5199.98 20 6413 20654.23 27.29 7656.49 20286.36
10 0.00 20281.15 0.08 6411.31

  * *

All values in all columns except for “Num” are the average results of 10
independent runs for each order number.

  * *

Column “Num” indicates the number of instances solved to optimality among the
10 instances.

Table 2 shows the results for solving problem IDM-CM. In the column “Makespan
Limit: $c$”, we show the average makespan limit within which all parts must be
finished. Formulation M2 is solved to optimality (with the gap 0.00) for all
the 10 instances under each order number within 0.11 second, shown as the bold
numbers. The LP relaxation results show M2 is tight, which can explain for its
good computational efficiency.

Notably, due to the makespan limit, the minimum cost obtained in solving
problem IDM-C in Table 1 cannot be obtained in solving problem IDM-CM.
Nevertheless, all these optimal solutions satisfy the makespan limit. While
the ECM does not support to set makespan limit, the results here show that our
integration provides the capability to set many practical production
requirements. Furthermore, for all the 10 order numbers, both the cost and
makespan of ECM solutions are larger than those of the M2 solutions. This
indicates by using our method, it is practical for the machine shop to achieve
solutions that lead to both smaller cost and shorter makespan.

###  6.4 Results for Problem IDM-CD

Table 3: Numerical results of problem IDM-CD by solving models M3-1, M3-2-L
and M3-3-L.

[t]  Order # M3-1 M3-2-L M3-3-L Obj: Cost Num Gap (%) Time (s) Obj: Cost Num
Gap (%) Time (s) Obj: Cost Num Gap (%) Time (s) 2 1920.68 10 0.00 0.02 1920.68
10 0.00 0.02 1920.68 10 0.00 0.03 4 4463.55 10 0.00 0.08 4463.55 10 0.00 0.05
4463.55 10 0.00 0.07 6 5412.14 10 0.00 0.25 5412.14 10 0.00 0.07 5412.14 10
0.00 0.14 8 7293.43 10 0.00 0.27 7293.43 10 0.00 0.10 7293.43 10 0.00 1.83 10
11129.40 10 0.00 1.02 11129.40 10 0.00 0.39 11129.40 10 0.01 25.43 12 12700.74
10 0.01 14.17 12700.74 10 0.00 2.66 12700.74 5 3.14 508.27 14 12615.67 8 0.44
163.11 12615.67 10 0.01 10.80 12615.67 1 6.05 601.27 16 13232.52 4 2.52 485.71
13232.52 10 0.01 45.09 13232.52 1 6.06 602.53 18 18001.38 0 6.47 600.30
18001.38 10 0.01 184.39 18001.38 0 8.06 602.33 20 21394.91 0 8.76 600.60
21394.91 8 0.76 380.35 21394.91 0 8.76 602.19

  * *

All values in all columns except for “Num” are the average results of 10
independent runs for each order number.

  * *

Column “Num” indicates the number of instances solved to optimality among the
10 instances.

Table 3 presents the numerical results for solving problem IDM-CD with
formulations M3-1, M3-2-L and M3-3-L, in terms of the objective value, number
of the instances solved to optimality out of 10, gap to the lower bound
returned by CPLEX and the solution time. Obviously, the solution quality of
M3-2-L is the better than that of M3-1 and further better than that of M3-3-L.
This can be seen from the number of instances solved to optimality. Among all
the 100 instances, M3-2-L solves 98 to optimality except for 2 instances at
$n=20$ orders, while that of M3-1 and M3-3-L are 72 and 66, respectively. The
same conclusions can be drawn from the gap and the solution time. For example,
at $n=18$, M3-2-L takes 0.01 second to solve all the 10 instances to
optimality, while M3-1 takes 600.30 seconds yet still has gap 6.47%, and
M3-3-L takes 602.33 seconds with gap 8.06%. Note that time above 600 seconds
means the solving procedure is terminated due to time limit of CPLEX.

It is interesting to see although with different gaps for the three models,
the objectives are all the same. In our experience, it can often be the case
that CPLEX has found the optimal solution in 5 minutes, though it takes much
longer time to certify its optimality. Thus, formulations with fast increasing
lower bounds can certify the optimality in shorter time.

As a closing remark for this section, the results with cost, makespan and due
date consideration validate the capability our methods to integrate
constraints for describing practical requirements for the operational
excellence of machine tool.

##  7 Conclusion, Discussion and Future Work

In summary, we have identified and studied a significant problem of the
integration of discrete-event dynamics and machining dynamics for the
operational excellence of machine tool for the first time. The machining
stability, various cutting speeds and tool life were incorporated into the
production of a series of parts from multiple orders in the queue of machine
tool. Based on a state-of-the-art logistic regression model for tool life
prediction, a new cost function was developed for orders with cutting speed
choices as the variables. Then the cost was minimized via choosing different
stable cutting speeds.

We have demonstrated a set of models and algorithms by considering various
practical production scenarios, including cost, makespan and due date. When
considering the integrated optimization of cost and makespan, we developed
integer optimization models and exact algorithm. Then, we considered the
integrated optimization of cost and due date. Three categories of mixed-
integer optimization models were presented based on the choices of decision
variables. The linearization technique was applied to obtain the models. The
results showed while satisfying the makespan and due date requirements, our
methods achieved on average 3% - 5% reduction of the cost, which has the
potential for cost saving in practice.

Future work includes the deeper integration from both the discrete-event
dynamics and the machining dynamics. Besides cutting speed, other machining
parameters can be as decision variables, including axial depth of cut and feed
rate. Thus, we may have nonlinear and nonconvex constraints from the stability
lobe diagram. Also, the stochastic tool life prediction can be considered. In
the meantime, the complexities of the problems for makespan and due date
constraints are one future direction. The practical extensions of the problems
are to consider different candidate speeds for orders, various part types of
an order, and multiple machine tools in the machine shop.

## Acknowledgements

This work has been supported in part by the University of Tennessee Knoxville
under the Graduate Advancement Training and Education (GATE) program of
Science Alliance, and the Southeastern Advanced Machine Tools Networks
(SEAMTN).

## References

  * [1] Y. Altintas, G. Stepan, E. Budak, T. Schmitz, and Z. M. Kilic.  Chatter stability of machining operations.  Journal of Manufacturing Science and Engineering, 142(11), 2020\.
  * [2] X. Beudaert, P.-Y. Pechard, and C. Tournier.  5-axis tool path smoothing based on drive constraints.  International Journal of Machine Tools and Manufacture, 51(12):958–965, 2011.
  * [3] W. Ferry and Y. Altintas.  Virtual five-axis flank milling of jet engine impellers—part ii: feed rate optimization of five-axis flank milling.  Journal of Manufacturing Science and Engineering, 130(1), 2008.
  * [4] J. M. Framinan and P. Perez-Gonzalez.  Order scheduling with tardiness objective: Improved approximate solutions.  European Journal of Operational Research, 266(3):840–850, 2018\.
  * [5] F. Glover and E. Woolsey.  Converting the 0-1 polynomial programming problem to a 0-1 linear program.  Operations Research, 22(1):180–182, 1974.
  * [6] M. F. Gomez and T. L. Schmitz.  Displacement-based dynamometer for milling force measurement.  Procedia Manufacturing, 34:867–875, 2019.
  * [7] S. Gurel and M. S. Akturk.  Considering manufacturing cost and scheduling performance on a cnc turning machine.  European Journal of Operational Research, 177(1):325–343, 2007\.
  * [8] J. Karandikar, A. Honeycutt, T. Schmitz, and S. Smith.  Stability boundary and optimal operating parameter identification in milling using bayesian learning.  Journal of Manufacturing Processes, 56:1252–1262, 2020.
  * [9] J. Karandikar, T. McLeay, S. Turner, and T. Schmitz.  Tool wear monitoring using naive bayes classifiers.  The International Journal of Advanced Manufacturing Technology, 77(9):1613–1626, 2015.
  * [10] J. Karandikar, T. Schmitz, and S. Smith.  Physics-guided logistic classification for tool life modeling and process parameter optimization in machining.  Journal of Manufacturing Systems, 59:522–534, 2021.
  * [11] R. K. Kayan and M. S. Akturk.  A new bounding mechanism for the cnc machine scheduling problems with controllable processing times.  European Journal of Operational Research, 167(3):624–643, 2005\.
  * [12] H. Ma, H. K. Lee, Z. Shi, and J. Li.  Workforce allocation in motorcycle transmission assembly lines: A case study on modeling, analysis, and improvement.  IEEE Robotics and Automation Letters, 5(3):4164–4171, 2020.
  * [13] T. No, M. Gomez, R. Copenhaver, J. U. Perez, C. Tyler, and T. L. Schmitz.  Force and stability modeling for non-standard edge geometry endmills.  Journal of Manufacturing Science and Engineering, 141(12):121002, 2019.
  * [14] M. Qin, R. Wang, Z. Shi, L. Liu, and L. Shi.  A genetic programming-based scheduling approach for hybrid flow shop with a batch processor and waiting time constraint.  IEEE Transactions on Automation Science and Engineering, 18(1):94–105, 2019.
  * [15] T. L. Schmitz and K. S. Smith.  Machining dynamics: Frequency response to improved productivity, 2nd edition.  Springer, page 303, 2018.
  * [16] Z. Shi, Z. Huang, and L. Shi.  Customer order scheduling on batch processing machines with incompatible job families.  International Journal of Production Research, 56(1-2):795–808, 2018\.
  * [17] Z. Shi, H. Ma, M. Ren, T. Wu, and J. Y. Andrew.  A learning-based two-stage optimization method for customer order scheduling.  Computers & Operations Research, 136:105488, 2021.
  * [18] Z. Shi, L. Wang, P. Liu, and L. Shi.  Minimizing completion time for order scheduling: Formulation and heuristic algorithm.  IEEE Transactions on Automation Science and Engineering, 14(4):1558–1569, 2015.
  * [19] F. W. Taylor.  On the art of cutting metals, volume 23.  American society of mechanical engineers, 1906.
  * [20] C. T. Tyler and T. L. Schmitz.  Analytical process damping stability prediction.  Journal of Manufacturing Processes, 15(1):69–76, 2013.
  * [21] C. T. Tyler, J. R. Troutman, and T. L. Schmitz.  A coupled dynamics, multiple degree of freedom process damping model, part 1: Turning.  Precision Engineering, 46:65–72, 2016.
  * [22] Z. Vagnorius, M. Rausand, and K. Sørby.  Determining optimal replacement time for metal cutting tools.  European Journal of Operational Research, 206(2):407–416, 2010\.
  * [23] L. Wang, Z. Shi, and L. Shi.  A novel quadratic formulation for customer order scheduling problem.  In 2013 IEEE International Conference on Automation Science and Engineering (CASE), pages 576–580. IEEE, 2013.
  * [24] A. Zaretalab, H. S. Haghighi, S. Mansour, and M. S. Sajadieh.  A mathematical model for the joint optimization of machining conditions and tool replacement policy with stochastic tool life in the milling process.  The International Journal of Advanced Manufacturing Technology, 96(5):2319–2339, 2018.
  * [25] C. Zhang, Z. Shi, Z. Huang, Y. Wu, and L. Shi.  Flow shop scheduling with a batch processor and limited buffer.  International Journal of Production Research, 55(11):3217–3233, 2017.