Feature Request: Multi-Objective Optimization

[SterlingButters]

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_keff but 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.

[paulromano]

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]

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:

• SPSA (Simultaneous Perturbation Stochastic Approximation) https://www.jhuapl.edu/SPSA/
• Compromise Programming https://ageconsearch.umn.edu/bitstream/35695/1/sp99du01.pdf https://pdfs.semanticscholar.org/3f01/6f0b4854abac408e36434aa8df7ba1dc2d71.pdf

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))

[SterlingButters]

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 s

[SterlingButters]

Contents 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 s

[SterlingButters]

In 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]

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)))**p

to

    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)))**p

to 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.

[pranavkantgaur]

Hi @SterlingButters ,

SPSA assumes that the derivatives (df / d_theta) are not available upfront, however, we can use derivative tallies in openMC to estimate say d_k_eff / d_ppm_boron which may significantly reduce the runtime (SPSA's approximate_gradient method) . I have added a demo notebook for this here: openmc-dev/openmc-notebooks#31. Kindly share your comments.

Relevant references:

1. https://dspace.mit.edu/bitstream/handle/1721.1/106690/969775837-MIT.pdf (section 4.2) [Adds tally derivative as additional constraint in curve-fitting]
2. https://www.sciencedirect.com/science/article/pii/S0306454914005404 [Adds synthetic points based on taylor expansion using derivative tallies for curve-fitting]