paulromano
commented
6 years ago I'm no expert in optimization. If you develop something that you think would be generally useful to others, we would welcome a pull request though.
Open SterlingButters opened 6 years ago
paulromano
commented
6 years ago I'm no expert in optimization. If you develop something that you think would be generally useful to others, we would welcome a pull request though.
SterlingButters
commented
5 years ago So my previous, goal was a bit unrealistic (ok, very) - I was trying to see if an optimization algorithm could predict optimal fuel loading maps for the whole core. 1) I should have started out with something a little more rudimentary and 2) I shouldn't have even thought about being able to obtain results in a non-HPC environment. However, I do have some logic to produce that I would like to share. I'm no expert in optimization myself but I think the combination of the following "tools" would make for some interesting results:
The only thing that has me concerned is that the "objectives" referenced in the corresponding links are functions of the same variable x while in the case below, the compromise_function is multivariate however, I can't see why, mathematically, this wouldn't work.
# Create the model. `ppm_Boron` will be the parametric variable.
def build_model(ppm_Boron, enrichment, fuel_Radius):
# Create the pin materials
fuel = openmc.Material(name='1.6% Fuel')
fuel.set_density('g/cm3', 10.31341)
fuel.add_element('U', 1., enrichment=enrichment)
fuel.add_element('O', 2.)
zircaloy = openmc.Material(name='Zircaloy')
zircaloy.set_density('g/cm3', 6.55)
zircaloy.add_element('Zr', 1.)
water = openmc.Material(name='Borated Water')
water.set_density('g/cm3', 0.741)
water.add_element('H', 2.)
water.add_element('O', 1.)
# Include the amount of boron in the water based on the ppm,
# neglecting the other constituents of boric acid
water.add_element('B', ppm_Boron * 1E-6)
# Instantiate a Materials object
materials = openmc.Materials([fuel, zircaloy, water])
# Create cylinders for the fuel and clad
fuel_outer_radius = openmc.ZCylinder(R=fuel_Radius)
clad_outer_radius = openmc.ZCylinder(R=fuel_Radius + 0.45720 - 0.39218)
# Create boundary planes to surround the geometry
min_x = openmc.XPlane(x0=-0.63, boundary_type='reflective')
max_x = openmc.XPlane(x0=+0.63, boundary_type='reflective')
min_y = openmc.YPlane(y0=-0.63, boundary_type='reflective')
max_y = openmc.YPlane(y0=+0.63, boundary_type='reflective')
# Create fuel Cell
fuel_cell = openmc.Cell(name='{}% Fuel'.format(enrichment))
fuel_cell.fill = fuel
fuel_cell.region = -fuel_outer_radius
# Create a clad Cell
clad_cell = openmc.Cell(name='{}% Clad'.format(enrichment))
clad_cell.fill = zircaloy
clad_cell.region = +fuel_outer_radius & -clad_outer_radius
# Create a moderator Cell
moderator_cell = openmc.Cell(name='{}% Moderator'.format(enrichment))
moderator_cell.fill = water
moderator_cell.region = +clad_outer_radius & (+min_x & -max_x & +min_y & -max_y)
# Create root Universe
root_universe = openmc.Universe(name='root universe', universe_id=0)
root_universe.add_cells([fuel_cell, clad_cell, moderator_cell])
# Create Geometry and set root universe
geometry = openmc.Geometry(root_universe)
# Finish with the settings file
settings = openmc.Settings()
settings.batches = 300
settings.inactive = 20
settings.particles = 1000
settings.run_mode = 'eigenvalue'
# Create an initial uniform spatial source distribution over fissionable zones
bounds = [-0.63, -0.63, -10, 0.63, 0.63, 10.]
uniform_dist = openmc.stats.Box(bounds[:3], bounds[3:], only_fissionable=True)
settings.source = openmc.source.Source(space=uniform_dist)
# We dont need a tallies file so dont waste the disk input/output time
settings.output = {'tallies': False}
model = openmc.model.Model(geometry, materials, settings)
return model
def compromise_function(ppm, enrich, radius):
w_a = .95
w_b = (1-w_a)/3
w_c = (1-w_a)/3
w_d = (1-w_a)/3
p = 1
k_eff = build_model(ppm, enrich, radius).run()[0]
A = (w_a * ( (1 - k_eff) / (1 - .95) ))**p
B = (w_b * ( (0 - ppm) / (0 - 2500) ))**p
C = (w_c * ((0 - enrich) / (0 - 5) ))**p
D = (w_d * ((0 - radius) / (0 - 0.39218)))**p
L = (A + B + C + D)**(1/p)
return L
theta0 = {"ppm": 500, "enrich": 2.0, "radius": 0.45}
m = spsa.SPSA_minimization(compromise_function, theta0, 150)
minimum = m.run()
print("minimums =", minimum)
print("goal at minimum =", m.evaluate_goal(minimum))Contents of spsa (got from someone else's github repo so I will need to sort through logic later):
import utils
class SPSA_minimization:
def __init__(self, f, theta0, max_iter, constraints = None, options = {}):
"""
The constructor of a SPSA_minimization object.
We use the notations and ideas of the following articles:
• Spall JC (1998), Implementation of the Simultuaneous Perturbation
Algorithm for Stochastic Optimization, IEEE Trans Aerosp Electron
Syst 34(3):817–823
• Kocsis & Szepesvari (2006), Universal Parameter Optimisation in
Games based on SPSA, Mach Learn 63:249–286
Args:
f (function) :
The function to minimize.
theta0 (dict) :
The starting point of the minimization.
max_iter (int) :
The number of iterations of the algorithm.
constraints (function, optional) :
A function which maps the current point to the closest point
of the search domain.
options (dict, optional) :
Optional settings of the SPSA algorithm parameters. Default
values taken from the reference articles are used if not
present in options.
"""
# Store the arguments
self.f = f
self.theta0 = theta0
self.iter = 0;
self.max_iter = max_iter
self.constraints = constraints
self.options = options
# some attributes to provide an history of evaluations
self.previous_gradient = {}
self.rprop_previous_g = {}
self.rprop_previous_delta = {}
self.history_eval = array.array('d', range(1000))
self.history_theta = [theta0 for k in range(1000)]
self.history_count = 0
self.best_eval = array.array('d', range(1000))
self.best_theta = [theta0 for k in range(1000)]
self.best_count = 0
# These constants are used throughout the SPSA algorithm
self.a = options.get("a", 1.1)
self.c = options.get("c", 0.1)
self.alpha = options.get("alpha", 0.601) # theoretical alpha=0.601, must be <= 1
self.gamma = options.get("gamma", 0.101) # theoretical gamma=0.101, must be <= 1/6
self.A = options.get("A", max_iter / 10.0)
def run(self):
"""
Return a point which is (hopefully) a minimizer of the goal
function f, starting from point theta0.
Returns:
The point (as a dict) which is (hopefully) a minimizer of "f".
"""
k = 0
theta = self.theta0
while True:
k = k + 1
self.iter = k
if self.constraints is not None:
theta = self.constraints(theta)
print("theta = " + utils.pretty(theta))
c_k = self.c / (k ** self.gamma)
a_k = self.a / ((k + self.A) ** self.alpha)
gradient = self.approximate_gradient(theta, c_k)
#print(str(k) + " gradient = " + utils.pretty(gradient))
# if k % 1000 == 0:
# print(k + utils.pretty(theta) + "norm2(g) = " + str(utils.norm2(gradient)))
# print(k + " theta = " + utils.pretty(theta))
## For SPSA we update with a small step (theta = theta - a_k * gradient)
## theta = utils.linear_combinaison(1.0, theta, -a_k, gradient)
## For steepest descent we update via a constant small step in the gradient direction
mu = -0.01 / max(1.0, utils.norm2(gradient))
theta = utils.linear_combinaison(1.0, theta, mu, gradient)
## For RPROP, we update with information about the sign of the gradients
theta = utils.linear_combinaison(1.0, theta, -0.01, self.rprop(theta, gradient))
## We then move to the point which gives the best average of goal
(avg_goal , avg_theta) = self.average_best_evals(30)
theta = utils.linear_combinaison(0.98, theta, 0.02, avg_theta)
if (k % 100 == 0) or (k <= 1000) :
(avg_goal , avg_theta) = self.average_evaluations(30)
print("iter = " + str(k))
print("mean goal (all) = " + str(avg_goal))
print("mean theta (all) = " + utils.pretty(avg_theta))
(avg_goal , avg_theta) = self.average_best_evals(30)
print("mean goal (best) = " + str(avg_goal))
print("mean theta (best) = " + utils.pretty(avg_theta))
print("-----------------------------------------------------------------")
if k >= self.max_iter:
break
return theta
def evaluate_goal(self, theta):
"""
Return the evaluation of the goal function f at point theta.
We also store an history the 1000 last evaluations, so as to be able
to quickly calculate an average of these last evaluations of the goal
via the helper average_evaluations() : this is handy to monitor the
progress of our minimization algorithm.
"""
v = self.f(**theta)
# store the value in history
self.history_eval [self.history_count % 1000] = v
self.history_theta[self.history_count % 1000] = theta
self.history_count += 1
return v
def approximate_gradient(self, theta, c):
"""
Return an approximation of the gradient of f at point theta.
On repeated calls, the esperance of the series of returned values
converges almost surely to the true gradient of f at theta.
"""
if self.history_count > 0:
current_goal, _ = self.average_evaluations(30)
else:
current_goal = 100000000000000000.0
bernouilli = self.create_bernouilli(theta)
count = 0
while True:
# Calculate two evaluations of f at points M + c * bernouilli and
# M - c * bernouilli to estimate the gradient. We do not want to
# use a null gradient, so we loop until the two functions evaluations
# are different. Another trick is that we use the same seed for the
# random generator for the two function evaluations, to reduce the
# variance of the gradient if the evaluations use simulations (like
# in games).
state = random.getstate()
theta1 = utils.linear_combinaison(1.0, theta, c, bernouilli)
f1 = self.evaluate_goal(theta1)
random.setstate(state)
theta2 = utils.linear_combinaison(1.0, theta, -c, bernouilli)
f2 = self.evaluate_goal(theta2)
if f1 != f2:
break
count = count + 1
if count >= 100:
# print("too many evaluation to find a gradient, function seems flat")
break
# Update the gradient
gradient = {}
for (name, value) in theta.items():
gradient[name] = (f1 - f2) / (2.0 * c * bernouilli[name])
if (f1 > current_goal) and (f2 > current_goal):
print("function seems not decreasing")
gradient = utils.linear_combinaison(0.1, gradient)
# For the correction factor used in the running average for the gradient,
# see the paper "Adam: A Method For Stochastic Optimization, Kingma and Lei Ba"
beta = 0.9
correction = 1.0 / 1.0 - pow(beta, self.iter)
gradient = utils.linear_combinaison((1 - beta), gradient, beta, self.previous_gradient)
gradient = utils.linear_combinaison(correction, gradient)
# Store the current gradient for the next time, to calculate the running average
self.previous_gradient = gradient
# Store the best the two evals f1 and f2 (or both)
if (f1 <= current_goal):
self.best_eval [self.best_count % 1000] = f1
self.best_theta[self.best_count % 1000] = theta1
self.best_count += 1
if (f2 <= current_goal):
self.best_eval [self.best_count % 1000] = f2
self.best_theta[self.best_count % 1000] = theta2
self.best_count += 1
# Return the estimation of the new gradient
return gradient
def create_bernouilli(self, m):
"""
Create a random direction to estimate the stochastic gradient.
We use a Bernouilli distribution : bernouilli = (+1,+1,-1,+1,-1,.....)
"""
bernouilli = {}
for (name, value) in m.items():
bernouilli[name] = 1 if random.randint(0, 1) else -1
g = utils.norm2(self.previous_gradient)
d = utils.norm2(bernouilli)
if g > 0.00001:
bernouilli = utils.linear_combinaison(0.55 , bernouilli, \
0.25 * d / g, self.previous_gradient)
for (name, value) in m.items():
if bernouilli[name] == 0.0:
bernouilli[name] = 0.2
if abs(bernouilli[name]) < 0.2:
bernouilli[name] = 0.2 * utils.sign_of(bernouilli[name])
return bernouilli
def average_evaluations(self, n):
"""
Return the average of the n last evaluations of the goal function.
This is a fast function which uses the last evaluations already
done by the SPSA algorithm to return an approximation of the current
goal value (note that we do not call the goal function another time,
so the returned value is an upper bound of the true value).
"""
assert(self.history_count > 0) , "not enough evaluations in average_evaluations!"
if n <= 0 : n = 1
if n > 1000 : n = 1000
if n > self.history_count : n = self.history_count
sum_eval = 0.0
sum_theta = utils.linear_combinaison(0.0, self.theta0)
for i in range(n):
j = ((self.history_count - 1) % 1000) - i
if j < 0 : j += 1000
if j >= 1000 : j -= 1000
sum_eval += self.history_eval[j]
sum_theta = utils.sum(sum_theta, self.history_theta[j])
# return the average
alpha = 1.0 / (1.0 * n)
return (alpha * sum_eval , utils.linear_combinaison(alpha, sum_theta))
def average_best_evals(self, n):
"""
Return the average of the n last best evaluations of the goal function.
This is a fast function which uses the last evaluations already
done by the SPSA algorithm to return an approximation of the current
goal value (note that we do not call the goal function another time,
so the returned value is an upper bound of the true value).
"""
assert(self.best_count > 0) , "not enough evaluations in average_evaluations!"
if n <= 0 : n = 1
if n > 1000 : n = 1000
if n > self.best_count : n = self.best_count
sum_eval = 0.0
sum_theta = utils.linear_combinaison(0.0, self.theta0)
for i in range(n):
j = ((self.best_count - 1) % 1000) - i
if j < 0 : j += 1000
if j >= 1000 : j -= 1000
sum_eval += self.best_eval[j]
sum_theta = utils.sum(sum_theta, self.best_theta[j])
# return the average
alpha = 1.0 / (1.0 * n)
return (alpha * sum_eval , utils.linear_combinaison(alpha, sum_theta))
def rprop(self, theta, gradient):
# get the previous g of the RPROP algorithm
if self.rprop_previous_g != {}:
previous_g = self.rprop_previous_g
else:
previous_g = gradient
# get the previous delta of the RPROP algorithm
if self.rprop_previous_delta != {}:
delta = self.rprop_previous_delta
else:
delta = gradient
delta = utils.copy_and_fill(delta, 0.5)
p = utils.hadamard_product(previous_g, gradient)
print("gradient = " + utils.pretty(gradient))
print("old_g = " + utils.pretty(previous_g))
print("p = " + utils.pretty(p))
g = {}
eta = {}
for (name, value) in p.items():
if p[name] > 0 : eta[name] = 1.1 ## building speed
if p[name] < 0 : eta[name] = 0.5 ## we have passed a local minima: slow down
if p[name] == 0 : eta[name] = 1.0
delta[name] = eta[name] * delta[name]
delta[name] = min(50.0, delta[name])
delta[name] = max(0.000001, delta[name])
g[name] = gradient[name]
print("g = " + utils.pretty(g))
print("eta = " + utils.pretty(eta))
print("delta = " + utils.pretty(delta))
# store the current g and delta for the next call of the RPROP algorithm
self.rprop_previous_g = g
self.rprop_previous_delta = delta
# calculate the update for the current RPROP
s = utils.hadamard_product(delta, utils.sign(g))
print("sign(g) = " + utils.pretty(utils.sign(g)))
print("s = " + utils.pretty(s))
return sContents of utils:
### Helper functions
def norm2(m):
"""
Return the L2-norm of the point m
"""
s = 0.0
for (name, value) in m.items():
s += value ** 2
return math.sqrt(s)
def norm1(m):
"""
Return the L1-norm of the point m
"""
s = 0.0
for (name, value) in m.items():
s += abs(value)
return s
def linear_combinaison(alpha = 1.0, m1 = {},
beta = None, m2 = {}):
"""
Return the linear combinaison m = alpha * m1 + beta * m2.
"""
if m2 == {}:
m2 = m1
beta = 0.0
m = {}
for (name, value) in m1.items():
m[name] = alpha * value + beta * m2.get(name, 0.0)
return m
def difference(m1, m2):
"""
Return the difference m = m1 - m2.
"""
return linear_combinaison(1.0, m1, -1.0, m2)
def sum(m1, m2):
"""
Return the sum m = m1 + m2.
"""
return linear_combinaison(1.0, m1, 1.0, m2)
def hadamard_product(m1, m2):
"""
Return a vector where each entry is the product of the corresponding
entries in m1 and m2
"""
m = {}
for (name, value) in m1.items():
m[name] = value * m2.get(name, 0.0)
return m
def regulizer(m, lambd , alpha):
"""
Return a regulizer r = lamdb * ((1 - alpha) * norm1(m) + alpha * norm2(m))
Useful to transform non-convex optimization problems into pseudo-convex ones
"""
return lambd * ((1 - alpha) * norm1(m) + alpha * norm2(m))
def sign_of(x):
"""
The sign function for a real or integer parameter
Return -1, 0 or 1 depending of the sign of x
"""
return (x and (1, -1)[x < 0])
def sign(m):
"""
Return a vector filled by the signs of m, component by component
"""
result = {}
for (name, value) in m.items():
result[name] = sign_of(value)
return result
def sqrt(m):
"""
Return a vector filled by the square roots of m, component by component
"""
result = {}
for (name, value) in m.items():
result[name] = math.sqrt(value)
return result
def copy_and_fill(m, v):
"""
Return a vector with the same component as m, filled with v
"""
result = {}
for (name, foo) in m.items():
result[name] = v
return result
def pretty(m):
"""
Return a string representation of m
"""
s = "{ "
for name in sorted(m):
s += name + " ="
if m[name] >= 0: s += " "
s += str("%.4f" % m[name]) + "; "
s += "}"
return sIn terms of results, the outcome I had hoped for did not emerge. Currently,
k-effective (Collision) = 1.66724 +/- 0.00219
k-effective (Track-length) = 1.66747 +/- 0.00347
k-effective (Absorption) = 1.66739 +/- 0.00146
Combined k-effective = 1.66737 +/- 0.00143
Leakage Fraction = 0.00000 +/- 0.00000
...
iter = 125
mean goal (all) = -12.471952006214359
mean theta (all) = { enrich = 18.7208; ppm = 499.9879; radius = 0.2585; }
mean goal (best) = -12.45793831171946
mean theta (best) = { enrich = 18.5106; ppm = 499.9865; radius = 0.2631; }
-----------------------------------------------------------------
theta = { enrich = 19.0559; ppm = 500.0112; radius = 0.2527; }So for whatever reason, the enrichment and fuel radius variables are changing quite dramatically yet the "ideal" value I specified for k_eff was 1. Perhaps I need to take a closer look at the contraint provision in the spsa minimization, specify a higher weight for that objective or some other alternative. All thoughts are welcome!
SterlingButters
commented
5 years ago hmmm I dont think the minimum should be negative - rerunning with revised logic now, will post back with new results.
EDIT: This time the fuel radius was dramatically increased by changing
A = (w_a * ( (1 - k_eff) / (1 - .95) ))**p
B = (w_b * ( (0 - ppm) / (0 - 2500) ))**p
C = (w_c * ((0 - enrich) / (0 - 5) ))**p
D = (w_d * ((0 - radius) / (0 - 0.39218)))**pto
A = (w_a * ( (k_eff - 1) / (1 - .95) ))**p
B = (w_b * ( (ppm - 0) / (0 - 2500) ))**p
C = (w_c * ((enrich - 0) / (0 - 5) ))**p
D = (w_d * ((radius - 0) / (0 - 0.39218)))**pto get positive "distance" of compromise function.
As soon as the fuel radius get sufficiently large:
theta = { enrich = 1.9968; ppm = 999.9771; radius = 0.5340; }
the lattice become insensible and OpenMC can no longer track the particles throwing an error.
Confirmed, I need to look 1) into the SPSA logic 2) at the constraint provision 3) why this combination allows for variables to go "past" their ideal values.
With all the new Multi-Objective Optimization Algorithms coming out, I was wondering if it would be possible to implement one in OpenMC.
This would almost be synonymous with the
search_for_keffbut with multiple parameters. I tried to create my own adaptation using a combination of distance metrics (https://www.eng.uwo.ca/research/iclr/fids/publications/products/Fuzzy_set_ranking2.pdf (Page 13)) and SPSA yet I won't reach convergence while I am alive because I am running the full-pwr-core example and changing the enrichment of every fuel pin in every assembly (If you would like the files, I could provide them). I could change it to generate the enrichments symmetrically but I'm not there yet. Obviously, this would be the most intense optimization but other things could be attempted to be optimized that are represented by a single scalar quantity like boron concentration, temperature in a certain region, etc.I've looked at Platypus which contains a good number of evolutionary algorithms though it spits out a number of solutions instead just one so sorting through the solutions might take some extra thought.