David Buchan
05-18-2006, 12:14 PM
Hi guys,
I see some people were asking if anyone has made use of the simulated annealing routine. I can happily report that I've used it to great effect.
I'd like to ask a few questions though regarding the code itself, which, perhaps someone (or the author(s)) could answer, and I've written a brief description of the problem solved so people can see a successful real life application.
Success Story:
I've made use of the simulated annealing technique to solve for target exit irradiations when simulating CANDU reactor cores. A CANDU reactor is composed of a regular square lattice of parallel fuel channels, all inside a large cyclindrical vessel (calandria). Within each fuel channel are 12 or 13 fuel bundles, end-to-end. We fuel on-power by stuffing fresh fuel in one end, and taking the irradiated fuel out the other end (one bundle in - one bundle out). The rate at which we fuel affects the channel power since fresher fuel being pushed in one end and burned-up fuel being removed on the other end results in a positive reactivity insertion.
From a simulation perspective, to determine a time-averaged distribution of fuel bundle irradiation and flux in the core, we simulate long term reactor operation with refuelling. This is accomplished by iteratively solving for bundle irradiation, channel axial flux, core flux distribution, and xenon concentration, while conforming to any bundle-shift rules for fuelling. One parameter which affects the calculated axial distribution of bundle irradiations along a channel is the expected bundle exit irradiation at the time of refuelling. This is called the "target" exit irradiation and it must be specified by the user of the simulation code. It directly affects the frequency at which the channel is visited by the fuelling machine.
We'd like to model this fuelling process such that we fuel in our model at a rate that causes the channel powers in the model to match, as closely as possible, those real channel powers in the field. Obviously, there are other factors affecting channel power, but the fuelling rate is the parameter of interest in this study. This is actually a very frustrating problem because, when one varies target exit irradiations in one location in the core, channel powers pop up or depress elsewhere, much like poking a balloon. At the same time, you're trying to get the core multiplication constant as close to 1 as possible. The goal is to have simulated channel powers match measured historical channel powers, and at the same time, achieve the desired multiplication constant.
I divided the core up into regions, where all channels in a region are assigned the same target exit irradiation. Theoretically, at one extreme, each channel could be defined as an individual region while, at the other extreme, all channels could be included in a single region. I chose some best guesses for each region's target exit irrad and ran a simulation to get channel powers. I then perturbed those target exit irrads, one at a time, and re-ran the simulation to get perturbed channel powers. All of this is automated, so it only takes about 1 minute.
I calculated the average of channel power differences (measured-simulated)/measured for each region. Then I created a coupling coefficient matrix [a] and intercept [b], where I assumed a linear dependence between channel power differences in region i and the target exit irradiation in region j. Thus
[y]=[a]*[x] + [b]
where:
yi is average channel power difference in region i,
aij gives the effect on average channel power difference yi in region i of changing target exit irradiation xj in region j, and
bi is the intercept.
I defined a metric of success based upon these goals as follows. The metric is a single number calculated by summing together the absolute value of average channel power difference from each region, as well as the absolute value of multiplication constant (Keff)-1.0. A smaller value is better ("downhill"), a larger value is worse ("uphill"). This is the "function" value I am trying to minimize with simulated annealing. Obviously for this problem, there is no generic formula to compute the effect of a change in target exit irradiation on channel power at some arbitrary location in the core (other than another full blown simulation), and I don't have a Jacobian or Hessian available. So the simulated annealing method is attractive to me since it only requires function evaluations. Creating the linear matrix just alleviates the need to run the full simulation each time I need to evaluate the "function."
metric = sum-over-all-i(average of channel power differences in region i) + (keff-1.0)
We define an N-dimensional space, where N is the number of variables we must adjust to optimize our solution. In this case, N is the number of target exit irradiation values.
With five regions, simulated annealing solves the problem so fast, it's done by the time I get my finger off the Enter key. The channel power differences are about 1/100th of a percent and Keff is 1 with 9 zeros! A confirmatory simulation is then done to show whether my linear-dependence assumption was correct. Usually it's good enough. If not, I recompute the coupling coefficient matrix, starting my perturbations from the set of exit irrads I just got from simulated annealing, and then repeat the whole process.
This application does not seem to be very sensitive to annealing schedule.
Previously, I was doing this by simply generating random normal target exit irradiation values and selecting the best ones - a sort of pseudo Monte Carlo method. It was very slow, and with larger numbers of regions, it becomes hopelessly slow, since you'd need trillions of Monte Carlo guesses to find the good solutions. I'm not sure if this problem falls under the np-hard class, but it is certainly "hard" to solve by any other means. Simulated annealing makes it a piece of cake.
Questions about code:
1) In their original paper on the downhill simplex method (Nelder and Mead, 1965, Computer Journal), they define the centroid P-bar where i does NOT equal h. h is the subscript for the vertex with the highest (worst) function value. In amebsa.f, i=h IS included when computing the centroid (actually psum, but not divided by ndim till necessary), if I'm reading this code correctly. Why was that choice made? Does it matter?
2) I don't understand the implementation of the reflection formula in amotsa.f, or amotry.f for that matter. If I take Nelder and Mead's formulae, and change their requirements for alpha to be positive for reflections, to alpha being negative, I can use the same routine to do reflections, 1-D expansions, and 1-D contractions. alpha is negative for reflections, >1 for expansions, and between 0 and 1 for contractions. The formula in amotsa.f uses fac, fac1 and fac2. I would've expected fac2 to be just fac instead. Then my argument above makes sense and fac is just alpha, beta, or gamma as in Nelder and Mead, except alpha is negative. Of course, when I modify the code to do that, it gives wrong answers. So why do we need fac2 to be what it is in amotsa.f?
3) (possibly related to question 2) On page 403 of NR in FORTRAN, it reads, "These steps are called reflections, and they are constructed to conserve the volume of the simplex (hence maintain its nondegeneracy)." If I read Nelder and Mead correctly, they do not conserve the volume of the simplex (I think - at least, I don't recall them mentioning it). I have two parts to this question:
a) I guess I need more explanation of what degeneracy we're avoiding. Avoiding simplexes with different volumes that give the same centroid? Is that bad?
b) How is the volume conserved in the code? Is that what's going on in question 2?
Also, a comment: I think (could be wrong) that the book should point out that you must choose an iiter value that is >> than the number of dimensions of the problem. Otherwise, amebsa will return after a single n-dimensional (overall) contraction takes place, and no further reflections or 1-dim expansions and 1-dim contractions can take place prior to the next temperature reduction.
Thanks,
Dave Buchan
Toronto
I see some people were asking if anyone has made use of the simulated annealing routine. I can happily report that I've used it to great effect.
I'd like to ask a few questions though regarding the code itself, which, perhaps someone (or the author(s)) could answer, and I've written a brief description of the problem solved so people can see a successful real life application.
Success Story:
I've made use of the simulated annealing technique to solve for target exit irradiations when simulating CANDU reactor cores. A CANDU reactor is composed of a regular square lattice of parallel fuel channels, all inside a large cyclindrical vessel (calandria). Within each fuel channel are 12 or 13 fuel bundles, end-to-end. We fuel on-power by stuffing fresh fuel in one end, and taking the irradiated fuel out the other end (one bundle in - one bundle out). The rate at which we fuel affects the channel power since fresher fuel being pushed in one end and burned-up fuel being removed on the other end results in a positive reactivity insertion.
From a simulation perspective, to determine a time-averaged distribution of fuel bundle irradiation and flux in the core, we simulate long term reactor operation with refuelling. This is accomplished by iteratively solving for bundle irradiation, channel axial flux, core flux distribution, and xenon concentration, while conforming to any bundle-shift rules for fuelling. One parameter which affects the calculated axial distribution of bundle irradiations along a channel is the expected bundle exit irradiation at the time of refuelling. This is called the "target" exit irradiation and it must be specified by the user of the simulation code. It directly affects the frequency at which the channel is visited by the fuelling machine.
We'd like to model this fuelling process such that we fuel in our model at a rate that causes the channel powers in the model to match, as closely as possible, those real channel powers in the field. Obviously, there are other factors affecting channel power, but the fuelling rate is the parameter of interest in this study. This is actually a very frustrating problem because, when one varies target exit irradiations in one location in the core, channel powers pop up or depress elsewhere, much like poking a balloon. At the same time, you're trying to get the core multiplication constant as close to 1 as possible. The goal is to have simulated channel powers match measured historical channel powers, and at the same time, achieve the desired multiplication constant.
I divided the core up into regions, where all channels in a region are assigned the same target exit irradiation. Theoretically, at one extreme, each channel could be defined as an individual region while, at the other extreme, all channels could be included in a single region. I chose some best guesses for each region's target exit irrad and ran a simulation to get channel powers. I then perturbed those target exit irrads, one at a time, and re-ran the simulation to get perturbed channel powers. All of this is automated, so it only takes about 1 minute.
I calculated the average of channel power differences (measured-simulated)/measured for each region. Then I created a coupling coefficient matrix [a] and intercept [b], where I assumed a linear dependence between channel power differences in region i and the target exit irradiation in region j. Thus
[y]=[a]*[x] + [b]
where:
yi is average channel power difference in region i,
aij gives the effect on average channel power difference yi in region i of changing target exit irradiation xj in region j, and
bi is the intercept.
I defined a metric of success based upon these goals as follows. The metric is a single number calculated by summing together the absolute value of average channel power difference from each region, as well as the absolute value of multiplication constant (Keff)-1.0. A smaller value is better ("downhill"), a larger value is worse ("uphill"). This is the "function" value I am trying to minimize with simulated annealing. Obviously for this problem, there is no generic formula to compute the effect of a change in target exit irradiation on channel power at some arbitrary location in the core (other than another full blown simulation), and I don't have a Jacobian or Hessian available. So the simulated annealing method is attractive to me since it only requires function evaluations. Creating the linear matrix just alleviates the need to run the full simulation each time I need to evaluate the "function."
metric = sum-over-all-i(average of channel power differences in region i) + (keff-1.0)
We define an N-dimensional space, where N is the number of variables we must adjust to optimize our solution. In this case, N is the number of target exit irradiation values.
With five regions, simulated annealing solves the problem so fast, it's done by the time I get my finger off the Enter key. The channel power differences are about 1/100th of a percent and Keff is 1 with 9 zeros! A confirmatory simulation is then done to show whether my linear-dependence assumption was correct. Usually it's good enough. If not, I recompute the coupling coefficient matrix, starting my perturbations from the set of exit irrads I just got from simulated annealing, and then repeat the whole process.
This application does not seem to be very sensitive to annealing schedule.
Previously, I was doing this by simply generating random normal target exit irradiation values and selecting the best ones - a sort of pseudo Monte Carlo method. It was very slow, and with larger numbers of regions, it becomes hopelessly slow, since you'd need trillions of Monte Carlo guesses to find the good solutions. I'm not sure if this problem falls under the np-hard class, but it is certainly "hard" to solve by any other means. Simulated annealing makes it a piece of cake.
Questions about code:
1) In their original paper on the downhill simplex method (Nelder and Mead, 1965, Computer Journal), they define the centroid P-bar where i does NOT equal h. h is the subscript for the vertex with the highest (worst) function value. In amebsa.f, i=h IS included when computing the centroid (actually psum, but not divided by ndim till necessary), if I'm reading this code correctly. Why was that choice made? Does it matter?
2) I don't understand the implementation of the reflection formula in amotsa.f, or amotry.f for that matter. If I take Nelder and Mead's formulae, and change their requirements for alpha to be positive for reflections, to alpha being negative, I can use the same routine to do reflections, 1-D expansions, and 1-D contractions. alpha is negative for reflections, >1 for expansions, and between 0 and 1 for contractions. The formula in amotsa.f uses fac, fac1 and fac2. I would've expected fac2 to be just fac instead. Then my argument above makes sense and fac is just alpha, beta, or gamma as in Nelder and Mead, except alpha is negative. Of course, when I modify the code to do that, it gives wrong answers. So why do we need fac2 to be what it is in amotsa.f?
3) (possibly related to question 2) On page 403 of NR in FORTRAN, it reads, "These steps are called reflections, and they are constructed to conserve the volume of the simplex (hence maintain its nondegeneracy)." If I read Nelder and Mead correctly, they do not conserve the volume of the simplex (I think - at least, I don't recall them mentioning it). I have two parts to this question:
a) I guess I need more explanation of what degeneracy we're avoiding. Avoiding simplexes with different volumes that give the same centroid? Is that bad?
b) How is the volume conserved in the code? Is that what's going on in question 2?
Also, a comment: I think (could be wrong) that the book should point out that you must choose an iiter value that is >> than the number of dimensions of the problem. Otherwise, amebsa will return after a single n-dimensional (overall) contraction takes place, and no further reflections or 1-dim expansions and 1-dim contractions can take place prior to the next temperature reduction.
Thanks,
Dave Buchan
Toronto