Optimum Functionally Gradient Materials for Dental Implant Using Simulated Annealing

1. Introduction

Biomaterials should simultaneously satisfy many requirements and possess properties such as non-toxicity, corrosion resistance, thermal conductivity, strength, fatigue durability, biocompatibility and sometimes aesthetics. A single composition with a uniform structure may not satisfy all such requirements. Natural biomaterials often possess the structure of Functionally Graded Materials (FGMs) which enables them to satisfy these requirements. FGMs provide the structure with which synthetic biomaterials should essentially be formed. The size of biomaterial components is relatively small. In the case of dental applications, the components are generally smaller than 20 mm. This substantially reduces the difficulty of fabricating such materials due to a mismatch in thermal expansion which causes micro crack formation during the cooling cycle.

Biomaterials are essential for life and health in certain cases. They have a generally high added value for their size. Thus, biomaterials form one of the most important areas for the application of FGMs. It is an area for which FGMs, at the present time, are sufficiently developed for practical use. The dental implant is used for restoring the function of chewing and biting, and therefore eating, which is the most fundamental activity of human beings required for living. We are living in an era of longer life expectancy and thus, dental care becomes especially important for better quality of life in old age.

Implant may be classified to “implant’’ as an artificial bone for medical use and ‘‘dental implant’’ as an artificial tooth for dental use. The specified properties are slightly different depending on their use. The implants in orthopaedics are used mostly as structurally enforced artificial bone which is inserted inside the corpus. Medical implants lay more weight on strength, toughness, torque in mechanical properties and the specific problem of tribology and abrasion resistance in artificial joint. Dental implant is usually much smaller and is used to reconstruct the masticatory function when tooth root is completely lost or extracted.

Implant is placed in the jawbone in the manner to penetrate from the inside to the outside of the bone. The function is quite different at the inside of bone, outside and at the boundary. In the inside of jaw bone, bone affinity and stress relaxation are important and in the outside of bone, that is, in oral cavity, the sufficient strength is necessary. In the application of human body implant, FGM is usually composed of Collagen Hydroxyapatite (HAP) and titanium (Watari et al., 2004; Hedia, 2005).

HAP is indeed a principal component in human bones and related tissues. HAP inclusion in forming the dental implant material can bring about an enhanced biocompatibility with the native hosting tissues. The main advantages of using FGM dental implant are: 1) reduction of stress shielding effect on the surrounding bones that usually arises in the presence of fully metallic implants (Hedia, 2005), 2) improvement of biocompatibility with bone tissues (Watari et al., 2004), 3) preventing the thermal-mechanical failure at the interface of HAP coated metallic implants (Wang et al., 2007) and 4) meeting the biomechanical requirements at each region of the bone while enhance the bone remodeling, hereby maintaining the bone’s health status (Yang & Xiang, 2007). The latter is more related to volume fraction of FGM. The first three aspects of using FGM implant have been investigated in the previous studies (Watari et al., 2004; Hedia, 2005; Wang et al., 2007).

However, limited knowledge has been available in the effect on bone remodeling due to the use of FGM dental implants. The other issue needs to be systematically studied is how to devise an optimal FGM pattern for dental implant application. It has been widely accepted that a mating mechanical property to the host bone should be made in order to avoid stress shielding (Hedia & Mahmoud, 2004; Hedia et al., 2006) and promote osseointegration and bone remodeling (Chu et al., 2006; Yang & Xiang, 2007). However, there are few reports available which examine whether or not a mating property could result in the best remodeling consequence and ensure a long-term success.

Recently, optimization of FGM dental implant was studied by Lin et al. (2009) using the Response Surface Methodology (RSM) and the results show the incompatibility of properties with each other and the need for using multi-objective algorithms to overcome the problem. Another issue concerns the existing material engineering technology which may not allow us to make such mating pattern for individuals in a cost efficient way. As a result, how to optimally tailor FGM pattern for remodeling is of noteworthy implication in developing FGM implantation.

This chapter aims at extending a more realistic FGM design for dental implantation. Using Simulated Annealing (SA), the multi-objective optimization model was developed to optimize FGM gradient pattern for desirable on-going bone turnover outcome and mechanical responses. SA algorithm has shown great potential for solving optimization problems as they conduct global stochastic search. The multi-objective optimization problem was solved using SA and the results were compared with the RSM.

2. Properties of FGM dental implant

In this study, the configuration of FGM dental implant follows the patterns from literature (Wang et al., 2007; Yang & Xiang, 2007). The material gradient is governed by a power law with parameter m, as in Equation (1). The volume fractions of the two-phase composite FGM dental implant can be calculated from the following equations (Hedia, 2005; Wang et al., 2007; Yang & Xiang, 2007):

$$V_C={\left(\raisebox{1ex}{$y$}\!\left/ \!\raisebox{-1ex}{$h$}\right.\right)}^m$$E1

$$V_m=\left(1-V_C\right)$$E2

where V_c denotes the volumetric fraction of HAP/Col (ceramic), V_m denotes the volumetric fraction of titanium (metal), m is a constant to define the variation in material composition, y is the vertical position within the implant region and h is the total length of the implant. Fig. 1 shows the schematic view of FGM dental implant with graded material composition used in dentistry.

Figure 1.

Schematic view of FGM dental implant with graded material composition.

Accordingly, the Young modulus and Poisson ratio can be calculated as (Hedia, 2005):

$$E_0=E_c\left[\frac{E_c+\left(E_m-E_c\right)v_m^{2/3}}{E_c+\left(E_m-E_c\right)\left(v_m^{2/3}-v_m^{}\right)}\right]$$E3

$$v=v_m{V}_m+v_c{V}_c$$E4

where E₀ is the equivalent Young modulus at different regions of the implant, E_c is the Young modulus of HAP/Col, E_m is the Young modulus of titanium. v_c and v_m are the Poisson ratios for HAP/Col and titanium, respectively. The HAP/Col and titanium compositions vary according to the relative length of y/h, with respect to the material gradient m, meaning that m governs the variation in the volumetric fraction of the titanium to HAP/Col compositions. Referring to the properties of FGM implant, the values of E_c and E_m are kept within the range of E_c >>1 GPa and E_m >> 110 GPa, respectively (Hedia, 2005; Wang et al., 2007; Yang & Xiang, 2007).

Fig. 2 demonstrates the variation of mechanical properties including Young’s modulus and Poisson’s ration for diverse FGM pattern. From Fig. 2, the horizontal axis is the vertical position (y) along FGM dental implant, which is varied from 0 to 10 mm. By observing Fig. 2, y=0 mm indicates the region directly connected to the crown, where FGM has the richest content of titanium when m=10, while the highest content of collagen HAP is obtained when m=0.1. In other words, m=10 and m=0.1, respectively, give the highest and lowest gradients in the Young modulus and Poisson ratio in the region of the crown’s end. Therefore, altering m enables us to tailor the property gradient, thereby providing a means to optimizing the remodeling performance induced by the FGM dental implant.

Figure 2.

Variation in material properties for different FGM configuration: (a) Young’s modulus, (b) Poisson’s ration.

3. Bone remodeling calculations

The biomechanical environment changes considerably when using FGM dental implant. Consequently, the bone remodels itself to adapt to the new changes that is imposed on it by minimizing the difference between the new mechanical response and related equilibrium state. Strain energy density is one of the most important mechanical stimuli to explain the bone remodeling (Weinans et al., 1992). The mathematical equations of bone remodeling are given as follows (Huiskes et al., 1987; Weinans, 1992; Turner et al., 1997; Turner, 1998; Lin et al., 2008a, 2008b):

Bone apposition:

$$\frac{\partial \rho }{\partial t}=B(\frac{\eta }{\rho }-(1+\xi )\Xi )\begin{array}{cc}& if\begin{array}{c}\frac{\eta }{\rho }\end{array}\end{array}\succ (1+\xi )\Xi$$E5

Bone equilibrium:

$$\frac{\partial \rho }{\partial t}=0\begin{array}{cc}& if\begin{array}{c}(1-\xi )\Xi \le \frac{\eta }{\rho }\end{array}\end{array}\le (1+\xi )\Xi$$E6

Bone resorption:

$$\frac{\partial \rho }{\partial t}=B(\frac{\eta }{\rho }-(1-\xi )\Xi )\begin{array}{cc}& if\begin{array}{c}\frac{\eta }{\rho }\end{array}\end{array}\prec (1-\xi )\Xi$$E7

where η denotes the mechanical stimulus (i.e. strain energy density), ρ is the bone density, B is the remodeling constant set to 1gr/cm³ (Weinans et al., 1992), Ξ is the remodeling reference value equal to 0.004 J/g (Weinans et al., 1992; Turner, 1998), and ξ is the bandwidth of bone remodeling with an adapted value of 10% (Weinans et al., 1992). After the bone density values are calculated from the remodeling equations, the Young modulus of cortical and cancellous bones (in GPa) can be updated by using the following equations (Rho et al., 1995; O’Mahony et al., 2001):

$$E_{cortical}=-23.93+24\rho {}_{cortical}$$E8

$$E_{cancellous}=2.349{\rho }_{cancellous}^{2.15}$$E9

Equations (8) and (9) are utilized to update Young modulus after the densities are determined via the remodeling calculations. The internal bone remodeling system is formed using Equation (5) to Equation (9).

4. Design optimization problem

The bone remodeling provides quantitative data of changes in bone densities and the stiffness of dental apparatus. The former indicates how the bones react to the change in biomechanical environment in terms of the variation in bone morphology. The latter indicates how the bone remodeling alters the mechanical response, thereby stabilizing the implant and in turn strengthening the bone. In this research, the changes in bone densities and vertical displacement are taken as the direct measures of on-going performance of implantation.

From the biomechanical point of view, increase in surrounding bone density and decrease in the occlusal displacement indicate the positive sign to a long-term success in dental implantation. Thus, design of FGM gradient parameter (m) is expected to maximize the densities and minimize the displacement, which in a form of multi-objective optimization may be represented as:

$$\{\begin{array}{l}\min f_1\left(m\right)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$D_{cortical}$}\right.\hfill \\ \min f_2\left(m\right)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$D_{cancellous}$}\right.\hfill \\ \min f_3\left(m\right)=u\left(m\right)\hfill \\ subjectto0.1\le m\le 10\hfill \end{array}.$$E10

where f₁, f₂ and f₃ represent the objective functions, D_cortical and D_cancellous are the densities of cortical and cancellous bones, respectively and u(m) denotes the vertical displacement at the top of artificial crown. The objective functions f₁, f₂ and f₃ represent the condition of FGM dental implant at month 48 (Lin et al., 2009). These polynomial response functions are given as:

$$D_{cortical}\left(m\right)=-2e^{-6}{m}^6+8e^{-5}{m}^5-1e^{-3}{m}^4+6.2e^{-3}{m}^3-1.9e^{-2}{m}^2+1.76e^{-2}m+1.9297$$E11

$$D_{cancellous}\left(m\right)=2e^{-6}{m}^6-6e^{-5}{m}^5+6e^{-4}{m}^4-2.6e^{-3}{m}^3+3e^{-2}{m}^2-2.4e^{-3}m+1.1712$$E12

$$u\left(m\right)=7e^{-10}{m}^6-2e^{-8}{m}^5-3e^{-7}{m}^4-2e^{-6}{m}^3+7e^{-6}{m}^2-1e^{-5}m+4e^{-5}$$E13

The polynomial response functions were obtained using experimental tests on various quantities of material gradient on the FGM dental implant. Lin et al. (2009) proposed to adopt the RSM to construct the appropriate objective functions for a single design variable problem. Fig. 3 represents the cortical and cancellous densities, and displacement versus material gradient (m). As m increases, the displacement function u(m) is increased as shown in Fig. 3 c. The increase in m results in an increase in the cortical density as shown in Fig. 3 a. From Fig. 3b, the increase in m results the decrease in the density of cancellous.

The design objectives were to maximize both D_cortical and D_cancellous in order to determine the best possible material configuration for implant that will give the maximum amount of bone remodeling. At the same time, the objective function f₃ is minimized in order to reduce the downward implant displacement. An ideal situation would be to attain a consistent optimal material gradient, where the maximum bone turnover can be made and the displacement is kept to minimal. To explore such a multi-objective design, first, the two objective optimization problems of either cortical or cancellous densities versus the displacement were formulated as follows:

$$\{\begin{array}{l}\min f_1=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$D_{cortical}$}\right.\hfill \\ \min f_3=\max \left(u\left(m\right)\right)\hfill \\ subjectto0.1\le m\le 10\hfill \end{array}$$E14

and

$$\{\begin{array}{l}\min f_2=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$D_{cancellous}$}\right.\hfill \\ \min f_3=\max \left(u\left(m\right)\right)\hfill \\ subjectto0.1\le m\le 10\hfill \end{array}.$$E15

Figure 3.

Mass densities and displacement for FGM dental implant: (a) cortical density, (b) cancellous density, (c) displacement (u) (The horizontal axis is m).

Fig. 4 shows the trend of cortical density function (f₁), cancellous density function (f₂) and the displacement function (f₃). From Figs. 4b and 4c, it is clear that f₂ and f₃ have the same trend and behavior and it is expected to obtain one optimal solution. Based on Fig. 3c, we need to minimize the maximum of displacement in order to obtain the minimum value of displacement. The function f₃ is minimized by multiplying u(m) by -1. In other hand, Figs. 4a and 4c show the different trend and behavior which with increasing m, the cortical density function decreases while, the displacement function increase. Therefore, it is anticipated to obtain a range of optimal solutions. The problem investigated in this chapter, is taken from (Lin et al., 2009) where the results were obtained using RSM.

Figure 4.

Mass densities and displacement functions for FGM dental implant: (a) cortical density function (f₁), (b) cancellous density function (f₂), (c) displacement function (f₃) (The horizontal axis is m).

5. Multi-objective optimization

The multi-objective optimization has become an important research topic for scientists and researchers. This is mainly due to the multi-objective nature of real life problems. It is difficult to compare results of multi-objective methods to single objective techniques, as there is not a unique optimum in multi-objective optimization as in single objective optimization. Therefore, the best solution in multi-objective terms may need to be decided by the decision maker.

The increasing acceptance of SA for solving multi-objective optimization problems is due to their ability to: 1) find multiple solutions in a single run, 2) work without derivatives, 3) converge speedily to Pareto-optimal solutions with a high degree of accuracy, 4) handle both continuous functions and combinatorial optimization problems with ease and 5) are less susceptible to the shape or continuity of the Pareto front (Suman & Kumar, 2006).

5.1. Pareto-optimal solutions

The concept of the Pareto-optimal solutions was formulated by Vilfredo Pareto in the 19^th century (Rouge, 1896). Real life problems require simultaneous optimization of several incommensurable and often conflicting objectives. Usually, there is no single optimal solution. However, there may be a set of alternative solutions. These solutions are optimal in the wider sense that no other solutions in the search space are superior to each other when all the objectives are considered. They are known as Pareto-optimal solutions. When the objectives associated with any pair of non-dominated solutions are compared, it is found that each solution is superior with respect to at least one objective. The set of non-dominated solutions to a multi-objective optimization problem is known as the Pareto-optimal set (Zitzler & Thiele, 1998).

6. Simulated annealing

In 1953, Metropolis developed a method for solving optimization problems that mimics the way thermodynamic systems go from one energy level to another (Metropolis et al., 1953). He thought of this after simulating a heat bath on certain chemicals. This method is called Simulated Annealing (SA). Kirkpatrick et al. (1983) originally thought of using SA on a number of problems. The name and inspiration come from annealing in metallurgy, a technique involving heating and controlled cooling of a material to increase the size of its crystals and reduce their defects. The heat causes the atoms to become free from their initial positions (a local minimum of the internal energy) and wander randomly through states of higher energy.

The system is cooled and as the temperature is reduced the atoms migrate to more ordered states with lower energy. The final degree of order depends on the temperature cooling rate. The slow cooling process is characterized by a general decrease in the energy level for with occasional increase in energy. On the other hand, a fast cooling process, known as quenching, is characterized by a monotonic decrease in energy to an intermediate state of semi-order which is used as temperature schedule in this chapter.

At the final stages of the annealing process, the system’s energy reaches a much lower level than in rapid cooling (quenching). Annealing (slow cooling) therefore allows the system to reach lower global energy minimum than is possible using a quick quenching process, equivalent to a local energy minimum.

By analogy with this physical process, each step of the SA algorithm replaces the current solution by a random "nearby" solution, chosen with a probability that depends both on the difference between the corresponding function values and also on a global parameter T (temperature), that is gradually decreased during the process. The dependency is such that the current solution changes almost randomly when T is large, but increasingly "downhill" as T goes to zero (Fleischer, 1995). The allowance for "uphill" moves potentially saves the method from becoming stuck at local optima. Several parameters need to be included in an implementation of SA.

These are summarized by Davidson and Harel (1996):

• The set of configurations, or states, of the system, including an initial configuration (which is often chosen at random).

• A generation rule for new configurations, which is usually obtained by defining the neighborhood of each configuration and choosing the next configuration randomly from the neighborhood of the current one.

• The target, or cost function, to be minimized over the configuration space. (This is the analogue of the energy).

• The cooling schedule of the control parameter, including initial values and rules for when and how to change it (This is the analogue of the temperature and its decreases).

• The termination condition, which is usually based on the time and the values of the cost function and/or the control parameter.

SA is a popular optimization algorithm due to the simplicity of the model and its implementation. However, due to CPU time-consuming nature of standard SA, a fast temperature schedule to fulfill the required conditions is suggested. In fact, simulated annealing algorithm with the fast cooling process is called simulated quenching (SQ) which is used as an optimization method in this chapter to overcome the slow SA optimization process.

7. SA for multi-objective optimization

SA has been used as an optimization method for solving a wide range of combinatorial optimization problems. It has also been adapted for solving multi-objective problems due to its simplicity and capability of producing a desirable Pareto set of solutions. In addition, it is not susceptible to the shape of the Pareto set, since shape may be considered as a concern for mathematical programming techniques.

7.1. The method of Suppapitnarm and Parks (SMOSA)

The concept of archiving the Pareto-optimal solutions for solving multi-objective problems with SA has been used by Suppapitnarm et al. (2000). The Suppapitnarm Multi-Objective Simulated Annealing (SMOSA) enables the search to restart from an archived solution in a solution region, where each of the pair of non-dominated solutions may be superior with respect to at least one objective. Since SA only generates a single solution at a given iteration, an independent archive is required to record all non-dominated solutions found during search. Pioneering work in this area was first performed by Engrand (1997), and was further developed by Suppapitnarm et al. (2000).

7.1.1. Acceptance probability and archiving

A new acceptance probability formulation based on an annealing schedule with multiple temperatures (one for each objective function) was also proposed. The changes in each objective function values are compared with each other directly before archiving. This ensures that the moves to a non-dominated solution are accepted. It does not use any weight vector in the acceptance criteria. Hence, the acceptance probability step is given as:

$$P=\min \left(1,{\displaystyle \prod _{i=1}^N\exp \left\{\raisebox{1ex}{$-\Delta s_i$}\!\left/ \!\raisebox{-1ex}{$T_i$}\right.\right\}}\right)$$E16

where ΔS=(Z_i(Y)-Z_i(X)) and N is the number of objective functions, X is the current solution, Y is the generated solution, Z_i is the objective function, and T_i is the annealing temperature. Thus, the overall acceptance probability is the product of individual acceptance probabilities for each objective associated with a temperature T_i.

7.1.2. Annealing schedule

A new annealing schedule is developed to control the lowering of individual temperatures associated with each objective function. If the temperatures are lowered too fast the chance of accepting solutions reduces rapidly and large parts of the search space are never explored. In contrast, if the temperature is reduced too slowly, many redundant solutions which do not lead to non-dominated solutions are accepted and the Pareto-optimal set of solutions develops very slowly. The latter is particularly undesirable if objective function evaluations are expensive and/or if computation time is a important factor.

A statistical record of the values of each of the objective functions (f_i) is maintained. First, the temperatures are lowered after N_T1 iterations by setting each temperature to the standard deviation (σ) of the accepted values of f_i (T_i = σ_i). Thereafter, the temperatures based on the quenching schedule are updated after every N_T2 iterations or N_A acceptances as follows:

$$T_{i(k+1)}=T_{i(k)}\times \alpha \begin{array}{cc}& 0<{\alpha }_i<1\end{array}$$E17

where T_i is the temperature, k is the time index of annealing, and α_i is the cooling ratio of each objective function. The suitable values for N_T1 and N_T2 were chosen 1000 and 500 iterations, respectively (Suppapitnann, 1998).

7.1.3. Return to base strategy

In order to completely expose the trade-off between objective functions, the random selection of a solution from the archive, from which to recommence search, is systematically controlled using an intelligent return-to-base strategy. After the start of search, a return-to-base is first activated when the basic features of the trade-off between objectives have developed. It seems sensible that this take place when the temperatures are first lowered, i.e., after N_T1 iterations. Thereafter, the rate of return is, naturally, increased to intensify the search in the trade-off. The number of iterations N_Bi to be executed prior to the i^th return-to-base after the start of search is updated as given:

$$N_{Bi}=r_B{N}_{Bi-1}\begin{array}{cc}& i=2,3,4,\mathrm{...}\end{array}$$E18

where r_B is a constant parameter which varies between 0 and +1 and dictates the frequency of return. Recommendation values for r_B and N_B1 may be chosen as 0.9 and 2N_T2, respectively (Suppapitnann, 1998). In order to fully develop the trade-off, solutions that are more isolated from the rest of the trade-off solution should be favored in returns-to-base. The extreme solutions, those solutions that correspond to minimum values for each objective in the trade-off, also require special consideration. These solutions are almost invariably only just feasible, which makes the design space around them difficult to search.

For these reasons, a base set of candidate solutions was proposed which consists of a number of the most isolated of those solutions currently held in the archive and the M extreme solutions in the archive. Therefore, when a return-to-base is activated, the search diversifies into less well explored regions of the trade-off. To evaluate the degree of isolation for a solution, the following formula was proposed (Suppapitnarm et al., 2000):

$$I(X_j)={{\displaystyle \sum _{\begin{array}{l}i=1\\ i\ne j\end{array}}^{A_s}{\displaystyle \sum _{k=1}^M\left(\frac{\left|f_k(X_i)-f_k(X_j)\right|}{\left|f_{k\max }-f_{k\min }\right|}\right)}}}^2$$E19

where I(X_j) is the normalized value for distance in objective space for the j^th solution from all other archived solutions and X_j denotes the j^th archived solution. A_s and M are the total number of solutions and extreme solutions stored in the archive, respectively. f_kmax and f_kmin are the maximum and minimum values for k^th objective function (f_k), respectively. Each solution - except for the extreme solutions - is ranked in order to decrease isolation distance, thereby, establishing an ordered set with the most isolated solutions at its top and the least isolated solutions at the bottom.

7.1.4. Step size control

An improvement in SA performance may be gained by varying the maximum allowable step changes in each of the decision variables during perturbation between iterations (Parks, 1990). Hence, the value of each design variable is rescaled to U_ik such that it varies between - 1 and +1 at its lower and upper bounds, respectively. At the next iteration, U_i(k+1) is modified as given:

$$U_{i(k+1)}=U_{ik}+rand\times S_i$$E20

where rand is a uniformly distributed random number between -1 and +1, and S_i is the maximum (positive) step-size for each design variable. If the solution is accepted, S_i is updated using following equation:

$$S_i=S_i(0.9+0.21\times \left|rand\right|)$$E21

A suitable value for the upper bound for each S_i was set to 0.5 to permit, initially, a wide search around the current position. Accordingly, a lower bound of 0.0001 was chosen for the smallest possible value of each S_i to prevent stagnation during search (Parks, 1990). Therefore, the maximum step change in the design variables is monitored and is varied to reduce violation of the constraints.

7.1.5. The steps and flowchart of the SMOSA

The basic steps involved in the SMOSA algorithm for a problem having N objective functions and n decision variables are as follows (Suman, 2004):

1. Start with a randomly generated initial solution vector, X (an n×1 vector whose elements are decision variables) and evaluate all objective functions and put it into a Pareto set of solutions.

2. Generate a random perturbation and a new solution vector, Y, in the neighborhood of current solution vector, X, re-evaluate the objective functions and apply a penalty function approach to the corresponding objective functions, if necessary.

3. Compare the generated solution vector with all the solutions in the Pareto set and update the Pareto set, if necessary.

4. If the generated solution vector is archived, assign it as the current solution vector by putting X=Y and go to Step 7.

5. If the generated solution vector is not archived, accept it with the probability using Equation (16). If the generated solution is archived, assign it as the current solution vector by putting X=Y and go to Step 7.

6. If the generated solution vector is not archived, retain the earlier solution vector as the current solution vector and go to Step 7.

7. Periodically, restart with a randomly selected solution from the Pareto set. While periodically restarting with the archived solutions, Suppapitnarm et al. (2000) have recommended biasing towards the extreme ends of the trade-off surface.

8. Reduce the temperature using Equation (17) and annealing schedule.

9. Repeat Steps 2 to 8, until a predefined number of iterations is carried out.

In addition, the flowchart of SMOSA optimizer is illustrated in Fig. 5.

8. Optimization of the two objective functions using SMOSA

An ideal situation would be to attain a consistent optimal material gradient, where the maximum one turnover can be made and the displacement is kept to minimal. The problem described in previous sections was solved by multi-objective SA code written in MATLAB software programming and run on Pentium IV, 2500 GHz and 4 GB RAM. The SMOSA was run 5 times and the obtained averaged results were compared to results obtained from RSM.

The parameters which must be specified before running the algorithm are initial temperature, frozen state represented by the final temperature, cooling ratio (annealing), i.e. the rate at which the temperature is lowered between two cooling cycles, the randomly generated initial solution, and the lower and upper bounds for design variable (m). Table 1 shows the user parameters used for the SMOSA.

Figure 5.

Flowchart of the SMOSA optimizer.

User Parameters	Quantities
Initial temperature for each objective (Ti0)	100
Final temperature for f2 and f3	1e-6
Final temperature for f1 and f3	1e-8
Cooling ratio (αi)	0.98
Lower bound	0.1
Upper bound	10

Table 1.

Optimization parameters associated with the SMOSA.

The SMOSA was conducted on two objective functions, f₃ as displacement function and f₁ as cortical density function as in Equation (14). Fig. 6 represents the comparison of Pareto front for two objective functions using SMOSA and RSM (Lin et al., 2009). Comparing Figs. 6a and 6b, the same trend is seen in both figures. It can be seen that the increase in f₁, results in f₃ to decrease and vice versa. Hence, SMOSA confirms the results obtained by RSM. However, the data ranges obtained by SMOSA do not match exactly with the data range obtained by RSM.

The data range given by SMOSA for f₁ is from 0.51688 to 0.51778 as shown in Fig. 6b while, RSM gives values from 0.5169 to 0.5175 for f₁ (see Fig. 6a). Similarly, the data range given by SMOSA for f₃ is from -3.9054e-5 to -3.5907e-5 and for RSM for f₃ is from 3.5205e-5 to 3.603e-5. The date range for cortical density function given by SMOSA is almost 33% more than the range obtained by RSM in micro scale. This shows the efficiency of the SMOSA in providing more data range for cortical density than RSM.

The safety factor (acceptance capability) obtained by SMOSA is 33% higher than the one given by RSM for magnitude of cortical density function. This increasing of data range is interesting for FGM dental implant design and shows that SMOSA has outperformed RSM in cortical density function. SMOSA gives the data range for displacement almost 73% more than the values given by RSM in micro scale. That means SMOSA enables the designer with a wider choice for design.

In general, SMOSA method gives more selection of material gradient (m) for designing of the FGM dental implant but with higher displacement compared to the RSM. The order of quantities (micro-scale) shows the importance of accuracy in the optimization of the FGM implant. Hence, any improvement in quantities may be considerable although such improvements may seem to be negligible.

Fig. 7 illustrates the trend for cortical density and displacement functions with respect to m. By inspecting Fig. 7, the acceptable range for m is from 0.1 to 0.65. The magnitude of f₁ decreases as m increases which means the density of cortical (D_cortical) increases. In contrast, with the increase in m the quantity of vertical displacement (f₃) increases. This trend shows the non-dominated solution from the Pareto front plot so that the increase in one function results in the decrease of the other function and vice versa. Such data was extracted from the Pareto front depicted in Fig. 6 b.

Figure 6.

Pareto front for optimization of two objective functions (f₃ and f₁) using: (a) RSM, (b) SMOSA.

Fig. 8 demonstrates f₁, f₂ and f₃ with respect to the m. When the temperature is high at the beginning of the optimization process, all invalid moves were accepted by the acceptance probability (P), as shown in Fig. 8. By decreasing the temperature at each iteration, the probability of accepting invalid moves is reduced and therefore, only qualified points and valid moves (so called non-dominated solutions) were accepted.

As a next step, the multi-objective simulated annealing method was applied on other two objective functions as shown in Equation (15). The result indicates the Pareto front is a cluster of points that are gathered in one point representing an optimal solution. Going back to Figs. 4b and 4c for minimizing both objective functions (f₂ and f₃), we expect to have one optimal solution as the two functions have a similar trend to reach the optimal point.

In other words, if one moves from the right side of the Figs. 4 b and 4c to the left side, the cancellous density increases while the displacement decreases. This situation may be considered as an optimal state. The optimal material gradient is taken as 0.1 (m=0.1). The values for f₃ and f₂ are equal to -3.9068e-5 and 0.8540, respectively. Hence, the best material gradient for these two objective functions is given when m=0.1.

Figure 7.

Cortical density and displacement functions versus m using SMOSA for: (a) f₁, (b) f₃.

9. Optimization of the three objective functions using SMOSA

The results obtained from the optimization of two objective functions may not adequately address the full design requirements and therefore, the three objective functions were developed as defined in Equation (10). Fig. 9a shows the Pareto front surface for three objective functions using RSM. Fig. 9b illustrates the 3D plot of the Pareto front for three objective functions (f₁, f₂, and f₃) using SMOSA. From Fig. 9, the good harmony and similarity between the results obtained by SMOSA and RSM are depicted.

Fig. 10 represents various aspects of Fig. 9 b. This figure shows the trend for f₁, f₂, and f₃ in terms of the material gradient. The acceptable range for m varies from 0.1 to 0.65 (see Fig. 10). As shown in Fig. 10a, the maximum value for cortical density is obtained for m=0.63. In Fig. 10b, the increase in f₂ leads to increase in m, and hence, decrease in the cancellous density. This indicates that for material gradient 0.1 (m=0.1) cancellous density has the maximum value. In addition, by observing Fig. 10c, when m=0.1, the displacement has the minimum value. Finally, based upon the obtained results, the optimal range for material gradient is almost varying from 0.1 to 0.65 (0.1 ≤ m ≤ 0.65).

Figure 8.

Distribution of mass densities and displacement functions with respect to the m using SMOSA: (a) f₁, (b) f₂, (c) f₃.

Figure 9.

Pareto front for three objective functions (f₁, f₂ and f₃) using: (a) RSM, (b) SMOSA.

Figure 10.

Mass densities and displacement functions in optimization of three objectives in terms of material gradient using SMOSA for: (a) f₁, (b) f₂, (c) f₃.

10. Conclusions

This chapter presents some important findings for Functionally Graded Material (FGM) dental implant design in a power law configuration. This is vital for maintaining the overall health of the bone tissues. The research clearly suggests that, a better performance in bone turnover can be achieved by lowering the FGM dental implant material gradient. However, this will at the same time reduce the stiffness of implantation, consequently placing the bone-implant interface at higher risk of damage during the early healing stage.

The problem may be solved by the multi-objective optimization method. The Pareto front was determined using the Suppapitnarm Multi-Objective Simulated Annealing (SMOSA) optimization procedure. The results obtained from the SMOSA confirm the results obtained by the Response Surface Methodology (RSM), in addition to offering further improvements. The SMOSA optimized the objective functions on a wider data range than RSM and offered better results with respect to the cortical density function (almost 33% more than RSM). SMOSA optimization in this case, gives more selection of material gradient (m) for designing FGM dental implant compared to RSM. The material gradient varies from 0.1 to 0.65 given by SMOSA.

By considering the point that the scale in the FGM dental implant is in micro, the importance of accuracy in optimization of the FGM implant is understood. The design of FGM gradient parameter is expected to maximize the densities (cortical and cancellous) and minimize the displacement and plays a more important role in the design methodology. However, sacrifice may be made when the third criterion of displacement is introduced, which means that an optimal gradient m for bone remodeling may not be the best for stiffness. It is expected that the design methodology can produce more favorably patient specific implant, better improving the immediate and long-term restorative outcomes.

Acknowledgement

The authors would like to acknowledge for the Ministry of Higher Education of Malaysia and the University of Malaya, Kuala Lumpur, Malaysia for the financial support under UM.TNC2/IPPP/UPGP/628/6/ER013/2011A grant.