W kilku projektach, nad którymi pracowałem przyszło mi użyć już kilka razy algorytmu Symulowanego Wyżarzania (ang. Simulated Annealing). Symulowane wyżarzanie to rodzaj algorytmu heurystycznego, przeszukującego przestrzeń alternatywnych rozwiązań problemu w celu wyszukania rozwiązań najlepszych. Sposób działania symulowanego wyżarzania przypomina zjawisko wyżarzania w metalurgii.

Z racji dość uniwersalnego zastosowania, w nowym projekcie postanowiłem zdefiniować problem jako interfejs, i zaimplementować generalne rozwiązanie dla niego.

Na listingu poniżej klasa SimulatedAnnealing implementująca wspomniany algorytm:

[java]
package com.konri.science.sa;

import com.konri.science.utils.IProgressListener;

public class SimulatedAnnealing {

private ISimulatedAnnealingProblem problem;
private IProgressListener progressListener;

private double maxTemperature = 1000.0;
private double maxInternalSteps = 100.0;
private double temperatureDelta = 0.95;

public SimulatedAnnealing(ISimulatedAnnealingProblem problem) {
this.problem = problem;
}

public SimulatedAnnealing(ISimulatedAnnealingProblem problem, IProgressListener progressListener) {
this(problem);
this.progressListener = progressListener;
}

public SimulatedAnnealing(ISimulatedAnnealingProblem problem, double maxTemperature, double maxInternalSteps, double temperatureDelta) {
this(problem);
this.maxTemperature = maxTemperature;
this.maxInternalSteps = maxInternalSteps;
this.temperatureDelta = temperatureDelta;
}

public SimulatedAnnealing(ISimulatedAnnealingProblem problem, double maxTemperature, double maxInternalSteps, double temperatureDelta, IProgressListener progressListener) {
this(problem, maxTemperature, maxInternalSteps, temperatureDelta);
this.progressListener = progressListener;
}

@SuppressWarnings({“MethodWithMultipleLoops”})
public void solveProblem() {
double temperature = maxTemperature;
double minSolutionValue = problem.getSolutionValue();
do {
int innerCounter = 0;
do {
problem.makeCopyOfSolution();
problem.randomChangeOfSolution(temperature, maxTemperature);
double currentSolutionValue = problem.getSolutionValue();

if (currentSolutionValue < minSolutionValue)
minSolutionValue = currentSolutionValue;
else {
if (shouldTryWorseSolution(currentSolutionValue, minSolutionValue, temperature))
minSolutionValue = currentSolutionValue;
else
problem.restorePreviousSolution();
}

innerCounter++;
} while (innerCounter < (int) maxInternalSteps);
temperature *= temperatureDelta;
if (progressListener != null)
progressListener.progress(1 - temperature / maxTemperature, temperature);
} while (temperature > 0.01);
}

private static boolean shouldTryWorseSolution(double currentSolutionValue, double minSolutionValue, double temperature) {
double p = StrictMath.exp(-1.0 * ((currentSolutionValue – minSolutionValue) / temperature));
return StrictMath.random() > p;
}
}
[/java]

I jeszcze interfejs problemu:

[java]
package com.konri.science.sa;

public interface ISimulatedAnnealingProblem {

double getSolutionValue();

void randomChangeOfSolution(double currentTemperature, double maxTemperature);

void makeCopyOfSolution();

void restorePreviousSolution();
}
[/java]

Zdefiniowałem również interfejs pozwalający wyświetlać postęp algorytmu, lecz dla celów optymalizacji można go usunąć z kodu.

[java]
package com.konri.science.utils;

public interface IProgressListener {
void progress(double percent, double temperature);
}
[/java]

Przydał by się oczywiście jakiś test:

[java]
package com.konri.science.sa;

import junit.framework.TestCase;

public class SimulatedAnnealingTest extends TestCase {

private static final double SOLUTION_DELTA = 0.0001;

private MockSimpleSimulatedAnnealingProblem simpleProblem;
private MockComplexSimulatedAnnealingProblem complexProblem;
private SimulatedAnnealing simulatedAnnealing;

protected void setUp() throws Exception {
simpleProblem = new MockSimpleSimulatedAnnealingProblem();
complexProblem = new MockComplexSimulatedAnnealingProblem();

}

public void testSimpleSimulatedAnnealing() {
simulatedAnnealing = new SimulatedAnnealing(simpleProblem);
simpleProblem.setTargetNumber(3);
simulatedAnnealing.solveProblem();
assertEquals(simpleProblem.getTargetNumber(), simpleProblem.getCurrentNumber(), SOLUTION_DELTA);

simpleProblem.setTargetNumber(-3);
simulatedAnnealing.solveProblem();
assertEquals(simpleProblem.getTargetNumber(), simpleProblem.getCurrentNumber(), SOLUTION_DELTA);

simpleProblem.setTargetNumber(0);
simulatedAnnealing.solveProblem();
assertEquals(simpleProblem.getTargetNumber(), simpleProblem.getCurrentNumber(), SOLUTION_DELTA);
}

public void testComplexSimulatedAnnealing() {
simulatedAnnealing = new SimulatedAnnealing(complexProblem);
simulatedAnnealing.solveProblem();
assertEquals(0, complexProblem.getSolutionValue(), SOLUTION_DELTA);
}
}
[/java]

Oto dwa problemy testowe:

[java]
package com.konri.science.sa;

public class MockSimpleSimulatedAnnealingProblem implements ISimulatedAnnealingProblem {

private static final double NEGATIVE_FACTOR = 0.5;
private static final double MAXIMUM_CHANGE = 1.0;

private double targetNumber;
private double currentNumber = Math.random() * 2 – 1;
private double coppyOfSolution = currentNumber;

public final double getSolutionValue() {
return Math.abs(targetNumber – currentNumber);
}

public void randomChangeOfSolution(double currentTemperature, double maxTemperature) {
currentNumber += calculateChange(MAXIMUM_CHANGE) * (currentTemperature / maxTemperature);
}

public void makeCopyOfSolution() {
coppyOfSolution = currentNumber;
}

private static double calculateChange(double maximumChange) {
return (Math.random() – NEGATIVE_FACTOR) * maximumChange;
}

public void restorePreviousSolution() {
currentNumber = coppyOfSolution;
}

public final double getTargetNumber() {
return targetNumber;
}

public void setTargetNumber(double targetNumber) {
this.targetNumber = targetNumber;
}

public final double getCurrentNumber() {
return currentNumber;
}
}
[/java]

[java]
package com.konri.science.sa;

public class MockComplexSimulatedAnnealingProblem implements ISimulatedAnnealingProblem {

private static final double NEGATIVE_FACTOR = 0.5;
private static final double MAXIMUM_CHANGE = 1.0;

private static final double[][] input = {
{4.0, 7.0, 6.0},
{0.0, 8.0, 7.0},
{6.0, 6.0, 1.0},
{3.0, 9.0, 8.0},
{1.0, 9.0, 3.0}};

private static final double[] output = {18.600, 11.200, 8.600, 19.900, -0.700};

//From spreadsheet: 2.3; -1.4; 3.2
private final double[] params = {0.0, 0.0, 0.0};
private final double[] paramsCopy = {0.0, 0.0, 0.0};

public final double getSolutionValue() {
return calculateSquerError();
}

public void randomChangeOfSolution(double currentTemperature, double maxTemperature) {
int index = (int) StrictMath.round(StrictMath.random() * 2.0);
params[index] += calculateChange(MAXIMUM_CHANGE) * (currentTemperature / maxTemperature);
}

private static double calculateChange(double maximumChange) {
return (StrictMath.random() – NEGATIVE_FACTOR) * maximumChange;
}

public void restorePreviousSolution() {
copyParamsArray(paramsCopy, params);
}

public void makeCopyOfSolution() {
copyParamsArray(params, paramsCopy);
}

private void copyParamsArray(double[] source, double[] destination) {
System.arraycopy(source, 0, destination, 0, params.length);
}

private double calculateSquerError() {
double sum = 0.0;
for (int i = 0; i < 5; i++)
sum += calculateFunctionSqrError(i);

if (sum > 0.0)
return Math.sqrt(sum);
else
return 0.0;
}

private double calculateFunctionSqrError(int i) {
return StrictMath.pow(calculateFunctionError(i), 2.0);
}

private double calculateFunctionError(int i) {
return output[i] – calculateFunction(i);
}

private double calculateFunction(int i) {
return calculateOutput(i, 0) + calculateOutput(i, 1) + calculateOutput(i, 2);
}

private double calculateOutput(int i, int item) {
return params[item] * input[i][item];
}

public final double[] getParams() {
return params;
}
}
[/java]

Poniżej znajdują się natomiast kody źródłowe przykładowego zastosowania. Program SampleSaProblem poszukuje postaci wcześniej zdefiniowanej funkcji na podstawie znanych wektorów x i fx.
Być może nie są one zaimplementowane optymalnie, ale służą jedynie jako poglądowe przykłady:

[java]
package com.konri.apps;

import com.konri.science.sa.SimulatedAnnealing;
import com.konri.science.utils.IProgressListener;

import java.text.MessageFormat;

@SuppressWarnings({“AssignmentToCollectionOrArrayFieldFromParameter”})
public class SampleSaProblem {

public static void main(String[] args) {

Function function = new Function(
new double[]{0.3, 0.75, 1.3, 2.2},
new double[]{1.3, 1.75, 2.3}
);

double[] x = new double[400];
double[] fx = new double[400];

int i = 0;
for (double v = 0.0; v < 100.0; v += 0.25) {
x[i] = v;
fx[i] = function.eval(v);
i++;
}

System.out.println("+------------------------------------------------------------------------------");
System.out.println(MessageFormat.format("| Original f(x) = {0}", function));
System.out.println("+------------------------------------------------------------------------------");

final FunctionFinder functionFinder = new FunctionFinder(x, fx, 3, 1.0, 5.0, 3.0);
double svb = functionFinder.getSolutionValue();
System.out.println(MessageFormat.format("| Sqr. diff. before: {0}", svb));
System.out.println(MessageFormat.format("| F: {0}", functionFinder.getF()));

SimulatedAnnealing simulatedAnnealing = new SimulatedAnnealing(functionFinder, 100000000.0, 1000000.0, 0.99, new IProgressListener() {

private int ct;

@Override
public void progress(double percent, double temperature) {
ct++;
if (ct > 9) {
System.out.println(MessageFormat.format(“| dT: {0,number,00}%, SV: {1}, T: {2}, Ch: {3}”, percent * 100.0, functionFinder.getSolutionValue(), temperature, functionFinder.getChangesCt()));
ct = 0;
}
}
});

long timeStart = System.currentTimeMillis();
simulatedAnnealing.solveProblem();

double sva = functionFinder.getSolutionValue();
System.out.println(MessageFormat.format(“| Sqr. diff. after : {0}”, sva));
System.out.println(MessageFormat.format(“| F: {0}”, functionFinder.getF()));
System.out.println(MessageFormat.format(“| Time: {0} sec. Result: {1}”, (System.currentTimeMillis() – timeStart) / 1000L, sva / svb));
}
}
[/java]

[java]
package com.konri.apps;

import com.konri.science.sa.ISimulatedAnnealingProblem;

import java.util.Random;

@SuppressWarnings({“AssignmentToCollectionOrArrayFieldFromParameter”})
public class FunctionFinder implements ISimulatedAnnealingProblem {

private static final double MAXIMUM_CHANGE = 15.0;
private static final double NEGATIVE_FACTOR = 0.5;

private static final Random RANDOM = new Random();
private double maximumChange = MAXIMUM_CHANGE;

private long changesCt;

private Function f;
private Function fcopy;

private double[] x;
private double[] fx;

public FunctionFinder(double[] x, double[] fx, int n, double maximumChange, double maxw, double maxp) {
this.x = x;
this.fx = fx;
f = new Function(n, maxw, maxp);

fcopy = f.copy();
this.maximumChange = maximumChange;
}

public FunctionFinder(double[] x, double[] fx, int n, double maxw, double maxp) {
this(x, fx, n, MAXIMUM_CHANGE, maxw, maxp);
}

@Override
public double getSolutionValue() {
double delta = 0.0;
for (int i = 0; i < x.length; i++) {
double v = x[i];
double nfx = f.eval(v);
delta += (fx[i] - nfx) * (fx[i] - nfx);
}
return StrictMath.sqrt(delta / (double) x.length);
}

@Override
public void randomChangeOfSolution(double currentTemperature, double maxTemperature) {
changesCt++;
double delta = calculateChange(maximumChange) * (currentTemperature / maxTemperature);
if (RANDOM.nextBoolean()) {
int rand = RANDOM.nextInt(f.deg() + 1);
f.setW(rand, f.getW(rand) + delta);
} else {
int rand = RANDOM.nextInt(f.deg());
f.setP(rand, f.getP(rand) + delta);
}
}

public long getChangesCt() {
return changesCt;
}

@Override
public void makeCopyOfSolution() {
fcopy = f.copy();
}

@Override
public void restorePreviousSolution() {
f = fcopy.copy();
}

private static double calculateChange(double maximumChange) {
return (RANDOM.nextDouble() - NEGATIVE_FACTOR) * maximumChange;
}

public Function getF() {
return f;
}
}
[/java]

[java]
package com.konri.apps;

import java.text.MessageFormat;

public class Function {

private double[] w;
private double[] p;

public Function(int n, double maxw, double maxp) {
if (n <= 0)
throw new IllegalArgumentException("N must be grater than 0");

w = new double[n + 1];
p = new double[n];

for (int i = 0; i <= n; i++) {
w[i] = (StrictMath.random() - 0.5) * 2.0 * maxw;
if (i < n)
p[i] = (StrictMath.random() - 0.5) * 2.0 * maxp;
}
}

public Function(double[] w, double[] p) {
this(w.length - 1, 5.0, 3.0);
System.arraycopy(w, 0, this.w, 0, w.length);
System.arraycopy(p, 0, this.p, 0, p.length);
}

@Override
public String toString() {
StringBuilder res = new StringBuilder("f(x) = ");
res.append(format(getW(w.length - 1)));
res.append("x^");
res.append(format(getP(w.length - 2)));

for (int i = w.length - 2; i > 0; i–) {
if (getW(i) < 0.0) {
res.append(" - ");
res.append(format(-getW(i)));
} else {
res.append(" + ");
res.append(format(getW(i)));
}
res.append("x^").append(format(getP(i - 1)));
}

if (getW(0) < 0.0) {
res.append(" - ");
res.append(format(-getW(0)));
} else if (getW(0) > 0.0) {
res.append(” + “);
res.append(format(getW(0)));
}
return res.toString();
}

private static String format(double v) {
return MessageFormat.format(“{0,number,0.00}”, v);
}

public double eval(double x) {
double y = w[deg()];
for (int i = deg() – 1; i >= 0; i–)
y = y * x + w[i];

return y;
}

public int deg() {
return w.length – 1;
}

public double getW(int i) {
return w[i];
}

public double getP(int i) {
return p[i];
}

public void setW(int i, double x) {
w[i] = x;
}

public void setP(int i, double x) {
p[i] = x;
}

public Function copy() {
return new Function(w, p);
}
}
[/java]

Powyższy kod skompilujesz pod Java 1.6, w przypadku wersji 1.5 musisz usunąć anotacje @Override z metod implementujących interfejsy.
W razie czego proszę dać znać o błędach