### Online Optimization on the Macroscale  
In this subsection, online optimization was conducted to obtain the relative density distribution of the design domain. The optimization formulation and sensitivity analysis are given and calculated for the heat transfer problem. In the optimization process, an element's relative density was taken as the design variable:  
\[(\mathbb{P}) \begin{cases}\min_{\rho}&J_{\text{com}}=\mathbf{F}^{\text{T}}\mathbf{T}( \rho)\\ s.t.&\mathbf{K}\mathbf{T}=\mathbf{F}\\ &V(\rho)-\overline{V}\leq 0\\ &0<\rho_{\text{min}}\leq\rho\leq\rho_{\text{max}}\leq 1,\end{cases}  \tag{9}\]  
where \(J\) is the objection function, \(\rho\) is the relative density, \(\mathbf{F}\) is the external heat load vector, \(\mathbf{T}\) is the temperature vector, \(\mathbf{K}\) is the global heat conductivity matrix, \(V=\int_{\Omega}\rho d\Omega\) is the current material volume fraction, \(\overline{V}\) is the predefined material volume fraction, and \(\rho_{\text{min}}\) and \(\rho_{\text{max}}\) are the lower and upper limit of the relative density, respectively. \begin{table}
\begin{tabular}{c c c c} \hline \hline \(\mathbf{a}\) & \(\mathbf{b}\) & \(\mathbf{c}\) & \(\mathbf{d}\) \\ \hline
43.67 & \(-13.42\) & \(10.08\) & \(-2.11\) \\ \hline \hline \end{tabular}
\end{table} 
Table 2: The fitting parameters \(a,b,c\), and \(d\) in Equation (8). Figure 4: Mapping relationship: (**a**) the relative density and the average aspect ratio; (**b**) the relative density and the postprocessing equivalent heat conductivity. The OC optimality criteria [38] were employed in this work to solve the heat transfer problem:  
\[\rho_{e}^{\text{new}}= \begin{cases}\max(\rho_{\text{min}},\rho_{e}-\Delta)&\text{ if }\rho_{e}\Lambda_{e}^{\zeta}\leq\max(\rho_{\text{min}},\rho_{e}-\Delta)\\ \min(\rho_{\text{max}},\rho_{e}+\Delta)&\text{if }\rho_{e}\Lambda_{e}^{\zeta}\geq\min( \rho_{\text{max}},\rho_{e}-\Delta)\\ \rho_{e}\Lambda_{e}^{\zeta}&\text{otherwise}\end{cases} , \tag{10}\]  
where \(\Delta\) is a positive move limit, \(\zeta\) is a damping coefficient, and \(\Lambda_{e}\) is calculated by  
\[\Lambda_{e}=\frac{-\frac{\partial J}{\partial\rho_{e}}}{\lambda\frac{\partial V }{\partial\rho_{e}}}, \tag{11}\]  
where \(\lambda\) is the Lagrange multiplier. It is stressed that \(\rho_{\text{min}}\) should be a positive value in the optimization for numerical stability. As can be seen in the mapping relationship, when \(\rho_{\text{min}}\) was close to zero, the corresponding average aspect ratio would be relatively large and the elastic constants were relatively small. As a result, the optimization process would go through oscillation, which made the reconstruction of the Voronoi-GSCS very difficult. Consequently, in this paper, \(\rho_{\text{min}}=0.1\) was used. The finite element method was used to make the design domain discrete and to solve the heat equilibrium equation. The derivative of the objective function concerning an element's relative density is calculated by  
\[\frac{\partial J}{\partial\rho_{e}}=-T^{\text{T}}\frac{\partial\mathbf{K}}{ \partial\rho_{e}}T, \tag{12}\]  
where the derivative of the global heat conductivity matrix concerning an element's relative density is calculated by  
\[\frac{\partial\mathbf{K}}{\partial\rho_{e}}=\sum_{e=1}^{N_{e}}\mathbf{C}_{e}^ {\text{T}}\Bigg{(}\int_{Y_{e}}\mathbf{B}_{e}^{\text{T}}\frac{\partial\mathbf{ d}_{\text{post}}^{\text{H}}}{\partial\rho_{e}}\mathbf{B}_{e}\mathrm{d}Y \Bigg{)}\mathbf{C}_{e}, \tag{13}\]  
where \(N_{e}\) is the total number of elements, \(\mathbf{C}_{e}\) is the expansion matrix, \(\mathbf{B}\) is the strain matrix, and \(\frac{\partial\mathbf{d}_{\text{post}}^{\text{H}}}{\partial\rho_{e}}\) can be calculated by,  
\[\frac{\partial\mathbf{d}_{\text{post}}^{\text{H}}}{\partial\rho_{e}}=\left[ \begin{array}{cc}\frac{\partial(d_{xx}+d_{yy})}{2\partial\rho_{e}}&0\\ 0&\frac{\partial(d_{xx}+d_{yy})}{2\partial\rho_{e}}\end{array}