- Email: [email protected]
Design window search using continuous evolutionary algorithm and clustering––its application to shape design of microelectrostatic actuator
Computers and Structures 80 (2002) 2469–2481 www.elsevier.com/locate/compstruc
Design window search using continuous evolutionary algorithm and clustering––its application to shape design of microelectrostatic actuator Daisuke Ishihara a, Min Joong Jeong b, Shinobu Yoshimura
b,*
, Genki Yagawa
a
a
b
School of Engineering, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Department of Quantum Engineering and Systems Science, Institute of Environmental Studies, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Received 31 May 2001; accepted 25 July 2002
Abstract This paper describes a search method for multi-dimensional design window (DW), which is defined as an existing area of satisfactory solutions in a design parameter space. The method consists of the following two steps: (1) direct search for satisfactory design solutions using a continuous evolutionary algorithm (CEA) from a wide area of the design parameter space in a robust manner, and (2) identification of a precise structure of DW by clustering the detected satisfactory solutions with a modified K-means algorithm. The CEA is a kind of genetic algorithms modified to deal with continuous variables. The modified K-means clustering algorithm contains an explicit procedure to determine the optimum number of clusters. The proposed DW search method was implemented to an integrated computational aided engineering system for multi-disciplinary structural design developed by the authors, and then the method was applied to shape design of a microelectrostatic actuator for next generation high-density optical memory. A DW consisting of four clusters, i.e. four sub-DWs was obtained, and the features of the representative design solutions of the four subDWs were compared with each other in detail, and a final design solution was determined. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Structural design; Design window; Continuous evolutionary algorithm; K-means clustering; Finite element analysis; Microelectrostatic actuator
1. Introduction In a field of structural design, much research efforts have been put onto the development of an effective procedure to find an optimum solution, assuming that an objective function and all constraints are explicitly formulated in mathematical forms. In practical engineering situations, it is however natural that some of design requirements are hardly formulated in a mathe-
*
Corresponding author. Tel.: +81-3-3812-6960; fax: +81-35800-6876. E-mail address: [email protected] (S. Yoshimura).
matical form. Those unformulated requirements are taken into account in the last stage of design. Manufacturability and manufacturing cost may be some of such requirements. In these cases, it seems more meaningful to obtain satisfactory solutions that meet all the formulated design requirements rather than obtaining a single optimum solution. Thus we focus on a design window (DW) search method in this paper. The DW schematically indicates an existing area of satisfactory solutions within a design parameter space. Some DW search methods have been proposed and applied in several kinds of engineering problems [1–5]. In Ref. [1], a DW projected to a two-dimensional space is detected using a few points search method. In Ref. [2], a DW is detected using a simplified mathematical
0045-7949/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 2 9 3 - 6
2470
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
model approximating an original design problem. These approaches are too simple to apply to such general design problems that a relation between design parameters and design characteristics is nonlinear, discontinuous, multi-dimensional and multi-modal. Some of the present authors also studied a DW search method [3–5]. In Ref. [3], a whole area search method with a help of a multi-layer neural network was proposed, and applied to some engineering applications [4,5]. Here an orthogonal lattice is first generated in a design parameter space, and each lattice point is evaluated whether to be a satisfactory solution or not. The multi-layer neural network was used to approximate an actual relation between design parameters and design characteristics for the purpose of shortening DW searching time. Being similar to response surface methods, such an approximate method might introduce unallowable errors when solving highly nonlinear and discontinuous design problems [6–8]. Ref. [9] reported that errors were larger than 50% in some design problems. In this study, we newly propose a direct DW search method that is applicable to general design problems with nonlinear, discontinuous, multi-dimensional and/or multi-modal features. This is a two-steps method using a continuous evolutionary algorithm (CEA) and a modified K-means algorithm. At the first step, using the CEA, satisfactory solutions are searched out of a whole area of the design parameter space in a robust matter. The CEA is a kind of genetic algorithms (GAs) modified to deal with continuous variables and is good at solving illposed design problems [10]. After a sufficient number of evaluations, it is well expected that all the obtained satisfactory solutions may be a discrete expression of a continuous DW. Evolutionary algorithms (EAs), which have been named to represent GA [11], evolutionary programming [12], evolution strategy [13,14] and their combined algorithm [15], have appeared as robust optimization techniques in a last few decades. There are many studies on engineering applications of the EAs. The CEA, which was developed by Furukawa and Yagawa [10] is one of the GAs for real search space. In the past researches, the GAs have been employed to find a unique optimum solution. In the present study, the CEA is, on the other hand, employed to find numerous satisfactory solutions from a wide area of the design parameter space. The second step is to identify a precise structure of DW by clustering the detected satisfactory solutions with the modified K-means clustering algorithm (KMA). In the original KMA [16], each cluster center iteratively moves to form optimum clusters, fixing the number of clusters a priori. Here the objective function consists of the sum of distances between each cluster center and a location of each member belonging to the relevant cluster. As in many realistic engineering problems, the
number of clusters cannot be known a priori. Thus we newly develop the modified KMA with including an explicit mechanism to determine the optimum number of clusters during a clustering process. The proposed DW search method was implemented to an integrated computational aided engineering (CAE) system for multi-disciplinary structural design developed by some of the authors, and the method was applied to shape design of a microelectrostatic actuator for next generation high-density optical memory. A DW consisting of four clusters, i.e. four sub-DWs was obtained, and the features of the design solutions of the four subDWs were compared with each other in detail.
2. Design window search method 2.1. Definition of satisfactory design problem A satisfactory design problem is generally defined in the following. f ðXÞ is an objective function, gi ðXÞ is the ith constraint, and X (xi 2 Ri ) is a search vector. Ri is one-dimensional real space, and the domain to be searched is defined as R ¼ R1 R2 Rn . If f ðX 0 Þ is equal to or greater than a specified value, i.e. satisfactory value f0 , X 0 is considered as one of satisfactory solutions. The satisfactory value f0 is empirically defined by an engineer. As schematically illustrated in Fig. 1, the existing area of satisfactory solutions, i.e. DW may consist of multiple subgroups in multi-modal design problems. Since the illustration in this figure is twodimensional, it is easy to understand the situation. However it is very difficult to understand a precise structure of multiple solutions in multi-modal and multidimensional design problems. Thus we propose a new DW search method in this paper.
Fig. 1. Schematic illustration of DW consisting multiple subDWs.
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
2471
Fig. 2. Schematic flow of proposed two-steps DW search method.
2.2. Basic strategy We propose the following two-steps search method. Fig. 2 schematically shows its flow. First step: Detect satisfactory design solutions from all over a design parameter space by using the CEA. During this process, finite element (FE) analysis codes are directly used to evaluate a fitness function, i.e. an objective function, at each search point. This process is schematically illustrated in Fig. 2(a). It should be noted here that although a DW is basically continuous subspace within the multi-dimensional design parameter space, the DW is searched in a discrete manner, i.e. point by point by the CEA. Second step: In multi-modal problems, the DW may consist of multiple sub-DWs as explained in Section 2.1 and Fig. 1. Thus we cluster all the satisfactory design solutions detected in the first step into one to several subgroups by using a modified KMA. The one to severalsubgroups of the satisfactory solutions are regarded as sub-DWs. This process is also illustrated in Fig. 2(b).
2.3. Direct search of satisfactory solutions using CEA 2.3.1. Flow of CEA The detail of the CEA is described in Ref. [10]. Only a pooling procedure to pick up satisfactory solutions is additionally employed in the present study. The flow of the present CEA is illustrated in Fig. 3. At first, we define terms used in the algorithm. An individual means each search vector X ðpÞ (p ¼ 1; . . . ; N ), where N denotes the number of the individuals involved in a population. The population is defined as a school of individuals in each generation. A fitness of each individual X ðpÞ is defined as f ðX ðpÞ Þ, which corresponds to an objective function of a design problem. The main processes of the present CEA are briefly explained below: Initialization: The population is first initialized at random using the following expression:
Fig. 3. Flow of CEA with pooling procedure. ðpÞ
xi ð1 þ lÞx0i ;
ð1Þ
ðpÞ xi
is the ith component, i.e. ith design parameter, where of the pth individual, and x0i is the ith design parameter of a prototype. The prototype is an initial set of design parameters proposed empirically by an engineer, l is a random variable of a normal distribution C0 with a mean of 0. Evaluation and pooling: For each individual, the fitness f is calculated. If as f ðX ðpÞ Þ is equal to or greater than a satisfactory value f0 , the corresponding individual as X ðpÞ is pooled. Elite transfer: In the tth generation, k individuals with best fitness function values, which are called elite, are chosen from the population Pt , and are transferred to the next generation. Sexual reproduction and mutation: (N k) individuals are created by means of a sexual reproduction procedure. In each reproduction process, two different individuals X ðf Þ and X ðmÞ are chosen out of the population Pt at random, in which the probability of selection is proportional to the fitness f . Then they are recombined to create a single offspring X by the following expression: xi lc xfi þ ð1 lc Þxmi ;
ð2Þ
where the subscript i denotes the ith component of X, the superscripts f and m denote the components of the
2472
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
parental individuals X ðf Þ and X ðmÞ , and lc is a random valuable of a normal distribution Cc with a mean of 0. The offspring is next mutated with small probability and large magnitude by the following expression: xi
ð1 þ lm Þxi ;
ð3Þ
where xi is a randomly selected component, and lm is a random variable of a normal distribution Cm with a mean of 0. The above processes are repeated as many times as necessary to create (N k) new individuals. The newly created offspring individuals are merged with the previously hold elite individuals in order to form the (t þ 1)th population Ptþ1 . The evolution step is repeated until a criterion of convergence is satisfied. The criterion employed in this study is described in Section 2.3.3. 2.3.2. Detection of satisfactory solutions The CEA is a kind of multi-point search optimization algorithms for a continuous parameter space. When solving an ordinary optimization problem, a user is usually interested in only the one individual which has the best fitness function, i.e. the maximized objective function. In the process, a number of volumetric points are evaluated all over the design parameter space. Then numerous individuals whose fitness function values exceed the specified satisfactory value f0 can be found as well. These individuals are regarded as satisfactory design solutions. Since the multi-point search is almost at random performed over the whole design parameter space, it can be regarded that a finite number of the satisfactory solutions are the discrete expression of a continuous satisfactory solution subspace, i.e. DW. 2.3.3. Criterion of convergence The optimization process of the CEA is theoretically expected to converge to a unique optimum solution. However it often takes a long time to reach such a converged situation. In this paper, the following criterion of convergence is employed from a practical point of view. If the objective function is not improved during three successive generations, the population is re-initialized, preserving only the elite individuals. When the objective function is not improved even if this procedure is applied two times, we judge that the process is converged. The validity of this criterion of convergence has been confirmed in preliminary studies. 2.4. Clustering of satisfactory solutions 2.4.1. Modified K-means algorithm Let xi be the ith input pattern of a data set, where i ¼ 1; . . . ; n. n is the total number of data patterns. Each xi , consists of d design variables.
Step 1: Set K ¼ 2, and select the initial centers of K clusters. The mean vector, i.e. the center of kth cluster C ðkÞ is denoted as mðkÞ , where k ¼ 1; . . . ; K. Step 2: Compute the normalized Euclidean distance mðkÞ Di between xi and mðkÞ as follows: !1=2 d X mðkÞ ðkÞ 2 ðxij mj Þ ; ð4Þ Di ¼ j¼1 ðkÞ
where xij and mj are the jth components of xi and mðkÞ , respectively. Step 3: Assign xi , to the kth cluster C ðkÞ as one of elemðkÞ ments, if Di is the minimum among all the mðhÞ Di , h ¼ 1; . . . ; K. Step 4: Compute a within-cluster variation for each cluster in order to quantify similarity among all the elements. The within-cluster variation eðkÞ is defined as the sum of the Euclidean distances as follows: eðkÞ ¼
n X
mðkÞ
Di
:
ð5Þ
i¼1
Step 5: Compute the sum of the within-cluster variations for all the clusters as follows: Eq ¼
K X
eðkÞ ;
ð6Þ
k¼1
where Eq is the values evaluated at qth iteration. Step 6: Compute new cluster centers, and repeat steps 2 through 5 until the Eq value is minimized. The minimized value of Eq is called the variation level. The steps mentioned above are almost the same as a typical k-means algorithm, and it is necessary to define the number of clusters a priori. A proper cluster number is determined through the following steps. Step 7: Changing initial centers for d times, and repeat steps 2 through 6. Step 8: If each Eq obtained for different initial centers reaches the same minimum value, i.e. the same variation level, then we consider that K is the optimum number of clusters. Otherwise, we increase the number of clusters, i.e. K K þ 1, and go back to step 2. 2.4.2. Numerical examples To investigate fundamental performances of the clustering algorithm presented in the previous subsection, it is applied to solve two examples. Both examples have two design parameters and well separated data patterns for the purpose of visual verification of the method.
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Data Pattern Optimum Cluster Center
X2
Example 1. In the first example, the number of data patterns is n ¼ 75, and the number of trials for different initial centers is set to be d ¼ 4. Fig. 4 shows the profiles of the variation levels with respect to the iteration number. The figure clearly shows that the variation levels of different cluster numbers are well minimized, respectively. Fig. 5 shows the variation levels for different cluster numbers. The variation levels of K ¼ 2 and those of K ¼ 3 are the same, respectively, for different trials of initial centers. Then they start to differ at K ¼ 4. Since the variation level of K ¼ 3 is smaller than that of K ¼ 2, K ¼ 3 is judged to be the optimum number of clusters. Fig. 6 shows all the data patterns (open marks) with the three optimum centers (solid marks).
2473
X1
Example 2. The number of data patterns is n ¼ 35, and the number of trials of initial centers is d ¼ 3 in this example. Data patterns are also well separated, but one cluster has only a few patterns. Fig. 7 shows the variation levels for three trials of initial cluster centers with
Variation Level
10
8
Initial Cluster Center - 1 Initial Cluster Center - 2 Initial Cluster Center - 3
8
7
Variation Level
2 Clusters 3 Clusters 4 Clusters 5 Clusters
12
Fig. 6. Clustering result of Example 1.
6
5
4
3
6 2
3
4
5
Clusters Number 4 1
2
3
4
5
6
7
Fig. 7. Variation levels vs. cluster numbers of Example 2.
Iteration Number
Fig. 4. Variation levels vs. iteration numbers.
8
Initial Custers - 1 Initial Custers - 2 Initial Custers - 3 Initial Custers - 4
Variation Level
7
respect to different cluster numbers. In the figure, K ¼ 2 is judged to be the optimum because the variation levels for different cluster centers coincide among others only when K ¼ 2. Fig. 8 shows all the data patterns with the two optimum cluster centers.
3. Implementation of the proposed method into an integrated CAE system
6
5
4 2
3
4
Cluster Number
Fig. 5. Variation levels vs. cluster numbers.
5
In this study, we employ an integrated CAE system for multi-disciplinary structural design developed by some of the authors [17–19], where the CEA is utilized as an optimizer. Since the detailed description of each mechanism implemented in the system can be found in [17–19], only a brief description of the system is given here. This CAE system consists of (1) encapsulation mechanism of existing engineering codes, (2) welldefined common data representation, and (3) distributed
2474
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Fig. 9. Microelectrostatic actuator for optical memory.
Fig. 8. Clustering result of Example 2.
calculation framework based on mobile agent technology. The codes to be encapsulated in the system include FE analysis programs, optimization engines, mesh generators, and others. Since the system is designed and constructed based on an object-oriented approach, various codes can easily be integrated into it in a flexible manner. The system can also be executed in a distributed computer environment such as PC clusters. Although the CEA is a robust search algorithm, it might require a huge amount of evaluations of fitness function, i.e. an objective function, which inevitably includes time consuming FE calculations. Thus parallel and distributed evaluations are strongly required. The detail of the parallel issue of the present CAE system is described in [17–19]. Only its basic idea is described here. Running in a local CPU, the main body of the CEA defines individuals and sends each of them to remote CPUs. Each individual is then evaluated by using FE codes in each remote CPU, and only requested results are returned to the local CPU. The system configuration employed in the present paper and parallel efficiency are described in Section 5.1.
4. Application: structural design of microelectrostatic actuator for next generation optical memory The DW search method implemented to the integrated CAE system is applied to a structural design of a microelectrostatic actuator. The actuator is planned to be employed as memory pick-up mechanism for a super high-density optical memory system of next generation. The memory system is now under development [20]. 4.1. Operation principle Fig. 9 shows the birds-eye-view and top-view of a conceptual design of the actuator. It consists of lens,
comb type electrostatic actuator, supporting bars and flame. These components except lens are made of single crystal silicon through a photo-chemical etching process. The main body consisting of lens and actuator is hung to the flame by the four supporting bars, and is driven in the y-direction by electrostatic force acting between the fixed comb teeth and the movable ones. The actuator is to be designed to control one-directional nanometerorder displacement of the lens. A whole size of the actuator is about a few millimeters, while the smallest dimensions of components such as width of comb teeth and that of supporting bars are only a few micrometers. 4.2. Design requirements The following six kinds of design requirements (DR) are described for the actuator: DR1: DR2: DR3: DR4: DR5: DR6:
range of displacement, accuracy in positioning, strength, accuracy in chasing target memory, feedback control ability, manufacturability.
Table 1 summarizes the requirements DR1–DR5, where each requirement is expressed in terms of design characteristic parameter yi (i ¼ 1; . . . ; 5) with constraint. On the other hand, DR6 basically comes from difficulty in manufacturing micromachines. Considering restriction originated from a photo-chemical etching process, the shape should be as simple as possible. However it should be noted here that it is not so easy to formulate such a design requirement as DR6 in a mathematical form. 4.3. Design process A design process using the proposed DW search method is as follows. DR1–DR5 are first formulated as a satisfactory design problem, which is described in Section 4.4 . Second the proposed DW search method is
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
2475
Table 1 Design characteristic parameters and their restrictions y1 y2 y3 y4 y5
¼ ymax =Ymax > 1 ¼ ymin =Ymin < 1 ¼ sSmax =Sy < 1 ¼ a=A > 1 ¼ traise =Traise < 1
ymax : maximum displacement of lens center in y-direction, Ymax : required maximum displacement, ymin : minimum displacement of lens center in y-direction, Ymin : required minimum displacement, smax : maximum MisesÕs mean stress in structure, sy : yield stress of structure, s: safety factor, a: acceleration of lens center in y-direction, A: required acceleration.
applied to this problem and a DW consisting of multiple sub-DWs is detected. This process is described in Section 5.2. Finally a final design solution is selected out of the multiple sub-DWs, considering DR6. The final process is described in Section 5.3. Fig. 10. Typical mesh.
4.4. Formulation of satisfactory design problem This subsection describes a process to formulate DR1; . . . ; DR5 into a satisfactory design problem. 4.4.1. Evaluation of design requirements To evaluate DR1–DR4, three-dimensional static stress analyses are performed, where a unit load is applied to the movable comb teeth. To evaluate DR5, a dynamic response analysis considering air effects is performed. The equation of motion considered here is as follows: ð½M þ ½MF Þfag þ ½CF fvg þ ½K fug ¼ ff g;
ð7Þ
where [M] denotes a mass matrix of the actuator, [MF ] an added mass matrix caused due to the surrounding air, [CF ] an added viscosity matrix due to the surrounding air, and [K] a stiffness matrix of the actuator, fag denotes an acceleration vector, fvg a velocity vector, and fug a displacement vector of the actuator. [MF ] and [CF ] can be derived based on either theoretical analyses or numerical experiments. In this paper, we employed simplified theoretical estimation proposed in [21–23]. Eq. (7) is implemented as a FE code and integrated into the CAE system. 4.4.2. Analysis model A ten-noded tetrahedral mesh is employed, and the number of nodes is controlled to be less than a specified number. Such limitation of mesh size is necessary to achieve realistic design time. In this study, the upper limit number of nodes is defined as 5000, so that the evaluation time of FE analyses for one evaluation of fitness function is kept to be less than a few minutes. A typical mesh is shown in Fig. 10. Material properties of single crystal of silicon and those of air (under 1 N/m2 , 20 °C) are given in Table 2. Fig. 11 shows all the pa-
Table 2 Material properties Silicon[1 1 0] Mass density (kg/m3 ) YoungÕs modulus (GPa) PoissonÕs ratio Yield stress (GPa)
2328 170 0.2 2.0
Air (20 °C) Mass density (kg/m3 ) Kinematic viscosity (kg/m2 )
1.205 15:01 10 3
rameters defining the actuator shape. As a result of a preliminary sensitivity study, five most influential parameters b1 , b2 , k2 , k3 , and m1 are selected as design parameters to be optimized. 4.4.3. Objective function, satisfactory value and constraints Five kinds of satisfaction functions of fi , i ¼ 1; . . . ; 5 are empirically defined, corresponding to the five design requirements DR1–DR5 described in Table 1. Some of them are illustrated in Fig. 12. Here the score of 10 means that the requirement is fully satisfied, while the half score of 5 means that the requirement is satisfied and that the solution is acceptable. The principal aim of the present design problem is to find solutions that satisfy the conditions such that fi > 5, i ¼ 1; . . . ; 5. In order to apply the proposed DW search algorithm, this original problem to search the feasible solutions is converted to the following satisfaction problem: Objective function f ðb1 ; b2 ; k2 ; k3 ; m1 Þ ¼ minff1 ; . . . ; f5 g; Satisfactory value f0 ¼ 5: ð8Þ
2476
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Fig. 11. Shape parameters of actuator and their values of basic design.
MB RAM serving as a central server and the place for the execution of the DW search method.
Fig. 12. Examples of satisfaction function corresponding to design requirements DR1 and DR2.
Constraints on the design variables are defined so as to prevent invalid shapes.
5. Results and discussion 5.1. System performance In this subsection, to examine performance of the CAE system, the problem defined in Section 4.4.3 is solved several times by changing a set of experiments, population size and the number of hosts employed. 5.1.1. System configuration The computer system employed in this study consists of 25 PCs (Pentium 133 MHz, 64 MB RAM) used for remote evaluations and a PC Pentium Pro 200 MHz, 64
5.1.2. Calculation efficiency Parallel efficiency of the system is examined first. Fig. 13 shows the parallel efficiency in evaluating 16 design solutions on 1–16 hosts. Five design parameters were initialized at random using formula (1). Changing the number of hosts and the combination of design parameter deviation, 50 experiments were performed. The average parallel efficiency was calculated using the following formula: c ¼ t0 =ht
ð9Þ
in which h is the number of employed hosts, t is average evaluation time and t0 is average evaluation time on one host with the same design parameter deviation. As can be seen in Fig. 13, the parallel performance decreases to the level of 0.6 when each host evaluates only one design solution. When mesh size varies largely, performance decreases significantly. This is because the number of designs to be evaluated by each host is fixed in this calculation. To resolve such a situation, it may be effective to employ a dynamic workload balancing feature. It is also important to estimate total calculation time. An averaged total calculation time Tavg ðnÞ as a function of the number of hosts n is evaluated for a set of 19 experiments using various population sizes from 50 to 200 and various numbers of hosts from 9 to 25. The five
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
2477
Parallel efficiency
1.2 1 0.8 0.6 0.4 0.2 0 1
10 Number of hosts
100
Fig. 15. Transition of objective function.
Fig. 13. Parallel efficiency vs. number of hosts.
Averaged total calculation time
12 11 10 9 8 7 6 5 4 5
10
15
20
25
30
Number of hosts
Fig. 14. Averaged total calculation time vs. number of hosts.
design parameters were initialized at random using formula (1). Tavg ðhÞ is defined as follows: ðeÞ P Tavg ðhÞ ¼
T ðhÞ Ns
E e¼1
E
;
ð10Þ
where T is a total calculation time, N the number of generations, s the number of individuals in each generation, e denotes each experiment and E the total number of experiments. The result is shown in Fig. 14. As shown in the figure, Tavg ðnÞ changes almost linearly in accordance with the number of hosts.
Fig. 15 shows an evolution of the objective function using the CEA. The evolution is stopped at 40th generation, following the criterion described in Section 2.3.3. Total calculation time took about two days, in which calculation time per individual took about 130– 135 s. More than 80% of the searching time was consumed to perform FE analyses for the individuals. The objective function f exceeded the score of 5 at the 25th generation, and satisfactory solutions began to appear. Finally 179 satisfactory design solutions were obtained out of the total search points of 16,000 within the whole design parameter space. The highest score of f was 6.35. Next the modified K-means algorithm was applied to investigate the precise structure of the detected 179 satisfactory solutions. Fig. 16 shows calculated variation levels for different cluster numbers from K ¼ 2 to 5. The figure clearly shows that the variation levels for three different trials coincide with each other well in the case of K ¼ 4. Then K ¼ 4 was judged as optimum. The obtained four clusters, i.e. four sub-DWs, are named as cluster 1 (63), cluster 2 (46), cluster 3 (42), and cluster 4 (28), where each number in bracket denotes the number of solutions classified into each cluster. The solution
5.2. DW Search The basic design illustrated in Fig. 9, which was derived empirically, is employed as a reference design to generate the initial population. As described in Section 5.1.2, total calculation time decreases almost linearly in accordance with the number of hosts employed. 25 hosts are employed at most. The five design parameters described in Section 4.4.2 are initialized at random using formula (1). As described in Section 5.1.1, parallel efficiency increases in accordance with the number of individuals per host. The number of individuals in a population is then defined to be 200.
Fig. 16. Variation levels vs. cluster numbers.
2478
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
with the highest score of f ¼ 6:35 is included in cluster 3. In this study, a mean vector of ith cluster, i ¼ 1; . . . ; 4 is chosen as the typical satisfactory solution and denoted as Type i (i ¼ 1; . . . ; 4). Type 1; . . . ; Type 4 are called as candidates of a final design solution. Shapes of Type 1; . . . ; Type 4 and the basic design are illustrated in Fig. 17, while their corresponding values of dimensions, de-
sign characteristics and satisfaction functions are tabulated in Table 3. In a usual optimization process, only the unique solution with the highest value of the objective function f is adopted as the final design solution, while other satisfactory design solutions are thrown out. In this study, all of the detected 179 satisfactory solutions are con-
Fig. 17. Comparison of shapes among basic design and four candidate designs.
5.777 6.039 5.880 5.352
6.351
10.00 10.00 10.00 10.00
10.00 10.00 7.325 0.8649 21.16 0.4683 0.2464 0.9461 Best design using conventional optimization method 0.6824 0.0179 0.3565 0.0593
1.144 1.163 1.170 1.055 0.8113 0.9313 0.9147 0.6845 0.3328 0.3089 0.3162 0.3264 Type Type Type Type
Means of clusters Cluster center
1 2 3 4
0.7067 0.8660 0.6401 0.7681
0.0180 0.0230 0.0169 0.0189
0.0556 0.0527 0.0507 0.0557
0.30 0.10 0.285 0.015 0.20
1.232
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 6.436 6.629 6.689 5.548 0.9223 0.8961 0.9120 0.9648 16.86 18.34 17.69 13.81 0.4296 0.3727 0.4931 0.3315 0.2287 0.2325 0.2337 0.2109
0.00303 0.01517
y2 y1 m1 k3 k2 b2 b1 One of initial designs
10.00
10.00 0.0000 0.1350 )273.7 0.04891
y5 y4
f1
0.07584
f2
10.00
f3
10.00
f5
5.3. Detected solutions
f4
2479
sidered as members of a DW, and utilized effectively to determine the final solution. This selection process will be explained in the next subsection.
y3
Satisfaction function values Design characters Design parameters (mm)
Table 3 Comparisons of design parameters, design characters, satisfaction functions among an initial design and the obtained design solutions
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
The supporting bars and comb teeth of all the candidates are much longer in the x-direction and thinner in the z-direction than those of the basic design, respectively. As the results, all the candidates have less stiffness, i.e. more flexibility, than the basic design, and then attain sufficient displacement and smaller effects of air. Variations in their shapes result in those in their design characteristics as given in Table 3. It is confirmed from the results that the proposed method can search a DW even if the initial design is set far from the DW. Precise characteristics of each candidate are summarized as follows. Type 1 has longer and thinner supporting bars and comb teeth. Type 2 has shorter and thicker supporting bars and teeth. In Type 1 and Type 2, respectively, supporting bars and teeth of comb are almost the same in length. Type 3 has longer teeth than supporting bars, so that its fabrication process might be most complicated. On the contrary, Type 4 has shorter teeth than supporting bars, so that the shape of the flame can be a simple square. Considering the difficulty of fabrication processes of micromachines, it is desirable that the flame of the actuator has simple shape. Considering DR6, i.e. manufacturability, it is judged that Type 4 is the best, Type 2 is better, Type 1 is worse, and Type 3 is the worst. Satisfactory functions among the candidates are next compared. Here relative satisfactory functions f10 –f50 are defined. fi0 corresponds to the ith design characteristic yi , and the score of 10 is given to the best value of fi0 , while the score of 1 to the worst value. The radar charts of fi0 are shown in Fig. 18. It is well understood from the figure that Type1 is the most balanced and better solution, and that Type 2 has no weak points. On the other hand, Type 3 and Type 4 are superior in some specific characteristics. Through the above discussions, it is finally concluded that we should select Type 2 as a final design.
6. Conclusions The direct DW search method using the CEA and the modified K-means algorithm was proposed. To demonstrate its practical performance, the proposed method was implemented into the integrated CAE system for multi-disciplinary structural design running in a parallel computer environment such as PC clusters, and applied to the shape design of a microelectrostatic actuator. The actuator is planned to be employed as memory pick-up
2480
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Fig. 18. Comparison of relative satisfactory values (nos. 1–5 denote corresponding satisfaction values).
mechanism for a super high-density optical memory system of next generation. More than 80% of the searching time was consumed to perform FE analyses for the individuals. The basic characteristics of the employed parallel calculation were examined as follows. The performance decreases to the level of 0.6 when each host evaluate only one design solution. On the other hand, the total calculation time decreases almost linearly according to the increase of the number of hosts. In this paper, we employed 25 hosts to decrease the total calculation time, while employed 200 individuals in each generation to increase parallel efficiency. The developed system was applied to the structural design problem, which consisted of five design parameters and six design requirements. The five out of the six design requirements were formulated as a satisfactory design problem. 179 satisfactory design solutions were found out of the whole design parameter space, and then the solutions were classified into four clusters, i.e. four sub-DWs. Shapes and characteristics were examined in detail among representative design solutions taken from the four clusters. Finally we selected the best design solution, considering the sixth unformulated design requirement, i.e. manufacturability. Following the present study, the research group of Fujitsu Research Institute Tokyo has produced the prototype of a microelectrostatic actuator for optical memory.
Acknowledgements A part of this work was financially supported by the Promotion Society of Optical Industrial Technology.
The authors wish to thank Profs. K. Itao, H. Hosaka and T. Ohkubo at Univ. of Tokyo, and Mr. H. Komoriya at Fujitsu Research Institute Tokyo for valuable discussions.
References [1] Hasan MZ, Monde M, Mitsutake Y, Islam MA. Design window for shield thermal energy storage for pulsed tokamak fusion reactors. Fusion Eng Des 1998;41:547– 54. [2] Wang XM, Ye L, Mai YW, Galea SC. Designing for piezoelectric ceramic wafers bonded on structures using force transfer criteria. Smart Mater Struct 2000;9:157–62. [3] Mochizuki Y, Yoshimura S, Yagawa G. Automated system for structural design using design windows search approach: its application to fusion first wall design. Adv Eng Software 1997;28:103–13. [4] Lee JS, Yoshimura S, Yagawa G, Shibaike N. A CAE system for micromachines: its application to electrostatic micro wobble actuator. Sens Actuators A 1995;50:209–21. [5] Matsuda A, Yoshimuira S, Yagawa G, Hirata T, Nishioka M. KANSEI-based quantitative satisfactory design of artifacts using numeral networks (its application to subjective evaluation for handling and stability of vehicle). Trans Jpn Soc Mech Eng 1998;64C:1004–12 [in Japanese]. [6] Shi Q, Hagiwara I, Takashima F. Global optimization method to multiple local optimals with the response surface approximation methodology. Trans Jpn Soc Mech Eng 1999;65A:232–9 [in Japanese]. [7] Hagiwara I, Shi Q, Takashima F. Function approximation method for crash optimization using neural network (2nd report, application to vehicle component). Trans Jpn Soc Mech Eng 1998;64A:2441–7 [in Japanese]. [8] Kaminaga A, Suzuki K, Fujii D, Ohtsubo H. Response surface model using moving least squares method. OPTIS2000 2000;00-27:181–6 [in Japanese].
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481 [9] Nikolaidis E, Long L, Ling Q. Neural networks and response surface polynomials for design of vehicle joints. Comput Struct 2000;75(6):593–607. [10] Furukawa T, Yagawa G. Inelastic constitutive parameter identification using an evolutionary algorithm with continuous individuals. Int J Numer Meth Eng 1997;40(6): 1071–90. [11] Holand JH. Adaptation in natural and artificial systems. Michigan: The University of Michigan Press; 1975. [12] Fogal LJ, Oweens AJ, Walsh MJ. Artificial intelligence through simulated evolution. Willey; 1966. [13] Rechenberg I. Evolution strategie: Optimierung Tchnischer Systeme Nach Prinzipien der Bologischen Evolution. Frommann-Holzboog; 1973. [14] Schwefel HP. Numerical optimization of computer models. Willey; 1981. [15] Hoffmeister F, Back T. Genetic algorithms and evolution strategies and differences. Technical report. University of Dortmund, Germany 1992; Sys-1/92. [16] Markin B. Mathematical classification and clustering. Kluwer Academic Publishers; 1996.
2481
[17] Yoshimura S, Furukawa T, Kowalczyk T, Yagawa G. An automated CAE system for multidisciplinary structural design and its application to micro-accelerometer design. Trans Jpn Soc Mech Eng 1997;62C:2173–80 [in Japanese]. [18] Yoshimura S, Kowalczyk T. CAE in parallel distributed processing environment using JAVA. J Jpn Soc Comput Eng Sci 1998;3-2:94–101 [in Japanese]. [19] Kowalczyk T, Yoshimura S, Yagawa G. Integrated CAE system for micromachine design. JSCES Conf Comput Eng Sci 1999;2:547–50. [20] http://www.aist.go.jp/sangi/28.html. [21] Ishihara D, Yoshimura S, Yagawa G, Komoriya H. Automated design of micromachine for optical memory (application of integrated CAE system considering effects of fluid–structure interaction). JSME Des Syst ConfÕ99 1999;99-27:70–3 [in Japanese]. [22] Tomomasa T. Fluid dynamics. Baifukan; 1982. [in Japanese]. [23] Nagakura H, Kaneko N. The stability of a cantilever beam subjected to one-dimensional leakage flow. Trans Jpn Soc Mech Eng 1992;58C:352–9 [in Japanese].
Design window search using continuous evolutionary algorithm and clustering––its application to shape design of microelectrostatic actuator Daisuke Ishihara a, Min Joong Jeong b, Shinobu Yoshimura
b,*
, Genki Yagawa
a
a
b
School of Engineering, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Department of Quantum Engineering and Systems Science, Institute of Environmental Studies, University of Tokyo, 7-3-1 Hongo, Bunkyo-ku, Tokyo 113-8656, Japan Received 31 May 2001; accepted 25 July 2002
Abstract This paper describes a search method for multi-dimensional design window (DW), which is defined as an existing area of satisfactory solutions in a design parameter space. The method consists of the following two steps: (1) direct search for satisfactory design solutions using a continuous evolutionary algorithm (CEA) from a wide area of the design parameter space in a robust manner, and (2) identification of a precise structure of DW by clustering the detected satisfactory solutions with a modified K-means algorithm. The CEA is a kind of genetic algorithms modified to deal with continuous variables. The modified K-means clustering algorithm contains an explicit procedure to determine the optimum number of clusters. The proposed DW search method was implemented to an integrated computational aided engineering system for multi-disciplinary structural design developed by the authors, and then the method was applied to shape design of a microelectrostatic actuator for next generation high-density optical memory. A DW consisting of four clusters, i.e. four sub-DWs was obtained, and the features of the representative design solutions of the four subDWs were compared with each other in detail, and a final design solution was determined. Ó 2002 Elsevier Science Ltd. All rights reserved. Keywords: Structural design; Design window; Continuous evolutionary algorithm; K-means clustering; Finite element analysis; Microelectrostatic actuator
1. Introduction In a field of structural design, much research efforts have been put onto the development of an effective procedure to find an optimum solution, assuming that an objective function and all constraints are explicitly formulated in mathematical forms. In practical engineering situations, it is however natural that some of design requirements are hardly formulated in a mathe-
*
Corresponding author. Tel.: +81-3-3812-6960; fax: +81-35800-6876. E-mail address: [email protected] (S. Yoshimura).
matical form. Those unformulated requirements are taken into account in the last stage of design. Manufacturability and manufacturing cost may be some of such requirements. In these cases, it seems more meaningful to obtain satisfactory solutions that meet all the formulated design requirements rather than obtaining a single optimum solution. Thus we focus on a design window (DW) search method in this paper. The DW schematically indicates an existing area of satisfactory solutions within a design parameter space. Some DW search methods have been proposed and applied in several kinds of engineering problems [1–5]. In Ref. [1], a DW projected to a two-dimensional space is detected using a few points search method. In Ref. [2], a DW is detected using a simplified mathematical
0045-7949/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 4 5 - 7 9 4 9 ( 0 2 ) 0 0 2 9 3 - 6
2470
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
model approximating an original design problem. These approaches are too simple to apply to such general design problems that a relation between design parameters and design characteristics is nonlinear, discontinuous, multi-dimensional and multi-modal. Some of the present authors also studied a DW search method [3–5]. In Ref. [3], a whole area search method with a help of a multi-layer neural network was proposed, and applied to some engineering applications [4,5]. Here an orthogonal lattice is first generated in a design parameter space, and each lattice point is evaluated whether to be a satisfactory solution or not. The multi-layer neural network was used to approximate an actual relation between design parameters and design characteristics for the purpose of shortening DW searching time. Being similar to response surface methods, such an approximate method might introduce unallowable errors when solving highly nonlinear and discontinuous design problems [6–8]. Ref. [9] reported that errors were larger than 50% in some design problems. In this study, we newly propose a direct DW search method that is applicable to general design problems with nonlinear, discontinuous, multi-dimensional and/or multi-modal features. This is a two-steps method using a continuous evolutionary algorithm (CEA) and a modified K-means algorithm. At the first step, using the CEA, satisfactory solutions are searched out of a whole area of the design parameter space in a robust matter. The CEA is a kind of genetic algorithms (GAs) modified to deal with continuous variables and is good at solving illposed design problems [10]. After a sufficient number of evaluations, it is well expected that all the obtained satisfactory solutions may be a discrete expression of a continuous DW. Evolutionary algorithms (EAs), which have been named to represent GA [11], evolutionary programming [12], evolution strategy [13,14] and their combined algorithm [15], have appeared as robust optimization techniques in a last few decades. There are many studies on engineering applications of the EAs. The CEA, which was developed by Furukawa and Yagawa [10] is one of the GAs for real search space. In the past researches, the GAs have been employed to find a unique optimum solution. In the present study, the CEA is, on the other hand, employed to find numerous satisfactory solutions from a wide area of the design parameter space. The second step is to identify a precise structure of DW by clustering the detected satisfactory solutions with the modified K-means clustering algorithm (KMA). In the original KMA [16], each cluster center iteratively moves to form optimum clusters, fixing the number of clusters a priori. Here the objective function consists of the sum of distances between each cluster center and a location of each member belonging to the relevant cluster. As in many realistic engineering problems, the
number of clusters cannot be known a priori. Thus we newly develop the modified KMA with including an explicit mechanism to determine the optimum number of clusters during a clustering process. The proposed DW search method was implemented to an integrated computational aided engineering (CAE) system for multi-disciplinary structural design developed by some of the authors, and the method was applied to shape design of a microelectrostatic actuator for next generation high-density optical memory. A DW consisting of four clusters, i.e. four sub-DWs was obtained, and the features of the design solutions of the four subDWs were compared with each other in detail.
2. Design window search method 2.1. Definition of satisfactory design problem A satisfactory design problem is generally defined in the following. f ðXÞ is an objective function, gi ðXÞ is the ith constraint, and X (xi 2 Ri ) is a search vector. Ri is one-dimensional real space, and the domain to be searched is defined as R ¼ R1 R2 Rn . If f ðX 0 Þ is equal to or greater than a specified value, i.e. satisfactory value f0 , X 0 is considered as one of satisfactory solutions. The satisfactory value f0 is empirically defined by an engineer. As schematically illustrated in Fig. 1, the existing area of satisfactory solutions, i.e. DW may consist of multiple subgroups in multi-modal design problems. Since the illustration in this figure is twodimensional, it is easy to understand the situation. However it is very difficult to understand a precise structure of multiple solutions in multi-modal and multidimensional design problems. Thus we propose a new DW search method in this paper.
Fig. 1. Schematic illustration of DW consisting multiple subDWs.
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
2471
Fig. 2. Schematic flow of proposed two-steps DW search method.
2.2. Basic strategy We propose the following two-steps search method. Fig. 2 schematically shows its flow. First step: Detect satisfactory design solutions from all over a design parameter space by using the CEA. During this process, finite element (FE) analysis codes are directly used to evaluate a fitness function, i.e. an objective function, at each search point. This process is schematically illustrated in Fig. 2(a). It should be noted here that although a DW is basically continuous subspace within the multi-dimensional design parameter space, the DW is searched in a discrete manner, i.e. point by point by the CEA. Second step: In multi-modal problems, the DW may consist of multiple sub-DWs as explained in Section 2.1 and Fig. 1. Thus we cluster all the satisfactory design solutions detected in the first step into one to several subgroups by using a modified KMA. The one to severalsubgroups of the satisfactory solutions are regarded as sub-DWs. This process is also illustrated in Fig. 2(b).
2.3. Direct search of satisfactory solutions using CEA 2.3.1. Flow of CEA The detail of the CEA is described in Ref. [10]. Only a pooling procedure to pick up satisfactory solutions is additionally employed in the present study. The flow of the present CEA is illustrated in Fig. 3. At first, we define terms used in the algorithm. An individual means each search vector X ðpÞ (p ¼ 1; . . . ; N ), where N denotes the number of the individuals involved in a population. The population is defined as a school of individuals in each generation. A fitness of each individual X ðpÞ is defined as f ðX ðpÞ Þ, which corresponds to an objective function of a design problem. The main processes of the present CEA are briefly explained below: Initialization: The population is first initialized at random using the following expression:
Fig. 3. Flow of CEA with pooling procedure. ðpÞ
xi ð1 þ lÞx0i ;
ð1Þ
ðpÞ xi
is the ith component, i.e. ith design parameter, where of the pth individual, and x0i is the ith design parameter of a prototype. The prototype is an initial set of design parameters proposed empirically by an engineer, l is a random variable of a normal distribution C0 with a mean of 0. Evaluation and pooling: For each individual, the fitness f is calculated. If as f ðX ðpÞ Þ is equal to or greater than a satisfactory value f0 , the corresponding individual as X ðpÞ is pooled. Elite transfer: In the tth generation, k individuals with best fitness function values, which are called elite, are chosen from the population Pt , and are transferred to the next generation. Sexual reproduction and mutation: (N k) individuals are created by means of a sexual reproduction procedure. In each reproduction process, two different individuals X ðf Þ and X ðmÞ are chosen out of the population Pt at random, in which the probability of selection is proportional to the fitness f . Then they are recombined to create a single offspring X by the following expression: xi lc xfi þ ð1 lc Þxmi ;
ð2Þ
where the subscript i denotes the ith component of X, the superscripts f and m denote the components of the
2472
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
parental individuals X ðf Þ and X ðmÞ , and lc is a random valuable of a normal distribution Cc with a mean of 0. The offspring is next mutated with small probability and large magnitude by the following expression: xi
ð1 þ lm Þxi ;
ð3Þ
where xi is a randomly selected component, and lm is a random variable of a normal distribution Cm with a mean of 0. The above processes are repeated as many times as necessary to create (N k) new individuals. The newly created offspring individuals are merged with the previously hold elite individuals in order to form the (t þ 1)th population Ptþ1 . The evolution step is repeated until a criterion of convergence is satisfied. The criterion employed in this study is described in Section 2.3.3. 2.3.2. Detection of satisfactory solutions The CEA is a kind of multi-point search optimization algorithms for a continuous parameter space. When solving an ordinary optimization problem, a user is usually interested in only the one individual which has the best fitness function, i.e. the maximized objective function. In the process, a number of volumetric points are evaluated all over the design parameter space. Then numerous individuals whose fitness function values exceed the specified satisfactory value f0 can be found as well. These individuals are regarded as satisfactory design solutions. Since the multi-point search is almost at random performed over the whole design parameter space, it can be regarded that a finite number of the satisfactory solutions are the discrete expression of a continuous satisfactory solution subspace, i.e. DW. 2.3.3. Criterion of convergence The optimization process of the CEA is theoretically expected to converge to a unique optimum solution. However it often takes a long time to reach such a converged situation. In this paper, the following criterion of convergence is employed from a practical point of view. If the objective function is not improved during three successive generations, the population is re-initialized, preserving only the elite individuals. When the objective function is not improved even if this procedure is applied two times, we judge that the process is converged. The validity of this criterion of convergence has been confirmed in preliminary studies. 2.4. Clustering of satisfactory solutions 2.4.1. Modified K-means algorithm Let xi be the ith input pattern of a data set, where i ¼ 1; . . . ; n. n is the total number of data patterns. Each xi , consists of d design variables.
Step 1: Set K ¼ 2, and select the initial centers of K clusters. The mean vector, i.e. the center of kth cluster C ðkÞ is denoted as mðkÞ , where k ¼ 1; . . . ; K. Step 2: Compute the normalized Euclidean distance mðkÞ Di between xi and mðkÞ as follows: !1=2 d X mðkÞ ðkÞ 2 ðxij mj Þ ; ð4Þ Di ¼ j¼1 ðkÞ
where xij and mj are the jth components of xi and mðkÞ , respectively. Step 3: Assign xi , to the kth cluster C ðkÞ as one of elemðkÞ ments, if Di is the minimum among all the mðhÞ Di , h ¼ 1; . . . ; K. Step 4: Compute a within-cluster variation for each cluster in order to quantify similarity among all the elements. The within-cluster variation eðkÞ is defined as the sum of the Euclidean distances as follows: eðkÞ ¼
n X
mðkÞ
Di
:
ð5Þ
i¼1
Step 5: Compute the sum of the within-cluster variations for all the clusters as follows: Eq ¼
K X
eðkÞ ;
ð6Þ
k¼1
where Eq is the values evaluated at qth iteration. Step 6: Compute new cluster centers, and repeat steps 2 through 5 until the Eq value is minimized. The minimized value of Eq is called the variation level. The steps mentioned above are almost the same as a typical k-means algorithm, and it is necessary to define the number of clusters a priori. A proper cluster number is determined through the following steps. Step 7: Changing initial centers for d times, and repeat steps 2 through 6. Step 8: If each Eq obtained for different initial centers reaches the same minimum value, i.e. the same variation level, then we consider that K is the optimum number of clusters. Otherwise, we increase the number of clusters, i.e. K K þ 1, and go back to step 2. 2.4.2. Numerical examples To investigate fundamental performances of the clustering algorithm presented in the previous subsection, it is applied to solve two examples. Both examples have two design parameters and well separated data patterns for the purpose of visual verification of the method.
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Data Pattern Optimum Cluster Center
X2
Example 1. In the first example, the number of data patterns is n ¼ 75, and the number of trials for different initial centers is set to be d ¼ 4. Fig. 4 shows the profiles of the variation levels with respect to the iteration number. The figure clearly shows that the variation levels of different cluster numbers are well minimized, respectively. Fig. 5 shows the variation levels for different cluster numbers. The variation levels of K ¼ 2 and those of K ¼ 3 are the same, respectively, for different trials of initial centers. Then they start to differ at K ¼ 4. Since the variation level of K ¼ 3 is smaller than that of K ¼ 2, K ¼ 3 is judged to be the optimum number of clusters. Fig. 6 shows all the data patterns (open marks) with the three optimum centers (solid marks).
2473
X1
Example 2. The number of data patterns is n ¼ 35, and the number of trials of initial centers is d ¼ 3 in this example. Data patterns are also well separated, but one cluster has only a few patterns. Fig. 7 shows the variation levels for three trials of initial cluster centers with
Variation Level
10
8
Initial Cluster Center - 1 Initial Cluster Center - 2 Initial Cluster Center - 3
8
7
Variation Level
2 Clusters 3 Clusters 4 Clusters 5 Clusters
12
Fig. 6. Clustering result of Example 1.
6
5
4
3
6 2
3
4
5
Clusters Number 4 1
2
3
4
5
6
7
Fig. 7. Variation levels vs. cluster numbers of Example 2.
Iteration Number
Fig. 4. Variation levels vs. iteration numbers.
8
Initial Custers - 1 Initial Custers - 2 Initial Custers - 3 Initial Custers - 4
Variation Level
7
respect to different cluster numbers. In the figure, K ¼ 2 is judged to be the optimum because the variation levels for different cluster centers coincide among others only when K ¼ 2. Fig. 8 shows all the data patterns with the two optimum cluster centers.
3. Implementation of the proposed method into an integrated CAE system
6
5
4 2
3
4
Cluster Number
Fig. 5. Variation levels vs. cluster numbers.
5
In this study, we employ an integrated CAE system for multi-disciplinary structural design developed by some of the authors [17–19], where the CEA is utilized as an optimizer. Since the detailed description of each mechanism implemented in the system can be found in [17–19], only a brief description of the system is given here. This CAE system consists of (1) encapsulation mechanism of existing engineering codes, (2) welldefined common data representation, and (3) distributed
2474
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Fig. 9. Microelectrostatic actuator for optical memory.
Fig. 8. Clustering result of Example 2.
calculation framework based on mobile agent technology. The codes to be encapsulated in the system include FE analysis programs, optimization engines, mesh generators, and others. Since the system is designed and constructed based on an object-oriented approach, various codes can easily be integrated into it in a flexible manner. The system can also be executed in a distributed computer environment such as PC clusters. Although the CEA is a robust search algorithm, it might require a huge amount of evaluations of fitness function, i.e. an objective function, which inevitably includes time consuming FE calculations. Thus parallel and distributed evaluations are strongly required. The detail of the parallel issue of the present CAE system is described in [17–19]. Only its basic idea is described here. Running in a local CPU, the main body of the CEA defines individuals and sends each of them to remote CPUs. Each individual is then evaluated by using FE codes in each remote CPU, and only requested results are returned to the local CPU. The system configuration employed in the present paper and parallel efficiency are described in Section 5.1.
4. Application: structural design of microelectrostatic actuator for next generation optical memory The DW search method implemented to the integrated CAE system is applied to a structural design of a microelectrostatic actuator. The actuator is planned to be employed as memory pick-up mechanism for a super high-density optical memory system of next generation. The memory system is now under development [20]. 4.1. Operation principle Fig. 9 shows the birds-eye-view and top-view of a conceptual design of the actuator. It consists of lens,
comb type electrostatic actuator, supporting bars and flame. These components except lens are made of single crystal silicon through a photo-chemical etching process. The main body consisting of lens and actuator is hung to the flame by the four supporting bars, and is driven in the y-direction by electrostatic force acting between the fixed comb teeth and the movable ones. The actuator is to be designed to control one-directional nanometerorder displacement of the lens. A whole size of the actuator is about a few millimeters, while the smallest dimensions of components such as width of comb teeth and that of supporting bars are only a few micrometers. 4.2. Design requirements The following six kinds of design requirements (DR) are described for the actuator: DR1: DR2: DR3: DR4: DR5: DR6:
range of displacement, accuracy in positioning, strength, accuracy in chasing target memory, feedback control ability, manufacturability.
Table 1 summarizes the requirements DR1–DR5, where each requirement is expressed in terms of design characteristic parameter yi (i ¼ 1; . . . ; 5) with constraint. On the other hand, DR6 basically comes from difficulty in manufacturing micromachines. Considering restriction originated from a photo-chemical etching process, the shape should be as simple as possible. However it should be noted here that it is not so easy to formulate such a design requirement as DR6 in a mathematical form. 4.3. Design process A design process using the proposed DW search method is as follows. DR1–DR5 are first formulated as a satisfactory design problem, which is described in Section 4.4 . Second the proposed DW search method is
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
2475
Table 1 Design characteristic parameters and their restrictions y1 y2 y3 y4 y5
¼ ymax =Ymax > 1 ¼ ymin =Ymin < 1 ¼ sSmax =Sy < 1 ¼ a=A > 1 ¼ traise =Traise < 1
ymax : maximum displacement of lens center in y-direction, Ymax : required maximum displacement, ymin : minimum displacement of lens center in y-direction, Ymin : required minimum displacement, smax : maximum MisesÕs mean stress in structure, sy : yield stress of structure, s: safety factor, a: acceleration of lens center in y-direction, A: required acceleration.
applied to this problem and a DW consisting of multiple sub-DWs is detected. This process is described in Section 5.2. Finally a final design solution is selected out of the multiple sub-DWs, considering DR6. The final process is described in Section 5.3. Fig. 10. Typical mesh.
4.4. Formulation of satisfactory design problem This subsection describes a process to formulate DR1; . . . ; DR5 into a satisfactory design problem. 4.4.1. Evaluation of design requirements To evaluate DR1–DR4, three-dimensional static stress analyses are performed, where a unit load is applied to the movable comb teeth. To evaluate DR5, a dynamic response analysis considering air effects is performed. The equation of motion considered here is as follows: ð½M þ ½MF Þfag þ ½CF fvg þ ½K fug ¼ ff g;
ð7Þ
where [M] denotes a mass matrix of the actuator, [MF ] an added mass matrix caused due to the surrounding air, [CF ] an added viscosity matrix due to the surrounding air, and [K] a stiffness matrix of the actuator, fag denotes an acceleration vector, fvg a velocity vector, and fug a displacement vector of the actuator. [MF ] and [CF ] can be derived based on either theoretical analyses or numerical experiments. In this paper, we employed simplified theoretical estimation proposed in [21–23]. Eq. (7) is implemented as a FE code and integrated into the CAE system. 4.4.2. Analysis model A ten-noded tetrahedral mesh is employed, and the number of nodes is controlled to be less than a specified number. Such limitation of mesh size is necessary to achieve realistic design time. In this study, the upper limit number of nodes is defined as 5000, so that the evaluation time of FE analyses for one evaluation of fitness function is kept to be less than a few minutes. A typical mesh is shown in Fig. 10. Material properties of single crystal of silicon and those of air (under 1 N/m2 , 20 °C) are given in Table 2. Fig. 11 shows all the pa-
Table 2 Material properties Silicon[1 1 0] Mass density (kg/m3 ) YoungÕs modulus (GPa) PoissonÕs ratio Yield stress (GPa)
2328 170 0.2 2.0
Air (20 °C) Mass density (kg/m3 ) Kinematic viscosity (kg/m2 )
1.205 15:01 10 3
rameters defining the actuator shape. As a result of a preliminary sensitivity study, five most influential parameters b1 , b2 , k2 , k3 , and m1 are selected as design parameters to be optimized. 4.4.3. Objective function, satisfactory value and constraints Five kinds of satisfaction functions of fi , i ¼ 1; . . . ; 5 are empirically defined, corresponding to the five design requirements DR1–DR5 described in Table 1. Some of them are illustrated in Fig. 12. Here the score of 10 means that the requirement is fully satisfied, while the half score of 5 means that the requirement is satisfied and that the solution is acceptable. The principal aim of the present design problem is to find solutions that satisfy the conditions such that fi > 5, i ¼ 1; . . . ; 5. In order to apply the proposed DW search algorithm, this original problem to search the feasible solutions is converted to the following satisfaction problem: Objective function f ðb1 ; b2 ; k2 ; k3 ; m1 Þ ¼ minff1 ; . . . ; f5 g; Satisfactory value f0 ¼ 5: ð8Þ
2476
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Fig. 11. Shape parameters of actuator and their values of basic design.
MB RAM serving as a central server and the place for the execution of the DW search method.
Fig. 12. Examples of satisfaction function corresponding to design requirements DR1 and DR2.
Constraints on the design variables are defined so as to prevent invalid shapes.
5. Results and discussion 5.1. System performance In this subsection, to examine performance of the CAE system, the problem defined in Section 4.4.3 is solved several times by changing a set of experiments, population size and the number of hosts employed. 5.1.1. System configuration The computer system employed in this study consists of 25 PCs (Pentium 133 MHz, 64 MB RAM) used for remote evaluations and a PC Pentium Pro 200 MHz, 64
5.1.2. Calculation efficiency Parallel efficiency of the system is examined first. Fig. 13 shows the parallel efficiency in evaluating 16 design solutions on 1–16 hosts. Five design parameters were initialized at random using formula (1). Changing the number of hosts and the combination of design parameter deviation, 50 experiments were performed. The average parallel efficiency was calculated using the following formula: c ¼ t0 =ht
ð9Þ
in which h is the number of employed hosts, t is average evaluation time and t0 is average evaluation time on one host with the same design parameter deviation. As can be seen in Fig. 13, the parallel performance decreases to the level of 0.6 when each host evaluates only one design solution. When mesh size varies largely, performance decreases significantly. This is because the number of designs to be evaluated by each host is fixed in this calculation. To resolve such a situation, it may be effective to employ a dynamic workload balancing feature. It is also important to estimate total calculation time. An averaged total calculation time Tavg ðnÞ as a function of the number of hosts n is evaluated for a set of 19 experiments using various population sizes from 50 to 200 and various numbers of hosts from 9 to 25. The five
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
2477
Parallel efficiency
1.2 1 0.8 0.6 0.4 0.2 0 1
10 Number of hosts
100
Fig. 15. Transition of objective function.
Fig. 13. Parallel efficiency vs. number of hosts.
Averaged total calculation time
12 11 10 9 8 7 6 5 4 5
10
15
20
25
30
Number of hosts
Fig. 14. Averaged total calculation time vs. number of hosts.
design parameters were initialized at random using formula (1). Tavg ðhÞ is defined as follows: ðeÞ P Tavg ðhÞ ¼
T ðhÞ Ns
E e¼1
E
;
ð10Þ
where T is a total calculation time, N the number of generations, s the number of individuals in each generation, e denotes each experiment and E the total number of experiments. The result is shown in Fig. 14. As shown in the figure, Tavg ðnÞ changes almost linearly in accordance with the number of hosts.
Fig. 15 shows an evolution of the objective function using the CEA. The evolution is stopped at 40th generation, following the criterion described in Section 2.3.3. Total calculation time took about two days, in which calculation time per individual took about 130– 135 s. More than 80% of the searching time was consumed to perform FE analyses for the individuals. The objective function f exceeded the score of 5 at the 25th generation, and satisfactory solutions began to appear. Finally 179 satisfactory design solutions were obtained out of the total search points of 16,000 within the whole design parameter space. The highest score of f was 6.35. Next the modified K-means algorithm was applied to investigate the precise structure of the detected 179 satisfactory solutions. Fig. 16 shows calculated variation levels for different cluster numbers from K ¼ 2 to 5. The figure clearly shows that the variation levels for three different trials coincide with each other well in the case of K ¼ 4. Then K ¼ 4 was judged as optimum. The obtained four clusters, i.e. four sub-DWs, are named as cluster 1 (63), cluster 2 (46), cluster 3 (42), and cluster 4 (28), where each number in bracket denotes the number of solutions classified into each cluster. The solution
5.2. DW Search The basic design illustrated in Fig. 9, which was derived empirically, is employed as a reference design to generate the initial population. As described in Section 5.1.2, total calculation time decreases almost linearly in accordance with the number of hosts employed. 25 hosts are employed at most. The five design parameters described in Section 4.4.2 are initialized at random using formula (1). As described in Section 5.1.1, parallel efficiency increases in accordance with the number of individuals per host. The number of individuals in a population is then defined to be 200.
Fig. 16. Variation levels vs. cluster numbers.
2478
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
with the highest score of f ¼ 6:35 is included in cluster 3. In this study, a mean vector of ith cluster, i ¼ 1; . . . ; 4 is chosen as the typical satisfactory solution and denoted as Type i (i ¼ 1; . . . ; 4). Type 1; . . . ; Type 4 are called as candidates of a final design solution. Shapes of Type 1; . . . ; Type 4 and the basic design are illustrated in Fig. 17, while their corresponding values of dimensions, de-
sign characteristics and satisfaction functions are tabulated in Table 3. In a usual optimization process, only the unique solution with the highest value of the objective function f is adopted as the final design solution, while other satisfactory design solutions are thrown out. In this study, all of the detected 179 satisfactory solutions are con-
Fig. 17. Comparison of shapes among basic design and four candidate designs.
5.777 6.039 5.880 5.352
6.351
10.00 10.00 10.00 10.00
10.00 10.00 7.325 0.8649 21.16 0.4683 0.2464 0.9461 Best design using conventional optimization method 0.6824 0.0179 0.3565 0.0593
1.144 1.163 1.170 1.055 0.8113 0.9313 0.9147 0.6845 0.3328 0.3089 0.3162 0.3264 Type Type Type Type
Means of clusters Cluster center
1 2 3 4
0.7067 0.8660 0.6401 0.7681
0.0180 0.0230 0.0169 0.0189
0.0556 0.0527 0.0507 0.0557
0.30 0.10 0.285 0.015 0.20
1.232
10.00 10.00 10.00 10.00 10.00 10.00 10.00 10.00 6.436 6.629 6.689 5.548 0.9223 0.8961 0.9120 0.9648 16.86 18.34 17.69 13.81 0.4296 0.3727 0.4931 0.3315 0.2287 0.2325 0.2337 0.2109
0.00303 0.01517
y2 y1 m1 k3 k2 b2 b1 One of initial designs
10.00
10.00 0.0000 0.1350 )273.7 0.04891
y5 y4
f1
0.07584
f2
10.00
f3
10.00
f5
5.3. Detected solutions
f4
2479
sidered as members of a DW, and utilized effectively to determine the final solution. This selection process will be explained in the next subsection.
y3
Satisfaction function values Design characters Design parameters (mm)
Table 3 Comparisons of design parameters, design characters, satisfaction functions among an initial design and the obtained design solutions
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
The supporting bars and comb teeth of all the candidates are much longer in the x-direction and thinner in the z-direction than those of the basic design, respectively. As the results, all the candidates have less stiffness, i.e. more flexibility, than the basic design, and then attain sufficient displacement and smaller effects of air. Variations in their shapes result in those in their design characteristics as given in Table 3. It is confirmed from the results that the proposed method can search a DW even if the initial design is set far from the DW. Precise characteristics of each candidate are summarized as follows. Type 1 has longer and thinner supporting bars and comb teeth. Type 2 has shorter and thicker supporting bars and teeth. In Type 1 and Type 2, respectively, supporting bars and teeth of comb are almost the same in length. Type 3 has longer teeth than supporting bars, so that its fabrication process might be most complicated. On the contrary, Type 4 has shorter teeth than supporting bars, so that the shape of the flame can be a simple square. Considering the difficulty of fabrication processes of micromachines, it is desirable that the flame of the actuator has simple shape. Considering DR6, i.e. manufacturability, it is judged that Type 4 is the best, Type 2 is better, Type 1 is worse, and Type 3 is the worst. Satisfactory functions among the candidates are next compared. Here relative satisfactory functions f10 –f50 are defined. fi0 corresponds to the ith design characteristic yi , and the score of 10 is given to the best value of fi0 , while the score of 1 to the worst value. The radar charts of fi0 are shown in Fig. 18. It is well understood from the figure that Type1 is the most balanced and better solution, and that Type 2 has no weak points. On the other hand, Type 3 and Type 4 are superior in some specific characteristics. Through the above discussions, it is finally concluded that we should select Type 2 as a final design.
6. Conclusions The direct DW search method using the CEA and the modified K-means algorithm was proposed. To demonstrate its practical performance, the proposed method was implemented into the integrated CAE system for multi-disciplinary structural design running in a parallel computer environment such as PC clusters, and applied to the shape design of a microelectrostatic actuator. The actuator is planned to be employed as memory pick-up
2480
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481
Fig. 18. Comparison of relative satisfactory values (nos. 1–5 denote corresponding satisfaction values).
mechanism for a super high-density optical memory system of next generation. More than 80% of the searching time was consumed to perform FE analyses for the individuals. The basic characteristics of the employed parallel calculation were examined as follows. The performance decreases to the level of 0.6 when each host evaluate only one design solution. On the other hand, the total calculation time decreases almost linearly according to the increase of the number of hosts. In this paper, we employed 25 hosts to decrease the total calculation time, while employed 200 individuals in each generation to increase parallel efficiency. The developed system was applied to the structural design problem, which consisted of five design parameters and six design requirements. The five out of the six design requirements were formulated as a satisfactory design problem. 179 satisfactory design solutions were found out of the whole design parameter space, and then the solutions were classified into four clusters, i.e. four sub-DWs. Shapes and characteristics were examined in detail among representative design solutions taken from the four clusters. Finally we selected the best design solution, considering the sixth unformulated design requirement, i.e. manufacturability. Following the present study, the research group of Fujitsu Research Institute Tokyo has produced the prototype of a microelectrostatic actuator for optical memory.
Acknowledgements A part of this work was financially supported by the Promotion Society of Optical Industrial Technology.
The authors wish to thank Profs. K. Itao, H. Hosaka and T. Ohkubo at Univ. of Tokyo, and Mr. H. Komoriya at Fujitsu Research Institute Tokyo for valuable discussions.
References [1] Hasan MZ, Monde M, Mitsutake Y, Islam MA. Design window for shield thermal energy storage for pulsed tokamak fusion reactors. Fusion Eng Des 1998;41:547– 54. [2] Wang XM, Ye L, Mai YW, Galea SC. Designing for piezoelectric ceramic wafers bonded on structures using force transfer criteria. Smart Mater Struct 2000;9:157–62. [3] Mochizuki Y, Yoshimura S, Yagawa G. Automated system for structural design using design windows search approach: its application to fusion first wall design. Adv Eng Software 1997;28:103–13. [4] Lee JS, Yoshimura S, Yagawa G, Shibaike N. A CAE system for micromachines: its application to electrostatic micro wobble actuator. Sens Actuators A 1995;50:209–21. [5] Matsuda A, Yoshimuira S, Yagawa G, Hirata T, Nishioka M. KANSEI-based quantitative satisfactory design of artifacts using numeral networks (its application to subjective evaluation for handling and stability of vehicle). Trans Jpn Soc Mech Eng 1998;64C:1004–12 [in Japanese]. [6] Shi Q, Hagiwara I, Takashima F. Global optimization method to multiple local optimals with the response surface approximation methodology. Trans Jpn Soc Mech Eng 1999;65A:232–9 [in Japanese]. [7] Hagiwara I, Shi Q, Takashima F. Function approximation method for crash optimization using neural network (2nd report, application to vehicle component). Trans Jpn Soc Mech Eng 1998;64A:2441–7 [in Japanese]. [8] Kaminaga A, Suzuki K, Fujii D, Ohtsubo H. Response surface model using moving least squares method. OPTIS2000 2000;00-27:181–6 [in Japanese].
D. Ishihara et al. / Computers and Structures 80 (2002) 2469–2481 [9] Nikolaidis E, Long L, Ling Q. Neural networks and response surface polynomials for design of vehicle joints. Comput Struct 2000;75(6):593–607. [10] Furukawa T, Yagawa G. Inelastic constitutive parameter identification using an evolutionary algorithm with continuous individuals. Int J Numer Meth Eng 1997;40(6): 1071–90. [11] Holand JH. Adaptation in natural and artificial systems. Michigan: The University of Michigan Press; 1975. [12] Fogal LJ, Oweens AJ, Walsh MJ. Artificial intelligence through simulated evolution. Willey; 1966. [13] Rechenberg I. Evolution strategie: Optimierung Tchnischer Systeme Nach Prinzipien der Bologischen Evolution. Frommann-Holzboog; 1973. [14] Schwefel HP. Numerical optimization of computer models. Willey; 1981. [15] Hoffmeister F, Back T. Genetic algorithms and evolution strategies and differences. Technical report. University of Dortmund, Germany 1992; Sys-1/92. [16] Markin B. Mathematical classification and clustering. Kluwer Academic Publishers; 1996.
2481
[17] Yoshimura S, Furukawa T, Kowalczyk T, Yagawa G. An automated CAE system for multidisciplinary structural design and its application to micro-accelerometer design. Trans Jpn Soc Mech Eng 1997;62C:2173–80 [in Japanese]. [18] Yoshimura S, Kowalczyk T. CAE in parallel distributed processing environment using JAVA. J Jpn Soc Comput Eng Sci 1998;3-2:94–101 [in Japanese]. [19] Kowalczyk T, Yoshimura S, Yagawa G. Integrated CAE system for micromachine design. JSCES Conf Comput Eng Sci 1999;2:547–50. [20] http://www.aist.go.jp/sangi/28.html. [21] Ishihara D, Yoshimura S, Yagawa G, Komoriya H. Automated design of micromachine for optical memory (application of integrated CAE system considering effects of fluid–structure interaction). JSME Des Syst ConfÕ99 1999;99-27:70–3 [in Japanese]. [22] Tomomasa T. Fluid dynamics. Baifukan; 1982. [in Japanese]. [23] Nagakura H, Kaneko N. The stability of a cantilever beam subjected to one-dimensional leakage flow. Trans Jpn Soc Mech Eng 1992;58C:352–9 [in Japanese].