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