001package org.cpsolver.ifs.algorithms; 002 003import java.text.DecimalFormat; 004 005import org.cpsolver.ifs.assignment.Assignment; 006import org.cpsolver.ifs.heuristics.NeighbourSelection; 007import org.cpsolver.ifs.model.Model; 008import org.cpsolver.ifs.model.Neighbour; 009import org.cpsolver.ifs.model.Value; 010import org.cpsolver.ifs.model.Variable; 011import org.cpsolver.ifs.solution.Solution; 012import org.cpsolver.ifs.util.DataProperties; 013import org.cpsolver.ifs.util.JProf; 014import org.cpsolver.ifs.util.ToolBox; 015 016 017/** 018 * Simulated annealing. In each iteration, one of the given neighbourhoods is selected first, 019 * then a neighbour is generated and it is accepted with probability 020 * {@link SimulatedAnnealingContext#prob(double)}. The search is guided by the 021 * temperature, which starts at <i>SimulatedAnnealing.InitialTemperature</i>. 022 * After each <i>SimulatedAnnealing.TemperatureLength</i> iterations, the 023 * temperature is reduced by <i>SimulatedAnnealing.CoolingRate</i>. If there was 024 * no improvement in the past <i>SimulatedAnnealing.ReheatLengthCoef * 025 * SimulatedAnnealing.TemperatureLength</i> iterations, the temperature is 026 * increased by <i>SimulatedAnnealing.ReheatRate</i>. If there was no 027 * improvement in the past <i>SimulatedAnnealing.RestoreBestLengthCoef * 028 * SimulatedAnnealing.TemperatureLength</i> iterations, the best ever found 029 * solution is restored. <br> 030 * <br> 031 * If <i>SimulatedAnnealing.StochasticHC</i> is true, the acceptance probability 032 * is computed using stochastic hill climbing criterion, i.e., 033 * <code>1.0 / (1.0 + Math.exp(value/temperature))</code>, otherwise it is 034 * cumputed using simlated annealing criterion, i.e., 035 * <code>(value<=0.0?1.0:Math.exp(-value/temperature))</code>. If 036 * <i>SimulatedAnnealing.RelativeAcceptance</i> neighbour value 037 * {@link Neighbour#value(Assignment)} is taken as the value of the selected 038 * neighbour (difference between the new and the current solution, if the 039 * neighbour is accepted), otherwise the value is computed as the difference 040 * between the value of the current solution if the neighbour is accepted and 041 * the best ever found solution. <br> 042 * <br> 043 * Custom neighbours can be set using SimulatedAnnealing.Neighbours property that should 044 * contain semicolon separated list of {@link NeighbourSelection}. By default, 045 * each neighbour selection is selected with the same probability (each has 1 point in 046 * a roulette wheel selection). It can be changed by adding @n at the end 047 * of the name of the class, for example: 048 * <pre><code> 049 * SimulatedAnnealing.Neighbours=org.cpsolver.ifs.algorithms.neighbourhoods.RandomMove;org.cpsolver.ifs.algorithms.neighbourhoods.RandomSwapMove@0.1 050 * </code></pre> 051 * Selector RandomSwapMove is 10× less probable to be selected than other selectors. 052 * When SimulatedAnnealing.Random is true, all selectors are selected with the same probability, ignoring these weights. 053 * <br><br> 054 * When SimulatedAnnealing.Update is true, {@link NeighbourSelector#update(Assignment, Neighbour, long)} is called 055 * after each iteration (on the selector that was used) and roulette wheel selection 056 * that is using {@link NeighbourSelector#getPoints()} is used to pick a selector in each iteration. 057 * See {@link NeighbourSelector} for more details. 058 * <br> 059 * 060 * @version IFS 1.3 (Iterative Forward Search)<br> 061 * Copyright (C) 2014 Tomáš Müller<br> 062 * <a href="mailto:muller@unitime.org">muller@unitime.org</a><br> 063 * <a href="http://muller.unitime.org">http://muller.unitime.org</a><br> 064 * <br> 065 * This library is free software; you can redistribute it and/or modify 066 * it under the terms of the GNU Lesser General Public License as 067 * published by the Free Software Foundation; either version 3 of the 068 * License, or (at your option) any later version. <br> 069 * <br> 070 * This library is distributed in the hope that it will be useful, but 071 * WITHOUT ANY WARRANTY; without even the implied warranty of 072 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 073 * Lesser General Public License for more details. <br> 074 * <br> 075 * You should have received a copy of the GNU Lesser General Public 076 * License along with this library; if not see 077 * <a href='http://www.gnu.org/licenses/'>http://www.gnu.org/licenses/</a>. 078 * @param <V> Variable 079 * @param <T> Value 080 */ 081public class SimulatedAnnealing<V extends Variable<V, T>, T extends Value<V, T>> extends NeighbourSearch<V,T> { 082 private DecimalFormat iDF5 = new DecimalFormat("0.00000"); 083 private DecimalFormat iDF10 = new DecimalFormat("0.0000000000"); 084 private double iInitialTemperature = -1; 085 private double iMaximalTemperature = 1.5; 086 private double iMinimalTemperature = 0.001; 087 private double iCoolingRate = 0.95; 088 private double iReheatRate = -1; 089 private long iTemperatureLength = 25000; 090 private double iReheatLengthCoef = 5.0; 091 private double iRestoreBestLengthCoef = -1; 092 private boolean iStochasticHC = false; 093 private boolean iRelativeAcceptance = true; 094 private Double[] iCoolingRateAdjusts = null; 095 private int iTrainingValues = 10000; 096 private double iTrainingProbability = 0.00001; 097 private double iTimeBetweenCooldowns = 10.0; 098 099 /** 100 * Constructor. Following problem properties are considered: 101 * <ul> 102 * <li>SimulatedAnnealing.InitialTemperature ... initial temperature (default 1.5) 103 * <li>SimulatedAnnealing.TemperatureLength ... temperature length (number of iterations between temperature decrements, default 25000) 104 * <li>SimulatedAnnealing.CoolingRate ... temperature cooling rate (default 0.95) 105 * <li>SimulatedAnnealing.ReheatLengthCoef ... temperature re-heat length coefficient (multiple of temperature length, default 5) 106 * <li>SimulatedAnnealing.ReheatRate ... temperature re-heating rate (default (1/coolingRate)^(reheatLengthCoef*1.7)) 107 * <li>SimulatedAnnealing.RestoreBestLengthCoef ... restore best length coefficient (multiple of re-heat length, default reheatLengthCoef^2) 108 * <li>SimulatedAnnealing.StochasticHC ... true for stochastic search acceptance criterion, false for simulated annealing acceptance (default false) 109 * <li>SimulatedAnnealing.RelativeAcceptance ... true for relative acceptance (different between the new and the current solutions, if the neighbour is accepted), false for absolute acceptance (difference between the new and the best solutions, if the neighbour is accepted) 110 * <li>SimulatedAnnealing.Neighbours ... semicolon separated list of classes implementing {@link NeighbourSelection} 111 * <li>SimulatedAnnealing.AdditionalNeighbours ... semicolon separated list of classes implementing {@link NeighbourSelection} 112 * <li>SimulatedAnnealing.Random ... when true, a neighbour selector is selected randomly 113 * <li>SimulatedAnnealing.Update ... when true, a neighbour selector is selected using {@link NeighbourSelector#getPoints()} weights (roulette wheel selection) 114 * </ul> 115 * 116 * @param properties 117 * problem properties 118 */ 119 public SimulatedAnnealing(DataProperties properties) { 120 super(properties); 121 iInitialTemperature = properties.getPropertyDouble(getParameterBaseName() + ".InitialTemperature", iInitialTemperature); 122 iMaximalTemperature = properties.getPropertyDouble(getParameterBaseName() + ".MaximalTemperature", (iInitialTemperature <= 0.0 ? -1.0 : iMaximalTemperature)); 123 iMinimalTemperature = properties.getPropertyDouble(getParameterBaseName() + ".MinimalTemperature", (iInitialTemperature <= 0.0 ? -1.0 : 0.0)); 124 iReheatRate = properties.getPropertyDouble(getParameterBaseName() + ".ReheatRate", iReheatRate); 125 iCoolingRate = properties.getPropertyDouble(getParameterBaseName() + ".CoolingRate", iCoolingRate); 126 iRelativeAcceptance = properties.getPropertyBoolean(getParameterBaseName() + ".RelativeAcceptance", iRelativeAcceptance); 127 iStochasticHC = properties.getPropertyBoolean(getParameterBaseName() + ".StochasticHC", iStochasticHC); 128 iTemperatureLength = properties.getPropertyLong(getParameterBaseName() + ".TemperatureLength", iTemperatureLength); 129 iReheatLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".ReheatLengthCoef", iReheatLengthCoef); 130 iRestoreBestLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".RestoreBestLengthCoef", iRestoreBestLengthCoef); 131 iCoolingRateAdjusts = properties.getPropertyDoubleArry(getParameterBaseName() + ".CoolingRateAdjustments", null); 132 iTrainingValues = properties.getPropertyInt(getParameterBaseName() + ".TrainingValues", iTrainingValues); 133 iTrainingProbability = properties.getPropertyDouble(getParameterBaseName() + ".TrainingProbability", iTrainingProbability); 134 iTimeBetweenCooldowns = properties.getPropertyDouble(getParameterBaseName() + ".TimeBetweenCooldowns", iTimeBetweenCooldowns); 135 if (iReheatRate < 0) 136 iReheatRate = Math.pow(1 / iCoolingRate, iReheatLengthCoef * 1.7); 137 if (iRestoreBestLengthCoef < 0) 138 iRestoreBestLengthCoef = iReheatLengthCoef * iReheatLengthCoef; 139 } 140 141 @Override 142 public String getParameterBaseName() { return "SimulatedAnnealing"; } 143 144 @Override 145 public NeighbourSearchContext createAssignmentContext(Assignment<V, T> assignment) { 146 return new SimulatedAnnealingContext(); 147 } 148 149 public class SimulatedAnnealingContext extends NeighbourSearchContext { 150 private double iTemperature = 0.0; 151 private int iMoves = 0; 152 private double iAbsValue = 0; 153 private double iBestValue = 0; 154 private long iLastImprovingIter = -1; 155 private long iLastBestIter = -1; 156 private long iLastReheatIter = 0; 157 private long iLastCoolingIter = 0; 158 private int iAcceptIter[] = new int[] { 0, 0, 0 }; 159 private long iReheatLength = 0; 160 private long iRestoreBestLength = 0; 161 private int iTrainingIterations = 0; 162 private double iTrainingTotal = 0.0; 163 164 /** Setup the temperature */ 165 @Override 166 protected void activate(Solution<V, T> solution) { 167 super.activate(solution); 168 iTrainingTotal = 0.0; iTrainingIterations = 0; 169 iTemperature = iInitialTemperature; 170 iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength); 171 iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength); 172 iLastImprovingIter = -1; 173 iLastBestIter = -1; 174 } 175 176 protected double getCoolingRate(int idx) { 177 if (idx < 0 || iCoolingRateAdjusts == null || idx >= iCoolingRateAdjusts.length || iCoolingRateAdjusts[idx] == null) return iCoolingRate; 178 return iCoolingRate * iCoolingRateAdjusts[idx]; 179 } 180 181 /** 182 * Set initial temperature based on the training period 183 * @param solution current solution 184 */ 185 protected void train(Solution<V, T> solution) { 186 double value = iTrainingTotal / iTrainingIterations; 187 if (iStochasticHC) { 188 iInitialTemperature = value / Math.log(1.0 / 0.01 * iTrainingProbability - 1.0); 189 if (iMaximalTemperature <= 0.0) 190 iMaximalTemperature = value / Math.log(1.0 / iTrainingProbability - 1.0); 191 if (iMinimalTemperature < 0.0) 192 iMinimalTemperature = 0.1 / Math.log(1.0 / iTrainingProbability - 1.0); 193 } else { 194 iInitialTemperature = - value / Math.log(0.01 * iTrainingProbability); 195 if (iMaximalTemperature <= 0.0) 196 iMaximalTemperature = - value / Math.log(iTrainingProbability); 197 if (iMinimalTemperature < 0.0) 198 iMinimalTemperature = - 0.1 / Math.log(iTrainingProbability); 199 } 200 iTemperature = iInitialTemperature; 201 info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)) 202 + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " + 203 "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) + 204 ", Best=" + iDF2.format(solution.getBestValue()) + 205 " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)"); 206 info("Temperature set to " + iDF5.format(iTemperature) + " " + "(" + 207 "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " + 208 "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " + 209 "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " + 210 "p(+" + iDF2.format(value) + ")=" + iDF5.format(100.0 * prob(value)) + "%)"); 211 info("Maximal temperature set to " + iDF5.format(iMaximalTemperature) + " " + "(" + 212 "p(+0.1)=" + iDF2.format(100.0 * prob(0.1, iMaximalTemperature)) + "%, " + 213 "p(+1)=" + iDF2.format(100.0 * prob(1, iMaximalTemperature)) + "%, " + 214 "p(+10)=" + iDF5.format(100.0 * prob(10, iMaximalTemperature)) + "%, " + 215 "p(+" + iDF2.format(value) + ")=" + iDF5.format(100.0 * prob(value, iMaximalTemperature)) + "%)"); 216 info("Minimal temperature set to " + iDF5.format(iMinimalTemperature) + " " + "(" + 217 "p(+0.001)=" + iDF2.format(100.0 * prob(0.001, iMinimalTemperature)) + "%, " + 218 "p(+0.01)=" + iDF2.format(100.0 * prob(0.01, iMinimalTemperature)) + "%, " + 219 "p(+0.1)=" + iDF5.format(100.0 * prob(0.1, iMinimalTemperature)) + "%)"); 220 logNeibourStatus(); 221 double speed = ((double)iIter) / (JProf.currentTimeMillis() - iT0); // iterations per ms 222 iTemperatureLength = Math.round(speed * iTimeBetweenCooldowns * 1000.0); 223 iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength); 224 iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength); 225 info("Training speed was " + iDF2.format(1000.0 * speed) + " it/s, temperature length adjusted to " + iTemperatureLength + "."); 226 iIter = 0; iT0 = JProf.currentTimeMillis(); 227 iLastImprovingIter = -1; iLastBestIter = 0; iBestValue = solution.getBestValue(); 228 iAcceptIter = new int[] { 0, 0, 0 }; 229 iMoves = 0; 230 iAbsValue = 0; 231 } 232 233 /** 234 * Cool temperature 235 * @param solution current solution 236 */ 237 protected void cool(Solution<V, T> solution) { 238 info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)) 239 + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " + 240 "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) + 241 ", Best=" + iDF2.format(solution.getBestValue()) + 242 " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %), " + 243 "Step=" + iDF2.format(1.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iTemperatureLength)); 244 if (iLastImprovingIter < 0 || iIter > iLastImprovingIter + iTemperatureLength) { 245 // Do not update the temperature when solution was improved recently 246 iTemperature *= getCoolingRate(solution.getAssignment().getIndex() - 1); 247 info("Temperature decreased to " + iDF5.format(iTemperature) + " " + "(#moves=" + iMoves + ", rms(value)=" 248 + iDF2.format(Math.sqrt(iAbsValue / iMoves)) + ", " + "accept=-" 249 + iDF2.format(100.0 * iAcceptIter[0] / iMoves) + "/" 250 + iDF2.format(100.0 * iAcceptIter[1] / iMoves) + "/+" 251 + iDF2.format(100.0 * iAcceptIter[2] / iMoves) + "%, " 252 + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + 253 "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " + 254 "p(+1)=" + iDF5.format(100.0 * prob(1)) + "%, " + 255 "p(+10)=" + iDF10.format(100.0 * prob(10)) + "%)"); 256 iAbsValue = 0; 257 iAcceptIter = new int[] { 0, 0, 0 }; 258 iMoves = 0; 259 } 260 logNeibourStatus(); 261 iLastCoolingIter = iIter; 262 } 263 264 /** 265 * Reheat temperature 266 * @param solution current solution 267 */ 268 protected void reheat(Solution<V, T> solution) { 269 iTemperature *= iReheatRate; 270 info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)) 271 + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " + 272 "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) + 273 ", Best=" + iDF2.format(solution.getBestValue()) + 274 " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)"); 275 info("Temperature increased to " + iDF5.format(iTemperature) + " " 276 + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + "p(+0.1)=" 277 + iDF2.format(100.0 * prob(0.1)) + "%, " + "p(+1)=" + iDF5.format(10.0 * prob(1)) + "%, " + "p(+10)=" 278 + iDF10.format(100.0 * prob(10)) + "%)"); 279 logNeibourStatus(); 280 iLastReheatIter = iIter; 281 if (iTemperature > iMaximalTemperature) { 282 // Max temperature reached -> restore best solution and stop re-heating 283 restoreBest(solution); 284 iLastImprovingIter = -1; 285 } 286 iBestValue = solution.getBestValue(); 287 setProgressPhase("Simulated Annealing [" + iDF5.format(iTemperature) + "]..."); 288 } 289 290 /** 291 * restore best ever found solution 292 * @param solution current solution 293 */ 294 protected void restoreBest(Solution<V, T> solution) { 295 info("Best solution restored."); 296 solution.restoreBest(); 297 iLastBestIter = iIter; 298 } 299 300 /** 301 * Neighbour acceptance probability 302 * 303 * @param value 304 * absolute or relative value of the proposed change (neighbour) 305 * @return probability of acceptance of a change (neighbour) 306 */ 307 protected double prob(double value) { 308 return prob(value, iTemperature); 309 } 310 311 /** 312 * Neighbour acceptance probability 313 * 314 * @param value 315 * absolute or relative value of the proposed change (neighbour) 316 * @param temp 317 * current temperature 318 * @return probability of acceptance of a change (neighbour) 319 */ 320 protected double prob(double value, double temp) { 321 if (iStochasticHC) 322 return 1.0 / (1.0 + Math.exp(value / temp)); 323 else 324 return (value <= 0.0 ? 1.0 : Math.exp(-value / temp)); 325 } 326 327 /** 328 * True if the given neighbour is to be be accepted 329 * 330 * @param assignment 331 * current assignment 332 * @param neighbour 333 * proposed move 334 * @return true if generated random number is below the generated probability 335 */ 336 @Override 337 protected boolean accept(Assignment<V, T> assignment, Model<V, T> model, Neighbour<V, T> neighbour, double value, boolean lazy) { 338 iMoves ++; 339 iAbsValue += value * value; 340 double v = (iRelativeAcceptance ? value : (lazy ? model.getTotalValue(assignment) : value + model.getTotalValue(assignment)) - iBestValue); 341 if (iTrainingIterations < iTrainingValues) { 342 if (v <= 0.0) { 343 iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++; 344 return true; 345 } else { 346 iTrainingIterations ++; iTrainingTotal += v; 347 } 348 return false; 349 } 350 double prob = prob(v); 351 if (v > 0) { 352 iTrainingIterations ++; iTrainingTotal += v; 353 } 354 if (prob >= 1.0 || ToolBox.random() < prob) { 355 iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++; 356 return true; 357 } 358 return false; 359 } 360 361 /** 362 * Increment iteration counter, cool/reheat/restoreBest if necessary 363 */ 364 @Override 365 protected void incIteration(Solution<V, T> solution) { 366 super.incIteration(solution); 367 iIter++; 368 if (isMaster(solution)) { 369 if (iInitialTemperature <= 0.0) { 370 if (iTrainingIterations < iTrainingValues) { 371 setProgress(Math.round(100.0 * iTrainingIterations / iTrainingValues)); 372 return; 373 } else { 374 train(solution); 375 } 376 } 377 if (iLastBestIter >= 0 && iIter > iLastBestIter + iRestoreBestLength) 378 restoreBest(solution); 379 if (iTemperature < iMinimalTemperature) { // Minimal temperature reached, start going up again 380 // If the current solution is more than 0.1% worse, restore best 381 if ((solution.getModel().getTotalValue(solution.getAssignment()) - solution.getBestValue()) / Math.abs(solution.getBestValue()) >= 0.001) 382 restoreBest(solution); 383 // Set last improving iteration so that the re-heating would continue 384 iLastImprovingIter = iIter; 385 // Do the first reheat 386 reheat(solution); 387 } else if (iLastImprovingIter >= 0 && iIter > Math.max(iLastReheatIter, iLastImprovingIter) + iReheatLength) { 388 reheat(solution); 389 } else if (iIter > iLastCoolingIter + iTemperatureLength) { 390 cool(solution); 391 } 392 setProgress(Math.round(100.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iReheatLength)); 393 } 394 } 395 396 /** 397 * Memorize the iteration when the last best solution was found. 398 */ 399 @Override 400 public void bestSaved(Solution<V, T> solution) { 401 super.bestSaved(solution); 402 if (iLastImprovingIter < 0 || Math.abs(iBestValue - solution.getBestValue()) / Math.max(Math.abs(iBestValue), Math.abs(solution.getBestValue())) >= 0.0001) { 403 iLastImprovingIter = iIter; 404 iLastBestIter = iIter; 405 iBestValue = solution.getBestValue(); 406 } 407 } 408 } 409 410}