## 2 Formulation of the Shape Optimization Problem
### Probability Density Evolution Method  
In conservative stochastic systems, fundamental principles analogous to mass, momentum, and energy conservation exist--termed the probability conservation principle [14]. It states that for a stochastic system with time-invariant random factors, probability is preserved during the evolution of system states. Based on this principle, the _generalized probability density evolution equation_ is derived [15] as follows:  
\[\frac{\partial P_{\mathrm{Y}\theta}(\mathbf{y},\mathbf{\theta},\tau)}{\partial\tau}+ \sum_{i=1}^{r}\dot{Y}_{i}(\mathbf{\theta},\tau)\frac{\partial P_{\mathrm{Y}\theta} (\mathbf{y},\mathbf{\theta},\tau)}{\partial y_{i}}=0 \tag{5}\]  
In the formula, \(\tau\) is the system state evolution parameter (which can be regarded as a generalized time parameter), \(Y\) represents the system state (such as the displacement of a key point), \(\theta\) is a basic random variable (such as material yield strength, material _Young's modulus_, etc.), which is the cause of randomness, \(P_{\mathrm{Y}\theta}(y,\theta,\tau)\) is the joint probability distribution density of \((y,\theta)\) at time \(\tau\), and \(r\) is the dimension of the state variable \(Y\). When \(r=1\), only a specific physical quantity is examined (the physical quantity examined in this article is the displacement of the structure), and the above formula degenerates into the following:  
\[\frac{\partial P_{\mathrm{Y}\theta}(\mathbf{y},\mathbf{\theta},\tau)}{\partial\tau}+ \dot{Y}_{i}(\mathbf{\theta},\tau)\frac{\partial P_{\mathrm{Y}\theta}(\mathbf{y},\mathbf{ \theta},\tau)}{\partial y}=0 \tag{6}\]  
For static problems, a proportional loading mechanism is introduced, where the loading parameter (e.g., load magnitude or displacement) serves as \(\tau\).  
Since analytical solutions for most engineering problems are intractable, the _generalized probability density evolution equation_ is typically solved numerically. The procedure involves [16] the following:  
* The value range of the basic random parameter \(\Omega_{\theta}\) is divided, and the \(q\)th region is set as \(\Omega_{q}\). The probability of assigning each part is calculated as \(P_{q}\). In each divided unit, it is defined as follows: \[P_{y}^{(q)}(y,\tau)=\int_{\Omega_{q}}p_{\mathrm{Y}\theta}(y,\mathbf{\theta},\tau)d\mathbf{\theta}\] (7) Then, Equation (7) will become a series of new equations as follows: \[\frac{\partial P_{y}^{(q)}(\mathbf{y},\tau)}{\partial\tau}+\dot{Y}_{i}(\tau)\frac {\partial P_{y}^{(q)}(\mathbf{y},\tau)}{\partial y}=0,q=1,2,3\ldots s\] (8)
* Select representative points \(\mathbf{\theta}_{q}\) within each \(\Omega_{q}\). Solve the physical equations deterministically at these points to obtain \(\dot{Y}\big{(}\mathbf{\theta}_{q},\tau_{m}\big{)}\), where \(\tau_{m}=m\cdot\Delta\tau\) and \(\Delta\tau\) is the time step for solving physical equations.  
* Substitute the physical solutions into Equation (8) with initial conditions and solve numerically. The joint PDF \(p_{\gamma\mathbf{\theta}}\big{(}y_{j},\mathbf{\theta},\tau_{k}\big{)}\) is computed at discrete points \(y_{j}=y_{0}+j\cdot\Delta y\) (\(\Delta y\) is the discrete step of y) and \(\tau_{k}=k\cdot\Delta\overline{\tau}\) (\(\Delta\overline{\tau}\) is the time step for the _probability density evolution equation_), and \(\Delta\tau\) is usually smaller than \(\Delta\overline{\tau}\).
* Integrate the probability density solution \(p_{\gamma\mathbf{\theta}}\big{(}y_{j},\mathbf{\theta},\tau_{k}\big{)}\) obtained in the previous step with respect to \(\mathbf{\theta}\) to obtain the probability density solution \(p_{\gamma}\big{(}y_{j},\tau\big{)}\) of the system response y.  
An open-source code [17] is employed to solve the _generalized probability density evolution equation_, which implements the theoretical framework developed by [PERSON], [PERSON], and colleagues at Tongji University [14; 18; 19; 20; 21].