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.5; 085 private double iCoolingRate = 0.95; 086 private double iReheatRate = -1; 087 private long iTemperatureLength = 25000; 088 private double iReheatLengthCoef = 5.0; 089 private double iRestoreBestLengthCoef = -1; 090 private boolean iStochasticHC = false; 091 private boolean iRelativeAcceptance = true; 092 private Double[] iCoolingRateAdjusts = null; 093 private int iTrainingValues = 10000; 094 private double iTrainingProbability = 0.01; 095 096 /** 097 * Constructor. Following problem properties are considered: 098 * <ul> 099 * <li>SimulatedAnnealing.InitialTemperature ... initial temperature (default 1.5) 100 * <li>SimulatedAnnealing.TemperatureLength ... temperature length (number of iterations between temperature decrements, default 25000) 101 * <li>SimulatedAnnealing.CoolingRate ... temperature cooling rate (default 0.95) 102 * <li>SimulatedAnnealing.ReheatLengthCoef ... temperature re-heat length coefficient (multiple of temperature length, default 5) 103 * <li>SimulatedAnnealing.ReheatRate ... temperature re-heating rate (default (1/coolingRate)^(reheatLengthCoef*1.7)) 104 * <li>SimulatedAnnealing.RestoreBestLengthCoef ... restore best length coefficient (multiple of re-heat length, default reheatLengthCoef^2) 105 * <li>SimulatedAnnealing.StochasticHC ... true for stochastic search acceptance criterion, false for simulated annealing acceptance (default false) 106 * <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) 107 * <li>SimulatedAnnealing.Neighbours ... semicolon separated list of classes implementing {@link NeighbourSelection} 108 * <li>SimulatedAnnealing.AdditionalNeighbours ... semicolon separated list of classes implementing {@link NeighbourSelection} 109 * <li>SimulatedAnnealing.Random ... when true, a neighbour selector is selected randomly 110 * <li>SimulatedAnnealing.Update ... when true, a neighbour selector is selected using {@link NeighbourSelector#getPoints()} weights (roulette wheel selection) 111 * </ul> 112 * 113 * @param properties 114 * problem properties 115 */ 116 public SimulatedAnnealing(DataProperties properties) { 117 super(properties); 118 iInitialTemperature = properties.getPropertyDouble(getParameterBaseName() + ".InitialTemperature", iInitialTemperature); 119 iReheatRate = properties.getPropertyDouble(getParameterBaseName() + ".ReheatRate", iReheatRate); 120 iCoolingRate = properties.getPropertyDouble(getParameterBaseName() + ".CoolingRate", iCoolingRate); 121 iRelativeAcceptance = properties.getPropertyBoolean(getParameterBaseName() + ".RelativeAcceptance", iRelativeAcceptance); 122 iStochasticHC = properties.getPropertyBoolean(getParameterBaseName() + ".StochasticHC", iStochasticHC); 123 iTemperatureLength = properties.getPropertyLong(getParameterBaseName() + ".TemperatureLength", iTemperatureLength); 124 iReheatLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".ReheatLengthCoef", iReheatLengthCoef); 125 iRestoreBestLengthCoef = properties.getPropertyDouble(getParameterBaseName() + ".RestoreBestLengthCoef", iRestoreBestLengthCoef); 126 iCoolingRateAdjusts = properties.getPropertyDoubleArry(getParameterBaseName() + ".CoolingRateAdjustments", null); 127 iTrainingValues = properties.getPropertyInt(getParameterBaseName() + ".TrainingValues", iTrainingValues); 128 iTrainingProbability = properties.getPropertyDouble(getParameterBaseName() + ".TrainingProbability", iTrainingProbability); 129 if (iReheatRate < 0) 130 iReheatRate = Math.pow(1 / iCoolingRate, iReheatLengthCoef * 1.7); 131 if (iRestoreBestLengthCoef < 0) 132 iRestoreBestLengthCoef = iReheatLengthCoef * iReheatLengthCoef; 133 } 134 135 @Override 136 public String getParameterBaseName() { return "SimulatedAnnealing"; } 137 138 @Override 139 public NeighbourSearchContext createAssignmentContext(Assignment<V, T> assignment) { 140 return new SimulatedAnnealingContext(); 141 } 142 143 public class SimulatedAnnealingContext extends NeighbourSearchContext { 144 private double iTemperature = 0.0; 145 private int iMoves = 0; 146 private double iAbsValue = 0; 147 private double iBestValue = 0; 148 private long iLastImprovingIter = -1; 149 private long iLastReheatIter = 0; 150 private long iLastCoolingIter = 0; 151 private int iAcceptIter[] = new int[] { 0, 0, 0 }; 152 private long iReheatLength = 0; 153 private long iRestoreBestLength = 0; 154 private int iTrainingIterations = 0; 155 private double iTrainingTotal = 0.0; 156 157 /** Setup the temperature */ 158 @Override 159 protected void activate(Solution<V, T> solution) { 160 super.activate(solution); 161 iTrainingTotal = 0.0; iTrainingIterations = 0; 162 iTemperature = iInitialTemperature; 163 iReheatLength = Math.round(iReheatLengthCoef * iTemperatureLength); 164 iRestoreBestLength = Math.round(iRestoreBestLengthCoef * iTemperatureLength); 165 iLastImprovingIter = -1; 166 } 167 168 protected double getCoolingRate(int idx) { 169 if (idx < 0 || iCoolingRateAdjusts == null || idx >= iCoolingRateAdjusts.length || iCoolingRateAdjusts[idx] == null) return iCoolingRate; 170 return iCoolingRate * iCoolingRateAdjusts[idx]; 171 } 172 173 /** 174 * Set initial temperature based on the training period 175 * @param solution current solution 176 */ 177 protected void train(Solution<V, T> solution) { 178 double value = iTrainingTotal / iTrainingIterations; 179 if (iStochasticHC) { 180 iInitialTemperature = value / Math.log(1.0 / iTrainingProbability - 1.0); 181 } else { 182 iInitialTemperature = - value / Math.log(iTrainingProbability); 183 } 184 iTemperature = iInitialTemperature; 185 info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)) 186 + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " + 187 "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) + 188 ", Best=" + iDF2.format(solution.getBestValue()) + 189 " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)"); 190 info("Temperature set to " + iDF5.format(iTemperature) + " " + "(" + 191 "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " + 192 "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " + 193 "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " + 194 "p(+" + iDF2.format(value) + ")=" + iDF5.format(100.0 * prob(value)) + "%)"); 195 logNeibourStatus(); 196 iIter = 0; iLastImprovingIter = -1; iBestValue = solution.getBestValue(); 197 iAcceptIter = new int[] { 0, 0, 0 }; 198 iMoves = 0; 199 iAbsValue = 0; 200 } 201 202 /** 203 * Cool temperature 204 * @param solution current solution 205 */ 206 protected void cool(Solution<V, T> solution) { 207 iTemperature *= getCoolingRate(solution.getAssignment().getIndex() - 1); 208 info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)) 209 + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " + 210 "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) + 211 ", Best=" + iDF2.format(solution.getBestValue()) + 212 " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %), " + 213 "Step=" + iDF2.format(1.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iTemperatureLength)); 214 info("Temperature decreased to " + iDF5.format(iTemperature) + " " + "(#moves=" + iMoves + ", rms(value)=" 215 + iDF2.format(Math.sqrt(iAbsValue / iMoves)) + ", " + "accept=-" 216 + iDF2.format(100.0 * iAcceptIter[0] / iMoves) + "/" 217 + iDF2.format(100.0 * iAcceptIter[1] / iMoves) + "/+" 218 + iDF2.format(100.0 * iAcceptIter[2] / iMoves) + "%, " 219 + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + 220 "p(+0.1)=" + iDF2.format(100.0 * prob(0.1)) + "%, " + 221 "p(+1)=" + iDF2.format(100.0 * prob(1)) + "%, " + 222 "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%)"); 223 logNeibourStatus(); 224 iLastCoolingIter = iIter; 225 iAcceptIter = new int[] { 0, 0, 0 }; 226 iMoves = 0; 227 iAbsValue = 0; 228 } 229 230 /** 231 * Reheat temperature 232 * @param solution current solution 233 */ 234 protected void reheat(Solution<V, T> solution) { 235 iTemperature *= iReheatRate; 236 info("Iter=" + iIter / 1000 + (iLastImprovingIter < 0 ? "" : "k, NonImpIter=" + iDF2.format((iIter - iLastImprovingIter) / 1000.0)) 237 + "k, Speed=" + iDF2.format(1000.0 * iIter / (JProf.currentTimeMillis() - iT0)) + " it/s, " + 238 "Value=" + iDF2.format(solution.getModel().getTotalValue(solution.getAssignment())) + 239 ", Best=" + iDF2.format(solution.getBestValue()) + 240 " (" + iDF2.format(100.0 * solution.getModel().getTotalValue(solution.getAssignment()) / solution.getBestValue()) + " %)"); 241 info("Temperature increased to " + iDF5.format(iTemperature) + " " 242 + (prob(-1) < 1.0 ? "p(-1)=" + iDF2.format(100.0 * prob(-1)) + "%, " : "") + "p(+1)=" 243 + iDF2.format(100.0 * prob(1)) + "%, " + "p(+10)=" + iDF5.format(100.0 * prob(10)) + "%, " + "p(+100)=" 244 + iDF10.format(100.0 * prob(100)) + "%)"); 245 logNeibourStatus(); 246 iLastReheatIter = iIter; 247 if (iTemperature > iInitialTemperature) 248 iLastImprovingIter = -1; 249 iBestValue = solution.getBestValue(); 250 setProgressPhase("Simulated Annealing [" + iDF5.format(iTemperature) + "]..."); 251 } 252 253 /** 254 * restore best ever found solution 255 * @param solution current solution 256 */ 257 protected void restoreBest(Solution<V, T> solution) { 258 info("Best solution restored."); 259 solution.restoreBest(); 260 iLastImprovingIter = -1; 261 } 262 263 /** 264 * Neighbour acceptance probability 265 * 266 * @param value 267 * absolute or relative value of the proposed change (neighbour) 268 * @return probability of acceptance of a change (neighbour) 269 */ 270 protected double prob(double value) { 271 if (iStochasticHC) 272 return 1.0 / (1.0 + Math.exp(value / iTemperature)); 273 else 274 return (value <= 0.0 ? 1.0 : Math.exp(-value / iTemperature)); 275 } 276 277 /** 278 * True if the given neighbour is to be be accepted 279 * 280 * @param assignment 281 * current assignment 282 * @param neighbour 283 * proposed move 284 * @return true if generated random number is below the generated probability 285 */ 286 @Override 287 protected boolean accept(Assignment<V, T> assignment, Model<V, T> model, Neighbour<V, T> neighbour, double value, boolean lazy) { 288 iMoves ++; 289 iAbsValue += value * value; 290 double v = (iRelativeAcceptance ? value : (lazy ? model.getTotalValue(assignment) : value + model.getTotalValue(assignment)) - iBestValue); 291 if (iTrainingIterations < iTrainingValues) { 292 if (v <= 0.0) { 293 iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++; 294 return true; 295 } else { 296 iTrainingIterations ++; iTrainingTotal += v; 297 } 298 return false; 299 } 300 double prob = prob(v); 301 if (v > 0) { 302 iTrainingIterations ++; iTrainingTotal += v; 303 } 304 if (prob >= 1.0 || ToolBox.random() < prob) { 305 iAcceptIter[value < 0.0 ? 0 : value > 0.0 ? 2 : 1]++; 306 return true; 307 } 308 return false; 309 } 310 311 /** 312 * Increment iteration counter, cool/reheat/restoreBest if necessary 313 */ 314 @Override 315 protected void incIteration(Solution<V, T> solution) { 316 super.incIteration(solution); 317 iIter++; 318 if (isMaster(solution)) { 319 if (iInitialTemperature <= 0.0) { 320 if (iTrainingIterations < iTrainingValues) { 321 setProgress(Math.round(100.0 * iTrainingIterations / iTrainingValues)); 322 return; 323 } else { 324 train(solution); 325 } 326 } 327 if (iLastImprovingIter >= 0 && iIter > iLastImprovingIter + iRestoreBestLength) 328 restoreBest(solution); 329 if (iLastImprovingIter >= 0 && iIter > Math.max(iLastReheatIter, iLastImprovingIter) + iReheatLength) 330 reheat(solution); 331 if (iIter > iLastCoolingIter + iTemperatureLength) 332 cool(solution); 333 setProgress(Math.round(100.0 * (iIter - Math.max(iLastReheatIter, iLastImprovingIter)) / iReheatLength)); 334 } 335 } 336 337 /** 338 * Memorize the iteration when the last best solution was found. 339 */ 340 @Override 341 public void bestSaved(Solution<V, T> solution) { 342 super.bestSaved(solution); 343 if (Math.abs(iBestValue - solution.getBestValue()) >= 1.0) { 344 iLastImprovingIter = iIter; 345 iBestValue = solution.getBestValue(); 346 } 347 } 348 } 349 350}