## 2 Problem Formulation
### Reformulated Model Based on Chance-Constrained Method  
Since Model_1 includes uncertainty constraints from Equations (4)-(7), it can be classified as a mixed-integer linear programming (MILP) problem with two-sided uncertainty. Clearly, for Model_1 to be valid, all possible values of the processing time \(\widetilde{P}_{ik}\) must satisfy these constraints. However, in actual production scheduling, it is rare for the processing times of all operations to correspond to the worst-case scenario, making this constraint overly conservative. Therefore, the probability of constraint feasibility from chance-constrained programming is introduced to reformulate Equations (4)-(7) in Model_1. Taking Equation (4) as an example, its chance-constrained expression is as follows:  
\[\mathrm{P}\bigg{\{}S_{jk}+M\cdot\left(3-z_{ik}-z_{jk}-x_{ijk}\right)\geq S_{ik }+\widetilde{P}_{ik},\forall i,j;k\in K_{i}\cap K_{j}\bigg{\}}\geq 1-a \tag{11}\]  
where \(1-a\) represents the confidence level at which the constraint is satisfied. It indicates that, regardless of how the processing time \(\widetilde{P}_{ik}\) varies, the above constraint must be met with a probability of at least \(1-a\). The solution obtained under this constraint is the robust solution at a risk level \(a\) (where \(0\leq a\leq 1\)). For Equation (11), if the probability distribution function \(F_{\widetilde{P}_{ik}}\left(\cdot\right)\) of the uncertain processing time \(\widetilde{P}_{ik}\) is known, then Equation (11) can be equivalently transformed into  
\[S_{jk}+M\cdot\left(3-z_{ik}-z_{jk}-x_{ijk}\right)\geq S_{ik}+F_{\widetilde{P} _{ik}}^{-1}(1-a),\quad\forall i,j;k\in K_{i}\cap K_{j} \tag{12}\]  
where \(F_{\widetilde{P}_{ik}}^{-1}(1-a)\) represents the quantile function of \(\widetilde{P}_{ik}\). In previous research, we proposed kernel density estimation (KDE) for estimating uncertain parameters [24]. Given historical data of uncertain parameters, KDE learns the probability density function of the samples according to  
\[p(x)=\frac{1}{nh}\sum_{i=1}^{n}\varphi\left(\frac{x-x_{i}}{h}\right) \tag{13}\]  
where \(p(x)\) denotes the estimated density function of uncertain parameter \(x\); \(h\) is a positive parameter representing the bandwidth; \(n\) is the number of samples; and \(\varphi(\cdot)\) is a kernel function, which satisfies  
\[\varphi(u)\geq 0 \tag{14}\]  
\[\int\varphi(u)du=1 \tag{15}\]  
In our research, the Gaussian function commonly used in engineering problems is introduced as the kernel function, and the value of bandwidth \(h\) is determined through pilot experiments. Building on this work, the processing time of operations can be described through statistical analysis of their historical data. The range of \(\widetilde{P}_{ik}\) can be described as follows:  
\[\widetilde{P}_{ik}=P_{ik}+\xi_{ik}\hat{P}_{ik}^{max}\quad\xi_{i}\in[-\sigma_{ ik},\sigma_{ik}] \tag{16}\]  
\[\sigma_{ik}=\frac{\hat{P}_{ik}^{1-\eta}}{\hat{P}_{ik}^{max}} \tag{17}\]  
where \(\hat{P}_{ik}^{max}\) represents the upper limit of the magnitude of processing time fluctuations, and \(\hat{P}_{ik}^{1-\eta}\) denotes the upper limit of the magnitude of processing time fluctuations at a confidence level of \(1-\eta\). Both \(\hat{P}_{ik}^{max}\) and \(\hat{P}_{ik}^{1-\eta}\) can be obtained through statistical analysis of historical data, as shown in Figure 2.  
Figure 2: Uncertain interval of \(\widetilde{P}_{ik}\).