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()