Simulated Annealing Tutorial

Main.SimulatedAnnealing History

Hide minor edits - Show changes to output

June 21, 2020, at 04:47 AM by 136.36.211.159 -
Deleted lines 338-356:

----

(:html:)
 <div id="disqus_thread"></div>
    <script type="text/javascript">
        /* * * CONFIGURATION VARIABLES: EDIT BEFORE PASTING INTO YOUR WEBPAGE * * */
        var disqus_shortname = 'apmonitor'; // required: replace example with your forum shortname

        /* * * DON'T EDIT BELOW THIS LINE * * */
        (function() {
            var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true;
            dsq.src = 'https://' + disqus_shortname + '.disqus.com/embed.js';
            (document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq);
        })();
    </script>
    <noscript>Please enable JavaScript to view the <a href="https://disqus.com/?ref_noscript">comments powered by Disqus.</a></noscript>
    <a href="https://disqus.com" class="dsq-brlink">comments powered by <span class="logo-disqus">Disqus</span></a>
(:htmlend:)
January 29, 2020, at 07:43 PM by 12.45.189.171 -
Added lines 31-334:

(:toggle hide sapython button show="Show Python Source":)
(:div id=sapython:)
(:source lang=python:)
## Generate a contour plot
# Import some other libraries that we'll need
# matplotlib and numpy packages must also be installed
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
import random
import math

# define objective function
def f(x):
    x1 = x[0]
    x2 = x[1]
    obj = 0.2 + x1**2 + x2**2 - 0.1*math.cos(6.0*3.1415*x1) - 0.1*math.cos(6.0*3.1415*x2)
    return obj

# Start location
x_start = [0.8, -0.5]

# Design variables at mesh points
i1 = np.arange(-1.0, 1.0, 0.01)
i2 = np.arange(-1.0, 1.0, 0.01)
x1m, x2m = np.meshgrid(i1, i2)
fm = np.zeros(x1m.shape)
for i in range(x1m.shape[0]):
    for j in range(x1m.shape[1]):
        fm[i][j] = 0.2 + x1m[i][j]**2 + x2m[i][j]**2 \
            - 0.1*math.cos(6.0*3.1415*x1m[i][j]) \
            - 0.1*math.cos(6.0*3.1415*x2m[i][j])

# Create a contour plot
plt.figure()
# Specify contour lines
#lines = range(2,52,2)
# Plot contours
CS = plt.contour(x1m, x2m, fm)#,lines)
# Label contours
plt.clabel(CS, inline=1, fontsize=10)
# Add some text to the plot
plt.title('Non-Convex Function')
plt.xlabel('x1')
plt.ylabel('x2')

##################################################
# Simulated Annealing
##################################################
# Number of cycles
n = 50
# Number of trials per cycle
m = 50
# Number of accepted solutions
na = 0.0
# Probability of accepting worse solution at the start
p1 = 0.7
# Probability of accepting worse solution at the end
p50 = 0.001
# Initial temperature
t1 = -1.0/math.log(p1)
# Final temperature
t50 = -1.0/math.log(p50)
# Fractional reduction every cycle
frac = (t50/t1)**(1.0/(n-1.0))
# Initialize x
x = np.zeros((n+1,2))
x[0] = x_start
xi = np.zeros(2)
xi = x_start
na = na + 1.0
# Current best results so far
xc = np.zeros(2)
xc = x[0]
fc = f(xi)
fs = np.zeros(n+1)
fs[0] = fc
# Current temperature
t = t1
# DeltaE Average
DeltaE_avg = 0.0
for i in range(n):
    print('Cycle: ' + str(i) + ' with Temperature: ' + str(t))
    for j in range(m):
        # Generate new trial points
        xi[0] = xc[0] + random.random() - 0.5
        xi[1] = xc[1] + random.random() - 0.5
        # Clip to upper and lower bounds
        xi[0] = max(min(xi[0],1.0),-1.0)
        xi[1] = max(min(xi[1],1.0),-1.0)
        DeltaE = abs(f(xi)-fc)
        if (f(xi)>fc):
            # Initialize DeltaE_avg if a worse solution was found
            #  on the first iteration
            if (i==0 and j==0): DeltaE_avg = DeltaE
            # objective function is worse
            # generate probability of acceptance
            p = math.exp(-DeltaE/(DeltaE_avg * t))
            # determine whether to accept worse point
            if (random.random()<p):
                # accept the worse solution
                accept = True
            else:
                # don't accept the worse solution
                accept = False
        else:
            # objective function is lower, automatically accept
            accept = True
        if (accept==True):
            # update currently accepted solution
            xc[0] = xi[0]
            xc[1] = xi[1]
            fc = f(xc)
            # increment number of accepted solutions
            na = na + 1.0
            # update DeltaE_avg
            DeltaE_avg = (DeltaE_avg * (na-1.0) +  DeltaE) / na
    # Record the best x values at the end of every cycle
    x[i+1][0] = xc[0]
    x[i+1][1] = xc[1]
    fs[i+1] = fc
    # Lower the temperature for next cycle
    t = frac * t

# print solution
print('Best solution: ' + str(xc))
print('Best objective: ' + str(fc))

plt.plot(x[:,0],x[:,1],'y-o')
plt.savefig('contour.png')

fig = plt.figure()
ax1 = fig.add_subplot(211)
ax1.plot(fs,'r.-')
ax1.legend(['Objective'])
ax2 = fig.add_subplot(212)
ax2.plot(x[:,0],'b.-')
ax2.plot(x[:,1],'g--')
ax2.legend(['x1','x2'])

# Save the figure as a PNG
plt.savefig('iterations.png')

plt.show()
(:sourceend:)
(:divend:)

(:toggle hide samatlab button show="Show MATLAB Source":)
(:div id=samatlab:)
(:source lang=matlab:)
clc;
clear;
close all;

%% Generate a contour plot
% Start location
x_start = [0.8, -0.5];

% Design variables at mesh points
i1 = -1.0:0.01:1.0;
i2 = -1.0:0.01:1.0;
[x1m, x2m] = meshgrid(i1, i2);
fm = 0.2 + x1m.^2 + x2m.^2 - 0.1*cos(6.0*pi*x1m) - 0.1*cos(6.0*pi*x2m);

% Contour Plot
fig = figure(1);
[C,h] = contour(x1m,x2m,fm);
clabel(C,h,'Labelspacing',250);
title('Simulated Annealing');
xlabel('x1');
ylabel('x2');
hold on;

%% Simulated Annealing

% Number of cycles
n = 50;
% Number of trials per cycle
m = 50;
% Number of accepted solutions
na = 0.0;
% Probability of accepting worse solution at the start
p1 = 0.7;
% Probability of accepting worse solution at the end
p50 = 0.001;
% Initial temperature
t1 = -1.0/log(p1);
% Final temperature
t50 = -1.0/log(p50);
% Fractional reduction every cycle
frac = (t50/t1)^(1.0/(n-1.0));
% Initialize x
x = zeros(n+1,2);
x(1,:) = x_start;
xi = x_start;
na = na + 1.0;
% Current best results so far
% xc = x(1,:);                         
fc = f(xi);
fs = zeros(n+1,1);
fs(1,:) = fc;
% Current temperature
t = t1;
% DeltaE Average
DeltaE_avg = 0.0;
for i=1:n
    disp(['Cycle: ',num2str(i),' with Temperature: ',num2str(t)])
    xc(1) = x(i,1);                       
    xc(2) = x(i,2);                       
    for j=1:m
        % Generate new trial points
        xi(1) = xc(1) + rand() - 0.5;
        xi(2) = xc(2) + rand() - 0.5;
        % Clip to upper and lower bounds
        xi(1) = max(min(xi(1),1.0),-1.0);
        xi(2) = max(min(xi(2),1.0),-1.0);
        DeltaE = abs(f(xi)-fc);
        if (f(xi)>fc)
            % Initialize DeltaE_avg if a worse solution was found
            %  on the first iteration
            if (i==1 && j==1)
                DeltaE_avg = DeltaE;
            end
            % objective function is worse
            % generate probability of acceptance
            p = exp(-DeltaE/(DeltaE_avg * t));
            %            % determine whether to accept worse point
            if (rand()<p)
                % accept the worse solution
                accept = true;
            else
                % don't accept the worse solution
                accept = false;
            end
        else
            % objective function is lower, automatically accept
            accept = true;
        end

        % accept                           
        if (accept==true)
            % update currently accepted solution
            xc(1) = xi(1);
            xc(2) = xi(2);
            fc = f(xc);
                                           
            xa(j,1) = xc(1);               
            xa(j,2) = xc(2);               
            fa(j) = f(xc);                 
            % increment number of accepted solutions
            na = na + 1.0;
            % update DeltaE_avg
            DeltaE_avg = (DeltaE_avg * (na-1.0) +  DeltaE) / na;
        else
            fa(j) = 0.0;                   
        end
    end
    % cycle)
    fa_Min_Index = find(nonzeros(fa) == min(nonzeros(fa)));
                                           
    if isempty(fa_Min_Index) == 0         
        x(i+1,1) = xa(fa_Min_Index,1); 
        x(i+1,2) = xa(fa_Min_Index,2); 
        fs(i+1)  = fa(fa_Min_Index);   
    else                                   
        x(i+1,1) = x(i,1);                 
        x(i+1,2) = x(i,2);                 
        fs(i+1)  = fs(i);                 
    end                                   
    % Lower the temperature for next cycle
    t = frac * t;
    fa = 0.0;                             
end
% print solution
disp(['Best candidate: ',num2str(xc)])
disp(['Best solution: ',num2str(fc)])
plot(x(:,1),x(:,2),'r-o')
saveas(fig,'contour','png')

fig = figure(2);
subplot(2,1,1)
plot(fs,'r.-')
legend('Objective')
subplot(2,1,2)
hold on
plot(x(:,1),'b.-')
plot(x(:,2),'g.-')
legend('x_1','x_2')

% Save the figure as a PNG
saveas(fig,'iterations','png')
(:sourceend:)

'''Save as f.m'''

(:source lang=matlab:)
function obj = f(x)
x1 = x(1);
x2 = x(2);
obj = 0.2 + x1^2 + x2^2 - 0.1*cos(6.0*pi*x1) - 0.1*cos(6.0*pi*x2);
end
(:sourceend:)
(:divend:)
Changed line 7 from:
Attach:sim_annealing_boltzmann.png
to:
{$P = \exp\left(\frac{-\Delta E}{k_b T}\right)$}
Changed lines 21-27 from:
You can download anneal.m and anneal.py files to retrieve example simulated annealing files in MATLAB and Python, respectively.
to:
!!!! Example Problem and Source Code

Find the minimum to the objective function

{$obj = 0.2 + x_1^2 + x_2^2 - 0.1 \, \cos \left( 6 \pi x_1 \right) - 0.1 \cos \left(6 \pi x_2\right)$}

by adjusting the values of {`x_1`} and {`x_2`}.
You can download anneal.m and anneal.py files to retrieve example simulated annealing files in MATLAB and Python, respectively.
May 31, 2018, at 09:03 PM by 10.5.113.123 -
Deleted lines 3-4:

!!!! Simulated Annealing
Changed line 11 from:
where  is the Boltzmann constant and T is the current temperature. By examining this equation we should note two things: the probability is proportional to temperature--as the solid cools, the probability gets smaller; and inversely proportional to  --as the change in energy is larger the probability of accepting the change gets smaller.
to:
where k'_b_' is the Boltzmann constant and T is the current temperature. By examining this equation we should note two things: the probability is proportional to temperature--as the solid cools, the probability gets smaller; and inversely proportional to  --as the change in energy is larger the probability of accepting the change gets smaller.
March 02, 2013, at 09:28 PM by 69.169.188.188 -
Changed lines 23-31 from:
Click the link below to use a Google Docs word processor to interactively and collaboratively develop the MATLAB or Python source code. These documents will only be unlocked for editing during class interactive sessions but they can be viewed at any time.

* [[https://docs.google.com/document/d/1pYurFx302tFh3J-i3O0yR7nFauVqiUPUu5eLJLHIKtc/edit?usp=sharing|MATLAB Collaborative Source Code]]
* [[https://docs.google.com/document/d/1oSdU7tETNkgXqpyEzuwNGPj6UkoeJWa51Sz1sSXs0pk/edit?usp=sharing|Python Collaborative Source Code]]

----

You can also download these files to retrieve the latest versions of the example simulated annealing files in MATLAB and Python
.
to:
You can download anneal.m and anneal.py files to retrieve example simulated annealing files in MATLAB and Python, respectively.
Deleted lines 26-27:

Please leave comments below if you would like to suggest improvements to the current code.
March 02, 2013, at 12:39 AM by 69.169.131.173 -
Changed line 7 from:
Simulated annealing copies a phenomenon in nature--the annealing of solids--to optimize a complex system.  Annealing refers to heating a solid and then cooling it slowly so that thermal equilibrium is maintained.  Atoms then assume a nearly globally minimum energy state.  In 1953 Metropolis created an algorithm to simulate the annealing process.  The algorithm simulates a small random displacement of an atom that results in a change in energy.  If the change in energy is negative, the energy state of the new configuration is lower and the new configuration is accepted.  If the change in energy is positive, the new configuration has a higher energy state; however, it may still be accepted according to the Boltzmann probability factor:
to:
Simulated annealing copies a phenomenon in nature--the annealing of solids--to optimize a complex system.  Annealing refers to heating a solid and then cooling it slowly.  Atoms then assume a nearly globally minimum energy state.  In 1953 Metropolis created an algorithm to simulate the annealing process.  The algorithm simulates a small random displacement of an atom that results in a change in energy.  If the change in energy is negative, the energy state of the new configuration is lower and the new configuration is accepted.  If the change in energy is positive, the new configuration has a higher energy state; however, it may still be accepted according to the Boltzmann probability factor:
March 02, 2013, at 12:37 AM by 69.169.131.173 -
Changed line 7 from:
Simulated annealing copies a phenomenon in nature--the annealing of solids--to optimize a complex system.  Annealing refers to heating a solid to a liquid state and then cooling it slowly so that thermal equilibrium is maintained.  Atoms then assume a nearly globally minimum energy state.  In 1953 Metropolis created an algorithm to simulate the annealing process.  The algorithm simulates a small random displacement of an atom that results in a change in energy.  If the change in energy is negative, the energy state of the new configuration is lower and the new configuration is accepted.  If the change in energy is positive, the new configuration has a higher energy state; however, it may still be accepted according to the Boltzmann probability factor:
to:
Simulated annealing copies a phenomenon in nature--the annealing of solids--to optimize a complex system.  Annealing refers to heating a solid and then cooling it slowly so that thermal equilibrium is maintained.  Atoms then assume a nearly globally minimum energy state.  In 1953 Metropolis created an algorithm to simulate the annealing process.  The algorithm simulates a small random displacement of an atom that results in a change in energy.  If the change in energy is negative, the energy state of the new configuration is lower and the new configuration is accepted.  If the change in energy is positive, the new configuration has a higher energy state; however, it may still be accepted according to the Boltzmann probability factor:
March 01, 2013, at 08:23 PM by 128.187.97.21 -
Added lines 36-39:

Attach:sim_annealing_contour.png

Attach:sim_annealing_iterations.png
March 01, 2013, at 08:18 PM by 128.187.97.21 -
Added lines 14-15:

Attach:sim_annealing_flow.png
March 01, 2013, at 08:14 PM by 128.187.97.21 -
Changed lines 21-22 from:
Click the link below to use a Google Docs word processor to interactively and collaboratively develop the MATLAB or Python source code.
to:
Click the link below to use a Google Docs word processor to interactively and collaboratively develop the MATLAB or Python source code. These documents will only be unlocked for editing during class interactive sessions but they can be viewed at any time.
Deleted line 24:
March 01, 2013, at 08:12 PM by 128.187.97.21 -
Changed line 23 from:
* [[https://docs.google.com/document/d/1oSdU7tETNkgXqpyEzuwNGPj6UkoeJWa51Sz1sSXs0pk/edit?usp=sharing|MATLAB Collaborative Source Code]]
to:
* [[https://docs.google.com/document/d/1pYurFx302tFh3J-i3O0yR7nFauVqiUPUu5eLJLHIKtc/edit?usp=sharing|MATLAB Collaborative Source Code]]
March 01, 2013, at 08:09 PM by 128.187.97.21 -
Changed line 32 from:
[[Attach:simulated_annealing.zip|simulated_annealing.png]]
to:
[[Attach:simulated_annealing.zip|Attach:simulated_annealing.png]]
March 01, 2013, at 08:06 PM by 128.187.97.21 -
Added lines 1-53:
(:title Simulated Annealing Tutorial:)
(:keywords simulated annealing, Python, MATLAB:)
(:description Tutorial on simulated annealing in optimization applications.:)

!!!! Simulated Annealing

Simulated annealing copies a phenomenon in nature--the annealing of solids--to optimize a complex system.  Annealing refers to heating a solid to a liquid state and then cooling it slowly so that thermal equilibrium is maintained.  Atoms then assume a nearly globally minimum energy state.  In 1953 Metropolis created an algorithm to simulate the annealing process.  The algorithm simulates a small random displacement of an atom that results in a change in energy.  If the change in energy is negative, the energy state of the new configuration is lower and the new configuration is accepted.  If the change in energy is positive, the new configuration has a higher energy state; however, it may still be accepted according to the Boltzmann probability factor:

Attach:sim_annealing_boltzmann.png
   
where  is the Boltzmann constant and T is the current temperature. By examining this equation we should note two things: the probability is proportional to temperature--as the solid cools, the probability gets smaller; and inversely proportional to  --as the change in energy is larger the probability of accepting the change gets smaller.

When applied to engineering design, an analogy is made between energy and the objective function.  The design is started at a high β€œtemperature”, where it has a high objective (we assume we are minimizing).  Random perturbations are then made to the design.  If the objective is lower, the new design is made the current design; if it is higher, it may still be accepted according the probability given by the Boltzmann factor.  The Boltzmann probability is compared to a random number drawn from a uniform distribution between 0 and 1; if the random number is smaller than the Boltzmann probability, the configuration is accepted. This allows the algorithm to escape local minima.

As the temperature is gradually lowered, the probability that a worse design is accepted becomes smaller.  Typically at high temperatures the gross structure of the design emerges which is then refined at lower temperatures. 

Although it can be used for continuous problems, simulated annealing is especially effective when applied to combinatorial or discrete problems.  Although the algorithm is not guaranteed to find the best optimum, it will often find near optimum designs with many fewer design evaluations than other algorithms.  (It can still be computationally expensive, however.)  It is also an easy algorithm to implement.

----

Click the link below to use a Google Docs word processor to interactively and collaboratively develop the MATLAB or Python source code.

* [[https://docs.google.com/document/d/1oSdU7tETNkgXqpyEzuwNGPj6UkoeJWa51Sz1sSXs0pk/edit?usp=sharing|MATLAB Collaborative Source Code]]
* [[https://docs.google.com/document/d/1oSdU7tETNkgXqpyEzuwNGPj6UkoeJWa51Sz1sSXs0pk/edit?usp=sharing|Python Collaborative Source Code]]


----

You can also download these files to retrieve the latest versions of the example simulated annealing files in MATLAB and Python.

[[Attach:simulated_annealing.zip|Download Simulated Annealing Example Files]]
[[Attach:simulated_annealing.zip|simulated_annealing.png]]

Please leave comments below if you would like to suggest improvements to the current code.

----

(:html:)
 <div id="disqus_thread"></div>
    <script type="text/javascript">
        /* * * CONFIGURATION VARIABLES: EDIT BEFORE PASTING INTO YOUR WEBPAGE * * */
        var disqus_shortname = 'apmonitor'; // required: replace example with your forum shortname

        /* * * DON'T EDIT BELOW THIS LINE * * */
        (function() {
            var dsq = document.createElement('script'); dsq.type = 'text/javascript'; dsq.async = true;
            dsq.src = 'https://' + disqus_shortname + '.disqus.com/embed.js';
            (document.getElementsByTagName('head')[0] || document.getElementsByTagName('body')[0]).appendChild(dsq);
        })();
    </script>
    <noscript>Please enable JavaScript to view the <a href="https://disqus.com/?ref_noscript">comments powered by Disqus.</a></noscript>
    <a href="https://disqus.com" class="dsq-brlink">comments powered by <span class="logo-disqus">Disqus</span></a>
(:htmlend:)