Multidisciplinary System Multidisciplinary ystem Design Optimization (MSDO) Design Optimization (MSDO) Simulated Annealing Lecture 5 March 2004 by Dr. Cyrus D. Jilla Professor Olivier de Weck ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 1 Today’s Topics Today’s Topics ? What are Heuristics? ? The Origin and Analogy of Simulated Annealing ? The Simulated Annealing Algorithm ? An Example: The Terrestrial Planet Finder Mission ? Sample Results ? Ways to Tailor the Simulated Algorithm ? Summary ? References ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 2 What is a Heuristic? What is a Heuristic? ? A Heuristic is simply a rule of thumb that hopefully will find a good answer. ? Why use a Heuristic? – Heuristics are typically used to solve complex (large, nonlinear, nonconvex (ie. contain many local minima)) multivariate combinatorial optimization problems that are difficult to solve to optimality. ? Unlike gradient-based methods (such as the simplex algorithm) in a convex trade space, heuristics are NOT guaranteed to find the true global optimal solution in a single objective problem, but should find many good solutions (the mathematician's answer vs. the engineer’s answer) ? Heuristics are good at dealing with local optima without getting stuck in them while searching for the global optimum. ? Reference: – Schulz, A.S., “Metaheuristics,” 15.057 Systems Optimization Course Notes, MIT, 1999. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 3 Types of Heuristics Types of Heuristics ? Heuristics Often Incorporate Randomization ? 2 Special Cases of Heuristics – Construction Methods ? Must first find a feasible solution and then improve it. – Improvement Methods ? Start with a feasible solution and just try to improve it. ? 3 Most Common Heuristic Techniques – Genetic Algorithms – Simulated Annealing – Tabu Search ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 4 Origin of Simulated Annealing (SA) Origin of Simulated Annealing (SA) ? Definition: A heuristic technique that mathematically mirrors the cooling of a material to a state of minimum energy. ? Origin: Applying the field of Statistical Mechanics to the field of Combinatorial Optimization (1983) ? Draws an analogy between the cooling of a metal and the solving of an optimization problem. ? Original Paper Introducing the Concept – Kirkpatrick, S., Gelatt, C.D., and Vecchi, M.P., “Optimization by Simulated Annealing,” Science, Volume 220, Number 4598, 13 May 1983, pp. 671-680. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 5 The Analogy The Analogy ? Statistical Mechanics: The behavior of systems with many degrees of freedom in thermal equilibrium at a finite temperature. ? Combinatorial Optimization: Finding the minimum of a given function depending on many parameters. ? Analogy: If a liquid material (ie. metal) cools and anneals too quickly, then the material will solidify into a sub-optimal configuration. If the liquid material cools slowly, the crystals within the material will solidify optimally into a state of minimum energy (ie. ground state). – This ground state corresponds to the minimum of the cost function in an optimization problem. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 6 4 Key Ingredients for Simulated Annealing 4 Key Ingredients for Simulated Annealing (Kirkpatrick et al, 1983) Kirkpatrick et al, 1983) ? A concise description of a configuration (ie. architecture) of the system (Design Vector). ? A random generator of rearrangements of the elements in a configuration (Neighborhoods). ? A quantitative objective function containing the trade-offs that have to be made (Simulation Model and Output Metric(s)). ? An annealing schedule of the temperatures and the length of times for which the system is to be evolved. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 7 The SA Algorithm The SA Algorithm ? Terminology: – Γ = Design Vector (ie. Design Architecture) – E = System Energy (ie. Objective Function Value) – T = System Temperature – ? = Difference in System Energy Between Two Design Vectors ? The Simulated Annealing Algorithm 1) Choose a random Γ , select the initial system temperature, and outline the cooling i (ie. annealing) schedule 2) Evaluate E(Γ ) using your simulation model i 3) Perturb Γ to obtain a neighboring Design Vector (Γ i+1 ) i 4) Evaluate E(Γ i+1 ) using your simulation model 5) If E(Γ i+1 )< E(Γ ), Γ i+1 is the new current solution i 6) If E(Γ i+1 )> E(Γ ), then accept Γ i+1 as the new current solution with a probability e (- ?/T) i where ? = E(Γ i+1 ) - E(Γ ). i 7) Reduce the system temperature according to the cooling schedule. 8) Terminate the algorithm. ? TPF Example: We will walk through each of the 8 steps. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 8 Terrestrial Planet Finder (TPF) Terrestrial Planet Finder (TPF) ? Mission Statement: “To study all aspects of planets ranging from their formation and development in disks of dust and gas around newly forming stars to the presence and features of those planets orbiting the nearest stars. Specifically, to conduct a search for Earth-like planets in star systems located within 15 parsecs of our solar system.” ? Primary Mission: To detect Earth-like planets around nearby stars, especially those in the habitable zone where liquid water is likely to exist – Bracewell Nulling interferometer ? Secondary Mission: To characterize approximately 50 of these Earth-like planets – Medium spectroscopy (50 planets) – Detailed spectroscopy (5 planets) ? To image astrophysical structures to within milli-arcsecond angular resolution (Michelson interferometer) requires longer baselines ? http://tpf.jpl.nasa.gov/ ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 9 Mathematical Formulation of the TPF Design Problem TPF Design Vector: G = [γ 1 γ 2 γ 3 γ 4 ] ? 640 unique TPF mission Variable Allowable Values architectures Heliocentric Orbital Altitude (AU) 1.0, 1.5, 2.0, 2.5, 3.0, γ 1 3.5, 4.0, 4.5, 5.0, 5.5 ? Complete enumeration Collector Connectivity/Geometry SCI-1D, SCI-2D, SSI γ 2 1D, SSI-2D determines accuracy # Collector Apertures 4, 6, 8, 10 γ 3 Diameter of Collector Apertures(m) 1, 2, 3, 4 γ 4 ? MDO methods limited to the Original Optimization Formulation evaluation of 48 design 5 ∑ F y () vectors (7.5% of the full- G :Objective Min y=1 factorial trade space) 5 ∑ ? y ()G y=1 Where :sConstraint toSubject y = year in the mission Isolation 2.5 ≤ θ r ≤ milli20 - arcsec Φ= cost O ≤ 10 ?6 Ψ= number of “images” Integrity Surveying SNR ≥ 5 θ r = angular resolution pySpectroscoMed. SNR ≥ 10 ?= null depth pySpectroscoDeep SNR ≥ 25 ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 10 SA Step 1 – Starting Out SA Step 1 – Starting Out ? 2 Methods to Select the Initial Design Vector Γ i – Select a random point (ie. random design architecture) in the trade space – Select a known point in the trade space (ie. a design architecture which has already been studied/demonstrated) ? Typically, a random Design Vector is used to begin the SA algorithm (Exception: Construction Problem with few known feasible solutions). ? After a “good” solution(s) has been found, this can be used as the starting point for a new run of the SA algorithm to try and find an even better solution (Assumption: You have now found the most favorable portion of the trade space). ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 11 SA Step 2 - Evaluate E(Γ ) Using Your ) Using Your i Simulation Model Design Vector Operations Systems GINA Spacecraft Bus/Payload Aperture Configuration Environment Dynamics, Optics, Control, & Structures Metrics Orbit 1 AU Interferometer Type SCI-1D Number of Apertures 4 Size of Apertures 4 m Total Cost $902M Total # of Images 1273 Cost/Image $709K ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 12 SA Step 3 - Perturb SA Step 3 - Perturb Γ To Obtain a Neighboring To Obtain a Neighboring i Design Vector (Γ i+1 ) ? What is a Neighbor? – Γ 2 is a neighbor of Γ 1 if Γ 2 can be obtained from a modification of Γ 1 . ? Neighborhood Degrees of Freedom (DOF) – The number of parameters that can vary between Γ 1 and Γ 2 for which they may still be considered neighbors. ? Neighbor Example: G 1 = [ SCIAU1 m4ap41D- ] G = [ SCIAU1 2m4ap1D- ] DOF = 1 2 ? The perturbation to a neighbor is usually done randomly. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 13 SA Step 4 - Evaluate E(Γ i+1 ) Using Your Simulation Model Design Vector Operations Systems GINA Spacecraft Bus/Payload Aperture Configuration Environment Dynamics, Optics, Control, & Structures Metrics Orbit 1 AU Interferometer Type SCI-1D Number of Apertures 4 Size of Apertures 2 m Total Cost $769M Total # of Images 769 Cost/Image $1000K ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 14 SA Steps 5 & 6 – Accepting a New SA Steps 5 & 6 – Accepting a New Current Solution Current Solution ? The system energy E(Γ) (ie. the Objective Function) is a quantitative measure of the “goodness” of a complex system. ? The current solution is the baseline design architecture from which we explore the trade space. It is the design vector we perturb to find a new neighbor. ? If E(Γ i+1 )< E(Γ ), Γ i+1 is the new current solution. i ? If E(Γ i+1 )> E(Γ ), then accept Γ i+1 as the new current solution with i a probability e (- ?/T) where ? = E(Γ i+1 ) - E(Γ ). i ? e (- ?/T) is called the Boltzman Factor ? This is known as the Metropolis Step – Metropolis,N., A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, "Equation of State Calculations by Fast Computing Machines", J. Chem. Phys.,21, 6, 1087-1092, 1953. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 15 SA Steps 5 & 6 – Accepting a New SA Steps 5 & 6 – Accepting a New Current Solution Cont’d Current Solution Cont’d ? Why move to a worse current solution? –To avoid getting trapped in a local optimum. ? Local Optimum – A solution is locally optimal if there is no neighbor who has a better objective function value. ? Global Optimum – A solution is globally optimal if there is no other solution in the entire feasible trade space that has a better objective function value. – Note: We are only talking about single objective problems. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 16 Current Solution Cont’d Current Solution Cont’d SA Steps 5 & 6 – Accepting a New ? Local Optima vs. the Global Optimum ? Example – f(x) = cos(14.5x – 0.3) + (x + 0.2)x – Min[f(x)] Global Optimum Local Optima ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 17 SA Steps 5 & 6 – Accepting a New SA Steps 5 & 6 – Accepting a New Current Solution Cont’d Current Solution Cont’d ? Local Optima vs. the Global Optimum Example ? Min [cos(14.5x – 0.3) + (x + 0.2)x] Subject to -2.2 ≥ x ≥ -3 Global Optimum Local Optimum Potential Neighbors Gradient Gradient Simulated Annealing ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 18 Current Solution Cont’d Current Solution Cont’d SA Steps 5 & 6 – Accepting a New ? Can there be a difference between the current solution and the best solution? Simulated Annealing Algorithm Exploration of the TPF Trade Space - DOF=3 800 1000 1200 1400 1600 1800 2000 2200 2400 L i fecycl e Cost ( $ m illions) Current Solution Path Best Solution Path $2M/Image $1M/Image $0.5M/Image $0.25M/Image 0 500 1000 1500 2000 2500 3000 3500 4000 Performance (total # images ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 19 SA Steps 5 & 6 – Accepting a New SA Steps 5 & 6 – Accepting a New Current Solution Cont’d Current Solution Cont’d ? Returning to the TPF Example – E(Γ 1 ) = $709K – E(Γ 2 ) = $1000K ? If E(Γ i+1 )< E(Γ i ), Γ i+1 is the new current solution. – Not the Case: ? $1000K > $709K ? E(Γ 2 ) > E(Γ 1 ) ? Therefore Γ 1 remains the current solution. ? But if E(Γ i+1 )> E(Γ i ), then accept Γ i+1 as the new current solution with a probability ? ? ? ? ? ? T ? Prob e ? ( ( where =? GE i + 1 )? GE ) i ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 20 SA Steps 5 & 6 – Accepting a New SA Steps 5 & 6 – Accepting a New Current Solution Cont’d Current Solution Cont’d ? ? = E(Γ i+1 ) - E(Γ ) i ? ? = $1000K - $709K = $291K ? How does the system temperature come into play? ? The probability of moving to a worse solution = e (- ?/T) –AsT →∞, Probability(moving to a worse solution ) → 1 –AsT → 0, Probability(moving to a worse solution ) → 0 ? T o = 1000 T reduce = 0.95 ? The probability with which we move to this worse solution with the hope of escaping a possible local minimum decreases over time, and is a function of how many iterations we have executed (ie. How cool the system is. Remember the thermodynamic analogy.). Iteration # T Prob e(-?/T) 1 950 0.74 10 599 0.62 20 358 0.44 30 215 0.26 40 129 0.10 50 77 0.02 ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 21 SA Steps 5 & 6 – Accepting a New SA Steps 5 & 6 – Accepting a New Current Solution Cont’d Current Solution Cont’d ? Another Key Point – The larger the difference ? in system energy between the two neighbors, the lower the probability of moving to that neighbor. T = 215 Iteration #30 $1600K ? 0.0006 Prob e(-?/T) T = 215 $800K 0.024 T = 215 $400K 0.16 T = 215 $200K 0.39 T = 215 $100K 0.63 T = 215 $50K 0.79 ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 22 SA Step 7 – Reducing System Temperature SA Step 7 – Reducing System Temperature ? The initial system temperature is typically large – on the order of magnitude of the expected range of the objective function. ? Two Most Common Cooling Schedules – Explore several perturbations at a given temperature, then reduce the temperature to a predetermined value, and repeat. – Reduce the temperature between each perturbation, but by a smaller amount. ? Creating a cooling schedule is more of an art than a science, and may involve trial and error. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 23 SA Step 8 - Terminating the SA Algorithm SA Step 8 - Terminating the SA Algorithm ? How do you terminate the simulated annealing algorithm? – Depends upon the cooling/annealing schedule. ? If spending multiple iterations at each temperature, terminate algorithm when T=0. ? If spending only one iteration at each temperature: – If a new “best solution” is not found after a given number of iterations, then terminate the algorithm. – Set upper bound on the total number of iterations allowed. Terminate the algorithm when this upper bound is reached. This should correspond ~ T=0. – Manually terminate the algorithm. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 24 Review of the SA Algorithm Review of the SA Algorithm ? Terminology: – Γ = Design Vector (ie. Design Architecture) – E = System Energy (ie. Objective Function Value) – T = System Temperature – ? = Difference in System Energy Between Two Design Vectors ? The Simulated Annealing Algorithm 1) Choose a random Γ , select the initial system temperature, and outline the i cooling (ie. annealing) schedule. 2) Evaluate E(Γ ) using your simulation model i 3) Perturb Γ to obtain a neighboring Design Vector (Γ i+1 ) i 4) Evaluate E(Γ i+1 ) using your simulation model 5) If E(Γ i+1 )< E(Γ ), Γ i+1 is the new current solution i 6) If E(Γ i+1 )> E(Γ ), then accept Γ i+1 as the new current solution with a i probability e (- ?/T) where ? = E(Γ i+1 ) - E(Γ ). i 7) Reduce the system temperature according to the cooling schedule. 8) Terminate the algorithm. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 25 SA Results on TPF Trade Space 2300 2400 2500 2600 2700 2800 2900 3000 1150 1200 1250 1300 1350 1400 Zoom in of the TPF System Trade Space Performance (total # of images) Lifecycle Cost ($millions) $0.55M/Image $0.5M/Image $0.45M/Image $0.4M/Image 0 500 1000 1500 2000 2500 3000 3500 4000 800 1000 1200 1400 1600 1800 2000 2200 2400 TPF System Trade Space Performance (total # of images) Lifecy cle Cost ($millions) $2M/Image $1M/Image $0.5M/Image $0.25M/Image ? SA Algorithm finds most cost effective families of design architectures after evaluating only a small fraction of the global trade space. ? Simulated Annealing has been successfully applied to very large combinatorial optimization problems where the total number of variables may range into the tens of thousands. – The Traveling Salesman Problem ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 26 Performance Performance Degrees of Freedom (DOF) and SA ? Objective: “What DOF allows simulated annealing to search the DSS trade space most efficiently and why?” ? Recall Neighborhood Definition: G 1 = [ SCIAU1 m4ap41D- ] G = [ SCIAU1 2m4ap1D- ] DOF = 1 2 DOF=1 DOF=3 DOF=5 Initial Initial Initial $ Solution Optimal Solution Optimal Solution $ Solution Optimal Solution $ Solution Performance Performance Performance ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 27 Degrees of Freedom (DOF) and SA Performance Performance 1000 900 Simulated Annealing Performance vs. Trial Size 1DOF 2DOF 3DOF 4DOF Performance of the Simulated Annealing Algorithm in Finding the C P I MSE ($K/Image) 2 800 Minimum CPI TPF Design Architecture 700 as a Function of DOF 600 500 400 300 0 200 400 600 800 1000 Trial Size ? Observations: – Setting DOF=2 appears to yield the most efficient exploration of the TPF trade space. – DOF=Min. → Overly constrictive neighborhood – DOF=Max. → Degenerates to a random search ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 28 Degrees of Freedom (DOF) and SA Degrees of Freedom (DOF) and SA Performance Performance Simulated Annealing Algorithm Exploration of the TPF Trade Space - DOF=1 Simulated Annealing Algorithm Exploration of the TPF Trade Space - DOF=2 800 Current Solution Path Best Solution Path $2M/Image $1M/Image $0.5M/Image $0.25M/Image 0 500 1000 1500 2000 2500 3000 3500 4000 Performance (total # images Simulated Annealing Algorithm Exploration of the TPF Trade Space - DOF=3 800 0 500 1000 1500 2000 2500 3000 3500 4000 Performance (total # images Simulated Annealing Algorithm Exploration of the TPF Trade Space - DOF=4 2200 Current Solution Path Best Solution Path $2M/Image $1M/Image $0.5M/Image $0.25M/Image 2200 2400 2400 Lifecycle Cost ($millions) 2000 Lifecycle Cost ($millions) 2000 1800 1800 1600 1600 1400 1400 1200 1200 1000 1000 2200 Current Solution Path Best Solution Path $2M/Image $1M/Image $0.5M/Image $0.25M/Image 2200 2400 2400 Lifecycle Cost ($millions) 2000 Lifecycle Cost ($millions) 2000 1800 1800 1600 1600 1400 1400 1200 1200 1000 1000 800 800 Current Solution Path Best Solution Path $2M/Image $1M/Image $0.5M/Image $0.25M/Image 0 500 1000 1500 2000 2500 3000 3500 4000 0 500 1000 1500 2000 2500 3000 3500 4000 Performance (total # images Performance (total # images ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 29 Simulated Annealing Lecture Summary Simulated Annealing Lecture Summary ? Originated from statistical mechanics. Analogous to the cooling of a material to a state of minimum energy. ? 8-Step Algorithm ? Terrestrial Planet Finder (TPF) Mission Example ? There are many items one may tailor within the algorithm to affect its performance (initial system temp, cooling schedule, DOF within a neighborhood, etc.) ? The simulated annealing algorithm, like other heuristic techniques, is NOT guaranteed to find the global optimum; but it should find several good solutions. ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 30 References References ? Brooks, R.R., Iyenger, S.S., and Rai, S, “Comparison of Genetic Algorithms and Simulated Annealing for Cost Minimization in a Multisensor System,” Optical Engineering, Vol. 37(2), Feb. 1998, pp. 505-516. ? Cerny, V., "Thermodynamical Approach to the Traveling Salesman Problem: An Efficient Simulation Algorithm", J. Opt. Theory Appl., 45, 1, 41-51, 1985 . ? Jilla, C.D., Miller, D.W., and Sedwick, R.J., “Application of Multidisciplinary Design Optimization Techniques to Distributed Satellite Systems,” Journal of Spacecraft and Rockets,” Vol. 37, No.4, 2000, pp. 481-490. ? Jilla, C.D., and Miller, D.W., “Assessing the Performance of a Heuristic Simulated Annealing Algorithm for the Design of Distributed Satellite Systems,” Acta Astronautica, Vol. 48, No. 5-12, 2001, pp. 529-543. ? Kirkpatrick, S., Gelatt, C.D., and Vecchi, M.P., “Optimization by Simulated Annealing,” Science, Volume 220, Number 4598, 13 May 1983, pp. 671-680. ? Metropolis,N., A. Rosenbluth, M. Rosenbluth, A. Teller, E. Teller, "Equation of State Calculations by Fast Computing Machines", J. Chem. Phys.,21, 6, 1087-1092, 1953. ? Schulz, A.S., “Metaheuristics,” 15.057 Systems Optimization Course Notes, MIT, 1999. ? http://www-2.cs.cmu.edu/afs/cs.cmu.edu/project/anneal/www/tech_reports.html ? http://www.taygeta.com/annealing/simanneal.html ? Massachusetts Institute of Technology – Dr. Cyrus D. Jilla & Prof. Olivier de Weck Engineering Systems Division and Dept. of Aeronautics & Astronautics 31