## 3 Topology Optimisation with a Bi-Directional Bone Remodelling Algorithm
### Bi-Directional Evolutionary Structural Optimisation  
Evolutionary computation, inspired by biological evolution, utilises selection, reproduction, and variation mechanisms to find optimal solutions for complex problems [24]. Compared with traditional optimisation techniques, evolutionary methods offer greater robustness, exploration capabilities, and flexibility. These attributes make them well-suited for tackling intricate problems with challenging cost functions. In recent years, methods like Evolutionary Structural Optimisation (ESO), developed by [PERSON] and [PERSON] [25], have gained widespread application in structural optimisation problems. This iterative technique removes material deemed inefficient or redundant from a specific domain to achieve an optimal design. While the ESO method demonstrates widespread application and popularity, it does possess limitations that have motivated the exploration of alternative techniques [26]. To enhance the solution viability achieved through optimisation, the concept of a bidirectional algorithm emerged. This approach would not only remove material in low-stress areas but also introduce material to reinforce high-stress regions. This led to the development of the Bi-Directional Evolutionary Structural Optimisation (BESO) method, drawing inspiration from both ESO's material removal capabilities and the additive material functionalities of the Additive Evolutionary Structural Optimisation (AESO) method [27]. BESO facilitates a more comprehensive exploration of the design domain, potentially leading to superior identification of the global minimum. The elasto-static analysis step initiates each optimisation iteration and returns the displacement, strain, and stress fields. This allows for the calculation of the equivalent von Mises stress for each point, along with the cubic average of the entire von Mises stress field. This average serves as a reference to identify potential locations of abrupt stress variations. High stress values influence the cubic average, while low stress values have minimal impact. \[\overline{\sigma}(x_{l})=\sqrt{\frac{(\sigma_{xx}-\sigma_{yy})^{2}+(\sigma_{yy} -\sigma_{zz})^{2}+(\sigma_{zz}-\sigma_{xx})^{2}+6(\tau_{xy}^{2}+\tau_{yz}^{2}+ \tau_{xx}^{2})}{2}} \tag{27}\]  
\[\overline{\sigma}_{cube}=\sqrt[3]{\frac{1}{n_{Q}}\sum_{i=1}^{n_{Q}}\overline {\sigma}(\mathbf{x}_{l})^{3}} \tag{28}\]  
A penalty system is implemented to designate a specific parameter, mass density in this case, to each integration point. The penalty values typically range from \(\xi\in[10^{-3},1]\), where 1 signifies rewarded domains (solid material) and \(10^{-3}\) signifies penalised domains (removed material). The performance of the BESO procedure hinges on the reward ratio, RR, and the penalty ratio, PR. The procedure identifies the \(n_{Q}^{RR}=RR\cdot n_{Q}\) integration points with the highest \(\sigma(\mathbf{x}_{l})\) (stress) and the \(n_{Q}^{PR}=PR\cdot n_{Q}\) integration points with the lowest \(\sigma(\mathbf{x}_{l})\). The \(n_{Q}^{RR}\) points are assigned a reward of \(\xi(\mathbf{x}_{l})=1\), while the \(n_{Q}^{PR}\) points receive a penalty of \(\xi(\mathbf{x}_{l})=10^{-3}\). A penalty parameter, \(\eta(\mathbf{x}_{i})\), is then assigned to each node \(\mathbf{x}_{i}\). For each integration point \(\mathbf{x}_{l}\), the closest nodes update their penalty values using the expression \(\eta(\mathbf{x}_{i})^{new}=0.5(\eta(\mathbf{x}_{i})^{old}+\xi(\mathbf{x}_{l}))\). Once all nodes are updated, the penalty parameters for each integration point are recalculated to filter and smooth the selected values. This is achieved using the interpolation function \(\xi(\mathbf{x}_{l})=\sum_{i=1}^{n}\varphi_{i}(\mathbf{x}_{l})\cdot\eta( \mathbf{x}_{i})\), where \(n\) represents the number of nodes within the analysed influence domain of \(\mathbf{x}_{l}\), and \(\varphi(\mathbf{x}_{l})\) represents the shape function vector of the integration point \(\mathbf{x}_{l}\). After the first iteration, some \(\xi\) values might deviate from one, indicating the absence of material.