## 3 Mathematical Formulation
### Sample Average Approximation Methodology  
The two-stage chance-constrained optimization model in the above section can be solved by the sample average approximation approach. In this approach, the true distribution of each uncertain parameter is approximated by a set of independent samples by using a Monte Carlo simulation while the corresponding sample average function replaces the objective function. Many studies in the literature show that SAA effectively solves a chance-constrained optimization problem [28; 29; 34; 35; 36; 37]. However, it can be seen that the more uncertain parameters, the larger the size of the optimization problem and the longer the computing time. To overcome this challenge, a common technique to improve computational efficiency is using clustering methods such as K-means, Fuzzy C-means methods [38; 39]. This technique reduces a large number of initial samples into a small number of clusters, then these clusters are represented as a new set of samples and used to solve the SAA problem. Algorithm 1 outlines the main steps of the SAA algorithm combined K-means approach to solve a general chance-constrained optimization problem as follows:  
\[V=\min\{f(x)+\mathbb{E}(Q(y,\xi))\} \tag{25}\]  
subject to  
\[\Pr(G(x,y,\xi)\leq 0)\geq 1-\varepsilon \tag{26}\]  
where \(x\) is the first-stage variable, \(y\) is the second-stage variable, \(\xi\) is random input data, and \(\varepsilon\) is the risk level of the chance constraint in the problem (25). In the first step, we generate \(M\) independent sample sets of size \(N\). Then, the k-means clustering approach is applied to divide each set into \(N_{L}\) clusters. Each cluster's centroid will be considered a scenario in SAA algorithm with its probability is equal to the total probabilities of all samples inside the cluster. Consequently, the optimization problem in Equations (25) and (26) are reformulated as:  
\[V=\min\Biggl{\{}f(x)+\sum_{n=1}^{N_{L}}p_{n}Q(y_{n},\xi_{n})\Biggr{\}} \tag{27}\]  
subject to  
\[\sum_{n=1}^{N_{L}}p_{n}\mathbf{1}_{(0,\infty)}(G(x,y_{n},\xi_{n}))\leq\gamma \tag{28}\]  
where \(\mathbf{1}_{(0,\infty)}(G(x,y_{n},\xi_{n}))\) is equal to one if \(G(x,y_{n},\xi_{n})\geq 0\) and zero otherwise, \(p_{n}\) is \(n\)th centroid's probability (\(n=1,2,\ldots,N_{L}\)), and \(\gamma\) is the risk level of the SAA model (27) \((0\leq\gamma\leq\varepsilon)\). Based on the method presented in [29; 40], the SAA problem (27) is implemented \(S\) iterations. In each iteration, this model will be solved M times in correspondence to M sample sets so that we can obtain M first-stage optimal solutions \(x_{s}^{m}\) and optimal values \(V_{s}^{m}\) (\(1\leq s\leq S,\;1\leq m\leq M\)). As the optimal problem contains chance constraint, so we need to check the feasibility of the first-stage solution \(x_{s}^{m}\) by evaluating the \((1-\beta)\)-confidence upper bound of the chance constraint with \(N^{\prime}\) samples \((Nt\gg N)\) as follows:  
\[U(x_{s}^{m})=g(x_{s}^{m})+\Phi^{-1}(1-\beta)\sqrt{\frac{g(x_{s}^{m})(1-g(x_{s} ^{m}))}{N^{\prime}}} \tag{29}\]  
where \(g(x_{s}^{m})=\Pr(G(x_{s}^{m},y,\xi)\leq 0)=\frac{1}{N^{\prime}}\sum_{n=1}^{N^{ \prime}}\mathbf{1}_{(0,\infty)}(G(x_{s}^{m},y_{n},\xi_{n}))\); and \(\Phi^{-1}\) is the inverse normal distribution function. If \(U(x_{s}^{m})\) is less than or equal to the risk level \(\varepsilon\) then \(x_{s}^{m}\) is a feasible solution with the confidence level \((1-\beta)\). The conclusion in [29] shows that we can obtain the best solution if \(\gamma=\varepsilon/2\). Following [29; 40], the average of the \(L\)th smallest optimal values \(V_{s}^{L}\) obtained in \(S\) iterations can be treated as the lower bound \(L\) of the true optimal value, where \(L\) is calculated as in [29]. Moreover, the true optimal value's upper bound can be estimated by Equation (30). \[\widehat{U}=min_{ \begin{subarray}{c}1\leq m\leq M\\ 1\leq s\leq S\end{subarray}