Minimizing the Standard Deviation of the Thermal Load in the Spent Nuclear Fuel Cask Loading Problem

The paper assumes that, at the end of the operational period of a Spanish nuclear power plant, an Independent Spent Fuel Storage Installation will be used for long-term storage. Spent fuel assemblies are selected and transferred to casks for dry storage, with a series of imposed restrictions (e.g., limiting the thermal load). In this context, we present a variant of the problem of spent nuclear fuel cask loading in one stage (i.e., the fuel is completely transferred from the spent fuel pool to the casks at once), offering a multi-start metaheuristic of three phases. (1) A mixed integer linear programming (MILP-1) model is used to minimize the cost of the casks required. (2) A deterministic algorithm (A1) assigns the spent fuel assemblies to a specific region of a specific cask based on an MILP-1 solution. (3) Starting from the A1 solutions, a local search algorithm (A2) minimizes the standard deviation of the thermal load among casks. Instances with 1200 fuel assemblies (and six intervals for the decay heat) are optimally solved by MILP-1 plus A1 in less than one second. Additionally, A2 gets a Pearson’s coefficient of variation lower than 0.75% in less than 260s CPU (1000 iterations).


Introduction
Optimization problems are found at different levels in electric power systems, from both technological and operations management standpoints. In transport and distribution, not intending to make an exhaustive list, there are studies related to classical modeling and problem solving for electric energy distribution static planning [1][2][3]; other works, also on planning, propose a flexible distribution attending to the emerging variability of demand [4]. There are also many examples of optimization problems at the source point, from renewable generation (see for instance [5] for wind) to nuclear power.
A nuclear power plant is a complex system in which many optimization processes take place. The balance between safety requirements, competition in a liberalized electricity market, and access to limited human and material resources is resolved in nuclear power plants through an approach that prioritizes safety (see, for instance, INSAG-4 [6] and IAEA-TECDOC-1329 [7]).
Nuclear Power Plants (NPP) produce electricity from the heat released in a fission chain reaction; the so-called "nuclear fuel" is actually fissionable material. In Light Water Reactors (LWR), the type of reactor that constitutes the majority of the world's nuclear fleet and the total of the Spanish fleet [8], nuclear fuel is made of uranium (slightly enriched in the 235 U isotope), in the form of ceramic pellets of UO 2 protected and waterproofed by a metallic cladding (zirconium alloy) and grouped together, forming the so-called fuel assemblies (FA) [9]. The fuel used in the two units of Ascó NPP (the plant that has inspired the present study) is manufactured by ENUSA under a license agreement with Westinghouse, the owner of the technology [10].
Presently, the two Ascó units have fuel burn-up cycles of one and a half year. Each 18 months, some 60 FA (out of the 157 contained in the core of the reactor) are replaced (each one containing about 500 kg of uranium). Up to the end of its service, each Spent Fuel Assembly (SFA) has produced around 50 GWd/t of heat.
In the fission process, uranium atoms are transformed into radioactive fission products. Roughly speaking, in an SFA, some 5% of the initial uranium mass has undergone fission. A smaller part of the uranium is also transformed, through neutron capture nuclear reactions, into transuranic elements, which are also radioactive. The hazard of spent fuel is due to the potentially harmful effects of ionizing radiation. Thus, it is necessary to isolate the spent fuel for a long time to avoid the dispersion of the radioactive substances. An ultimate storage solution is proposed in repositories excavated in deep geological formations [11], in a way similar to that followed with chemical contaminants such as mercury [12].
The decay of radioactive substances in the spent fuel causes the so-called "decay heat", which decreases as the activity of those substances declines, but which requires careful management of SFAs to prevent them from being damaged by excess temperature. After being removed from the core of the reactor, spent fuel is temporarily stored in the fuel pool of each power unit. Pool water provides shielding and cooling. Only when the power of the SFA has been sufficiently reduced can it be considered for dry storage in appropriate containers (either temporarily or permanently in a deep geological repository). The lack of political consensus for a definitive solution in Spain made it necessary to propose a Central Interim Storage Facility (CISF) [13], which was never built (now, it is proposed for 2028 [14]).
Over time, the spent fuel pools of the Spanish NPP have been reaching the limit of their capacity. In the case of the two units of Ascó, the pool in each of them can house 1421 FA, which is a figure that includes the 157 positions needed to discharge the entire reactor core if necessary. At the end of 2018, there were 1160 FAs in the Ascó 1 pool and 1104 in the Ascó 2 pool. The plant has 608 SFAs in an Independent Spent Fuel Storage Installation (ISFSI) in 19 HI-STORM containers [15], each hosting 32 SFAs [16].
In Spain, ENRESA (a state-owned company) is responsible for the decommissioning of the nuclear installations, as well as for the ultimate disposal of nuclear and radioactive wastes [17]. At the end of the exploitation period, it is in the interest of the plant operator to be able to transfer ownership of the installation to ENRESA in the shortest possible time. In order to do it, previously, all the fuel must be completely removed from the spent fuel pool [18]. At this moment, the ISFSI can be used as a logistics storage, housing spent fuel containers prior to transportation to the CISF. If the CISF is not yet available, the ISFSI should be usable as a full-scope storage (perhaps even in the long-term).
Thus, an optimization problem arises: it is necessary to remove all of the SFA from the pool within the minimum time and at the minimum cost [19]? Indeed, optimization problems can be found in a wide variety of operations in NPPs. For instance, in the spent fuel pool, where the SFAs are stored in racks, problems arise in optimizing the reliability of the design of the racks and the arrangement of the assemblies in the event of their possible displacement [20], and, as the capacity of the pool is limited to a specific number of assemblies, other management problems of spent nuclear fuel storage arise linked to this limitation [21]. At the right time, SFAs are removed from the pool and relocated in casks, often for long-term storage. The optimization here ranges from determining the SNF removal allocation strategy for a fleet of plants [22], to the design of the containers (e.g., what capacity, how much cooling time needed, etc.) [23], to the strategies to be followed in the loading of the casks.
As for the latter, Zerovnik et al. [24] propose combinatorial methods to optimize the cask loading, aiming at final geological repository. Considering as well the final deposition, Ranta and Cameron [25] used a two-phase heuristic method in order to minimize the thermal load of the canister with the largest load. Gomes et al. [26] propose a methodology, using genetic algorithms, to determine the optimal time to load the casks for dry interim storage, taking into account several cost factors. Spencer et al. [27] have develop a method (embedding greedy randomized adaptive search procedures in a memetic algorithm) to optimize the cask loading in a multi-objective search of configurations that minimize the number of casks, the thermal load, and the time when the transportation requirements are met.
In the present work, the cask loading problem is solved with two hierarchical optimization criteria: (i) minimize the cost of the regionalized casks necessary to relocate the spent fuel located in the pool at a given moment [28] using mixed integer linear programming [29]; and (ii) minimize the standard deviation of the thermal load in the set of casks using a multi-start metaheuristic with local search procedures [30].
In Section 2, we describe some dry storage systems for spent nuclear fuel. Section 3 is devoted to the description, modeling, and resolution of the spent nuclear fuel cask loading problem. Section 4 illustrates the procedure proposed in a case study inspired in Ascó NPP. Finally, Section 5 provides the conclusions derived from this work and a brief description of future research lines.

Dry Storage Systems for Spent Nuclear Fuel
Different solutions exist for the interim dry storage of spent nuclear fuel; as a comparison, the 14 closed sites in the USA use 11 different storage systems from four different suppliers [31].
In Spain, ENRESA has adopted different solutions. Some plants are using dual-purpose containers (transport and storage) whose structure provides shielding, confinement, and heat dissipation. Ascó's two units are using the multi-purpose HI-STORM system, which is described below [15,16].
The spent fuel storage and transport system in Ascó is made up of different elements: the multi-purpose canister (MPC-32), the storage module (HI-STORM), the transport overpack (HI-STAR), and the transfer cask (HI-TRAC), which allows the MPC to be transferred to HI-STORM and HI-STAR [15] (see Figure 1).
Energies 2020, 13, x FOR PEER REVIEW 4 of 29 Figure 1. Components of the HI-STORM system for the storage and transportation of spent nuclear fuel. The HI-TRAC transfer cask allows the MPC to be transferred to the HI-STORM and HI-STAR [15]. Components of the HI-STORM system for the storage and transportation of spent nuclear fuel. The HI-TRAC transfer cask allows the MPC to be transferred to the HI-STORM and HI-STAR [15].
The time required to completely remove the spent fuel from the pool depends on the minimum cooling time necessary for the last SFA to be loaded into a cask. This cooling period (that is, the time elapsed since the unload of the SFA from the reactor core) depends, in turn, on the type of container used and the technical specifications approved by the regulator for it.
To guarantee the integrity of the SFAs, the temperature of the metal cladding that surrounds the UO 2 pellets must be limited. This fact imposes requirements of sufficient heat dissipation in the cask design. However, also, it imposes a limit to the thermal power that the SFAs can generate inside the cask. For instance, in the HI-STORM system, under uniform charging conditions, each of the positions of the MPC-32 can house an SFA with a maximum power of 0.9375 kW at the time of loading (i.e., the total thermal power dissipated by the canister will be limited to 30 kW) [32]. The system allows regionalized loads; so, SFAs with higher decay heat can be located in the central region of the canister (region R a , formed by 12 positions), at the cost of limiting the maximum power per position in the peripheral region (region R b , 20 positions) and also limiting the total thermal load to a value less than 30 kW [33]. Figure 2 represents the positions assigned to each region.  In fact, the loading of a cask is constrained by the technical specifications in a more complex way, since they take into account, in addition to the cooling time, multiple fuel characteristics such as decay heat, burn-up, initial enrichment, the possibility that the assembly has some damage, or the possible presence of elements such as neutron sources, control rods, burnable poisons, etc. These conditions vary according to the type of cask, its use (storage or transport), and other specific limitations. In fact, the loading of a cask is constrained by the technical specifications in a more complex way, since they take into account, in addition to the cooling time, multiple fuel characteristics such as decay heat, burn-up, initial enrichment, the possibility that the assembly has some damage, or the possible presence of elements such as neutron sources, control rods, burnable poisons, etc. These conditions vary according to the type of cask, its use (storage or transport), and other specific limitations.
In what follows, it is assumed that at the end of the operational life of the plant, an ISFSI will exist that allows transfering all of the SFAs to it using the HI-STORM system (without imposing the present transport requirements, which are more restrictive than storage requirements). It is considered here that the storage requirements apply only to the decay heat of each SFA (calculated in turn from the burn-up, cooling time, and initial enrichment). The minimum time for the complete removal of all of the SFA from the pool is bounded by the moment when the hottest SFA fits into the requirements. The cost will depend on the number of containers to be used. So, the number of containers needed to remove all the SFA from the spent fuel pool needs to be minimized: this is the optimization problem addressed in the present study.

The Spent Nuclear Fuel Cask Loading Problem
In this section, we describe and formalize the spent nuclear fuel cask loading problem, and we propose a hybrid metaheuristic to solve it. Our method works in three phases using mixed integer linear programming (MILP), a deterministic algorithm, and local search heuristics.

Description of the Problem
At a given instant T, corresponding to the planned date for transferring in one go the contents of the spent fuel pool to the ISFSI, there are n SFAs, forming a set I (i = 1, . . . , n (n ≡ |I|)).
At the time of removal, T, each fuel assembly i ∈ I is characterized by a thermal power (decay heat), q i (T).
On the other end, there is an ISFSI ready to host all of the casks needed. Casks are regionalized, i.e., they are compartmentalized into regions with attributes of capacity and decay heat. Thus, there is a set of virtual regions J ( j = 0, 1, . . . , m (|J| ≡ m + 1)), in which each type of region ( j ∈ J) is characterized by the upper and lower bounds on the decay heat a single SFA may have q − j , q + j , and the number of SFAs in that region that can be held by a single cask, h j . j = 0 symbolizes the type of universal region capable of housing any fuel assembly without restrictions.
There is a set K (k = 0, 1, . . . , l (|K| ≡ l + 1)) of cask classes; k = 0 symbolizes the universal class, which is capable of accommodating any fuel assembly without restrictions at an arbitrary high cost γ 0 → ∞ . Each class is defined in relation to the maximum thermal load allowed in each region of the cask.
An adversity coefficient γ k is associated to each class of cask (k ∈ K), corresponding to the arbitrary cost of selecting a cask of class k ∈ K to house at least one SFA. We also assume for the class k ∈ K a maximum number of units available, d k (∀k). The total number of available casks is D = ∀k d k .
The relationship between a virtual region type ( j ∈ J) and a cask class (k ∈ K) is given by matrix B, whose coefficients b j,k (∀j∀k) take the value 1 when the region of type j ∈ J is associated to the cask of class k ∈ K, and 0 otherwise.
On the other hand, the relationship between the fuel assembly i ∈ I and the type of region j ∈ J is given by matrix A(T), whose coefficients, a i,j (T) (∀i∀j), are equal to 1 when the fuel assembly i ∈ I can be located in the region j ∈ J at time T, and 0 in any other case; i.e., a i,j (T) = 1 ⇔ q i (T) ∈ q − j , q + j . It follows from it that the relationship between the fuel assembly i ∈ I and the class of cask k ∈ K is given by matrix C(T) resulting from the boolean product of matrices A(T) and B; that is, C(T) = A(T) ⊗ B, thus fulfilling the following property: (a i,j (T) = 1 b j,k = 1) ⇒ c i,k (T) = 1 .
Under such conditions [28], the spent nuclear fuel cask loading problem consists of transferring the maximum number of the fuel assemblies, which are located in the pool of a nuclear power plant on a date T and expressed as time after the end of operation, into a set of regionalized casks of different classes, meeting the restrictions and minimizing costs (i.e., minimizing the total adversity).

Solving the Problem
As discussed above, our proposal solves the problem in three phases. The first phase uses MILP to allocate fuel assemblies to virtual regions with minimum cost (i.e., a minimum number of casks or minimum total adversity) [28]. In the second phase, a family of deterministic algorithms has been designed to assign each SFA to a specific cask actual region. Finally, in the third phase, a local search heuristic is proposed to minimize the standard deviation of the thermal load of the casks set.

Phase 1: MILP Model to Minimize the Cask Loading Cost
To formalize the proposed resolution methods, we use the following nomenclature. Parameters: T Time or expected date to remove in one go the SFA from the pool in the nuclear power plant, expressed as time after the end of operation. I Set of SFA located in the pool of the nuclear power plant at time T. Let the elements be i = 1, . . . , n (n ≡ |I|). J Set of virtual region types according to the upper and lower bounds on the decay heat of the SFAs. Let the regions be j = 0, 1, . . . , m with |J| ≡ m + 1. The universal region, j = 0, does not impose any limit; for practical purposes, the assigned elements to the universal region will remain in the pool. K Set of classes of casks: k = 0, 1, . . . , l with |K| ≡ l + 1. The type of cask considered here corresponds to MPC-32 described above. We will assume that there is a type of universal container, k = 0, associated with the universal region j = 0; for practical purpose, the elements assigned to the universal container will remain in the pool.
Decay heat of the fuel assembly i ∈ I at the time T. q + j Upper limit of a single SFA's decay heat in the virtual region j ∈ J. q − j Lower limit of a single SFA's decay heat in the virtual region type j ∈ J. Here, we will assume q − j = 0 ∀ j ∈ J. A(T), a i, j (T) Element-region compatibility matrix; technological coefficients a i, j equal to 1 if the fuel assembly i ∈ I can be located in the virtual region j ∈ J at time T, and 0 otherwise. Formally: Region-cask compatibility matrix; technological coefficients b j,k equal to 1 if the region of type j ∈ J is associated with the cask class k ∈ K, and 0 otherwise. C(T), c i,k (T) Element-cask compatibility matrix, which is defined as C(T) = A(T) ⊗ B; technological coefficients c i,k (T) are equal to 1 when the element i ∈ I can be placed in a cask of class k ∈ K at time T, and 0 otherwise. Therefore, it is true that: h j Storage capacity of a region of type j ∈ J measured by the number of SFAs a single cask can hold in that region at most.

H k
Maximum storage capacity of a cask of class k ∈ K measured by the number of fuel assemblies.
Number of available casks of the class k ∈ K at time T. The expected total cask availability is D = ∀k d k . γ k Adversity coefficient for a cask of class k ∈ K. The value γ k symbolizes the arbitrary cost if a cask of class k ∈ K is selected to hold SFAs. Variables: x i,j Binary variable that equals 1 if the fuel assembly i ∈ I is assigned to the virtual region j ∈ J, and 0 otherwise. y j Number of fuel assemblies assigned to the virtual region j ∈ J. z k Number of casks of class k ∈ K.

Γ(T)
Total adversity function. It represents the overall cost of the casks used to hold the spent nuclear fuel on the scheduled date, T, after the pool of the nuclear power plant is emptied.
Using these parameters and variables, we formulate the following mathematical model.
Subject to: x i,j ∈ {0, 1}∀i = 1, . . . , n ∀j = 0, 1, . . . , m In the MILP-1 model, the objective Function (1) represents the minimization of the cask loading cost or total adversity, whose value coincides with the total number of casks when all the fuel assemblies are removed from the pool (z 0 = 0) and the adversity coefficients are unitary for any class of cask (γ k = 1 ∀k = 1, . . . , l). Equalities (2) and (3) guarantee the assignment of all the fuel assemblies to any region, including the universal region j = 0. Equality (4) determines the number of fuel assemblies that must be assigned to each type of region j ∈ J. Constraint (5) lower bounds the number of casks of each class k ∈ K, while Constraint (6) upper bounds that value according to their availability. Condition (7) imposes that the decision variables (x i,j ) are binary and, finally, the non-negativity and integrity of the auxiliary variables y j and z k are given by (8) and (9), respectively.
Running the MILP-1 model offers three types of results: (i) an optimal assignment of the fuel assemblies to the various virtual regions (x * i,j ∈ {0, 1} ∀i∀j), (ii) the optimal number of SFA assigned to each region (y * j ∀j), and (iii) the optimal number of casks of each class (z * k ∀k) corresponding to the minimum cost.
Although the previous results are interesting, the location of the fuel assemblies in the real regions of the casks remains to be determined.

Phase 2: Assignment of the Fuel Assemblies to the Real Cask Regions
Any sequence of objects corresponding to the triplet of sets (I, J, K) can lead to a different solution. For example, if the sequence is built on the set I of fuel assemblies, then the number of different ways to assign is n! This can lead to an astronomical number of potentially different solutions for instances of realistic dimensions (e.g., 1164! in our case study). Based on this idea, we propose a family of constructive algorithms, A1[π(n)], which depends on the regulatory sequence π(n) of the assignment of fuel assemblies to virtual regions.
Starting from an optimal solution given by MILP-1(x * i,j ∀i∀j, y * j ∀j, z * k ∀k), the procedure A1[π(n)] creates a sequence π(n) of fuel assemblies (Statement 5 in A1[π(n)]) to establish the sequence of assignment of SFAs to casks. Once the procedure determines the class of the cask k * t that corresponds to the current element π t , the element is recorded in the matrix MPC32[·] (Statement 13 in A1[π(n)]), whose rows correspond to casks and whose columns correspond to positions within the cask. Subsequently, the heat loads Q mpc32 [·] in the casks are determined (Statement 16 in A1[π(n)]), and the elements are sequenced in the casks according to such loads (Statement 17 in A1[π(n)]). Finally, the results of the algorithm are returned (Statement 18 in A1[π(n)]). Algorithm 1. A1[π(n)]: Algorithms to assign fuel assemblies to regions according to π(n) // Create assignment sequence of SFAs 5. π(n) = (π 1 , π 2 , . . . , π n ) ← Generate_sequence(I, n, P κ ) 6.
// Find the virtual region assigned to element π t . Let j *

14.
End while 15. // Calculate residual heat loads and order elements in matrix MPC32 Note that Algorithm 1 A1[π(n)] can generate many solutions thanks to the fact that the starting point is a permutation π(n) of the fuel assemblies of the set I (Statement 5 in A1[π(n)]). Such permutations can be created using a family of rules P κ applied on the set I, which is initially ordered according to the SFA's identifier (i ∈ I). The family of rules P κ can be of diverse nature; they may consist of the same initial sequence (lexicographic order) or the simple sequence of the elements according to the decay heat or in more complex rules used for the generation of neighborhoods (exchange and insertion) in algorithms based on trajectories through the search domain.
Obviously, all the solutions generated by A1[π(n)] are optimal when the function objective is the minimization of the cost of the casks [28], as happens in Phases 1 and 2 of this work; however, when the objective function corresponds to another optimization criterion, such as the balancing of the heat load in the set of casks, A1[π(n)] can be used in, at least, two ways: (1) In the construction phase (or phases) of multi-start algorithms, v.gr. greedy randomized adaptive search procedure (GRASP) algorithms, to build initial solutions to be used as seeds for the next phase of trajectory-based improvement in the search domain. (2) In the systematic generation of solutions, to assist metaheuristics based on populations (v.gr. genetic and memetic algorithms).
In this work, Algorithm 1 A1[π(n)] constitutes one of the fundamental pillars for the construction of initial solutions, and the chosen way is (1).

Phase 3: Reduction of the Dispersion of the Thermal Load of the Casks
Phase 3 of the proposed procedure aims to reduce, by a predetermined number of iterations, the dispersion of the thermal load (measured in kW) in the set of casks. The dispersion metric to minimize chosen here is the standard deviation σ(Q, T); this is: where: Total optimal number of casks: z * tot = ∀k z * k . It is obtained by applying the MILP-1 procedure by setting values to the adversity coefficient γ k (∀k ∈ K). Q δ k Total thermal load due to the SFA loaded on date T in the δ k -th cask of class k: . These values are obtained after applying the Algorithm 1 A1[π(n)] by fixing a permutation π(n) of the elements of set I. Q avg Average value of the thermal load in the set of casks on date T: Standard deviation of the thermal load in the set of casks on date T.
In order to carry out Phase 3, the non-exhaustive descent Algorithm 2 A2[π(n)] is applied, based on local search optimization (see Algorithm 2 A2[π(n)] below). Its main characteristics are as follows.
The Algorithm 1 A1[π(n)] is applied to the solution X 0 , considering the arbitrary sequence π(n), which initially corresponds to π 0 (n): the lexicographic order of the SFA's identifiers (i ∈ I).
Thus, the extended solution X 0 mpc = X 0 , π(n), MPC32 z * k , 32 is obtained, which contains detailed information on the location of each element inside the casks z * k (∀k) in the matrix MPC32 z * k , 32 . • Solution X 0 mpc is fixed as the incumbent solution of the local search iterative procedure X c mpc ← X 0 mpc . In addition, X 0 mpc is fixed as the best solution (X * mpc ← X 0 mpc ) , and the standard deviation of this solution σ X * mpc is determined.
Phase 3.2: Improvement of the initial solution with A2[π(n)] • A tentative solutionX c mpc belonging to the restricted neighborhood of the incumbent solution X c mpc is generated-that is,X c mpc ∈ N X c mpc , P LS , S F -where P LS is the rule for generating neighboring solutions and S F is the set of possible solutions. Specifically, given the incumbent solution X c mpc , the generation rule P LS operates as follows: for all pairs of different casks (δ k , δ k ) and all pairs of positions in these casks (pos(δ k ), pos(δ k )), the tentative neighboring solutions X c mpc exchange the elements between these positions, within the set of possible solutions S F that is regulated by the SFA-cask compatibility matrix C(T).

•
If the tentative solutionX c mpc is a possible solution (X c mpc ∈ S F ) and σ X c mpc < σ X * mpc , then the solutionX c mpc becomes the best solution X * mpc ←X c mpc and also the new incumbent solution X c mpc ←X c mpc .

•
The above steps are repeated within Phase 3.2 as long as there is an improvement to the incumbent solution or the number of tentative solutions generated does not reach the preset valueît ∞ (v.gr. it ∞ = 1000 tentative solutions).
Algorithm 2. A2[π(n)]: Non-exhaustive local search descent algorithm based on π(n) Repeat for all X ∈ N X c mpc , P LS ,

Multi-Start Metaheuristics Assisted by Assignment Sequences
The proposed metaheuristic integrates the procedures of the three phases described above (see Schema MS X 0 , π(n) ). This metaheuristic can be classified within the category of multi-start algorithms, since it uses a multitude of initial solutions generated by the algorithm A1[π(n)] of Phase 2. In the generation of solutions, A1[π(n)] uses a current sequence string (π c (n)); this chain is progressively constructed by moderately altering the current sequence giving rise to the next ongoing permutation. The main characteristics of MS[·] are the following:

•
The Algorithm 3 metaheuristic MS X 0 , π(n) ) is initialized by executing Phase 3.1, starting with an initial solution X 0 offered by MILP-1 (Phase 1) and registering the sequence π(n) as an assignment current sequence; that is: π c (n) ← π(n) .

•
In a specific iteration, given an arbitrary current sequence, which we will denote as π c (n), the algorithms A1[π(n)] and A2[π(n)] are cascaded, making π c (n) ← π(n) , and thus obtaining a solution X * mpc (π c (n)) that is locally optimal.

•
In an iteration, to go from an ongoing assignment current sequence (π c (n)) to the next current sequence π c (n) , there are many alternatives. Specifically, in this work, we have chosen to apply a slight mutation on the first sequence π c (n). This mutation, characterized by the P m rule, consists of exchanging, in pairs, the fuel elements located in four positions of the current sequence chosen at random. • Thus, for each X 0 , the number of initial solutions locally optimized by MS X 0 , π(n) is equal to the number of sequences generated with the P m mutation rule. The MS[·] procedure ends when the number of mutated sequences reaches a maximum value equal to it π∞ (v.gr. it π∞ = 1000 mutate sequences).

Case Study
This section is dedicated to a case study inspired in the spent fuel pools of Ascó NPP (Tarragona, Spain). It contains the following points: (i) the model assumptions, (ii) the requirements for spent fuel cask loading, (iii) the procedures used, (iv) the data used, and (v) the results achieved.

Hypothesis for the Case
The following hypotheses are assumed for the case at hand: Hypothesis 1 (H1). Fuel assemblies are characterized by their type, burn-up (MWd/t), initial enrichment (% by mass), and the time elapsed since they are unloaded from the core (in years). From these data, the decay heat of each element can be determined at the time it is transferred from the pool to the cask.

Hypothesis 2 (H2).
There are no fuel assemblies damaged or affected by other additional constraints.

Hypothesis 3 (H3).
The elements associated with the fuel, such as neutron sources, control rods, consumable poisons, and others, are not considered in the model.

Hypothesis 4 (H4).
It is assumed that only one type of cask is available in terms of physical dimensions and number of positions for fuel assemblies, the MPC-32.

Hypothesis 5 (H5).
The fuel assemblies are loaded into the casks taking into account the requirements for storage in module-type HI-STORM storage containers. It is assumed that there is no capacity limit in the ISFSI. Hypothesis 6 (H6). The contents of the spent fuel pool at the end of the operation of the plant (T = 0) are assumed to be known. To carry out the case study, the real content of the spent fuel pool of NPP Ascó 1 is used here. The decay heat of each of the SFA at any given time T can be computed from their cooling time (this is T plus the time each SFA has previously spent in the pool) and the SFA's burn-up and initial enrichment.

Hypothesis 7 (H7).
The pool should be emptied using the minimum possible number of casks on the earliest possible date. The date has been chosen as T = 6 years in order to guarantee that even the hottest SFA meets the requirements.

Hypothesis 8 (H8).
It may be interesting that the thermal load of the casks is balanced. Therefore, we aim at minimizing the standard deviation of the thermal load of the set of casks.

Requirements for the Storage of Fuel Assemblies
The requirements to take into account for emptying the pool are as follows: R1. Casks are of the MPC-32 type, which are compartmentalized in 32 positions, each one containing a single fuel assembly (see Figure 2). R2. The positions of a MPC-32 are grouped into two regions: (R a ) a permissive region with a high limit for the decay heat formed by h a = 12 positions and (R b ) a restrictive region with a low limit for the decay heat formed by h b = 20 positions. R3. A fuel encapsulation position (concept for which we will use the Latin noun locus-loci) can house an SFA whose cooling time is greater than 5 years (in our case, with T = 6 years, this does not constitute any limitation). R4. The contents of each locus in any cask must have a decay heat less than or equal to the maximum allowed in the locus: q + a in region R a and q + b in region R b . R5. The maximum heat load allowed in a cask is equal to Q + MPC = 30 kW. (1) Set a value of the ratio of the maximum values of decay heat per locus ρ = q + a /q + b , such that 1 ≤ ρ ≤ 3 (ρ = 1 corresponts to the uniform loading, ρ = 3 is the maximum value allowed by the regulator for this ratio). (2) Determine the maximum allowed thermal load per locus, q + b = q + b (ρ, h a , h b ) for the restrictive region R b , from the expression: where f (ρ, h a , h b ) is a function bounded to 1 for ρ = 1 and decreasing with ρ. The function f (ρ, h a , h b ) for the MPC-32 is: (3) Determine the maximum allowed thermal load per locus, q + a (ρ, h a , h b ), for the permissive region R a :

Data and Other Characteristics of the Case
Obviously, the dimensions of the MILP-1 model are a function of the instances used. In this work, we use the ASCÓ#01_6.0 instance (see Table A1 in Appendix A) as data for decay heat. The main characteristics of the experiment are as follows: (1) The number of spent fuel assemblies is equal to n = 1164. The elements are identified with the indices i ∈ I corresponding to the n first natural numbers. (2) The expected date for the removal of elements from the pool is T = 6 years after the end of operation. The decay heat values q i (T) ∀i ∈ I listed in Table A1 of Appendix A correspond to T = 6 years. (3) The number of virtual regions of decay heat is equal to m = 6.
Let r j = q ∈ R : q ∈ q − j , q + j (∀j = 0, 1, . . . , 5) be those regions that are defined with closed intervals of decay heat per locus. The regions are specified from three values of the ratio of maximum thermal loads: ρ = 1, 2, 3. Let: For convenience, virtual regions have been numbered in descending order of the upper limit q + j ∀ j of each interval.
(1) The definition of virtual regions (r 0 to r 5 ) allows establishing the element-region compatibility matrix A(T), whose technological coefficients adopt the value 1 if the fuel assembly i ∈ I can be assigned to a virtual region of type j ∈ J at the time T, and 0 otherwise. Formally: . . , 1164; ∀j = 0, 1, . . . , 5 (2) The number of cask classes is l = 3 if we omit the neutral virtual region r 0 , which corresponds to the pool. The correspondence between the virtual regions (r 1 to r 5 ) and the real classes of the regionalized cask (CR 1 a CR 3 ) is established through the region-cask compatibility matrix B, whose coefficients b j,k are equal to 1 when the region r j ( j = 1, . . . , 5) is associated with cask CR k (k = 1, . . . , 3), and they are 0 otherwise. In our case, matrix B is as follows: Thus, a cask CR 1 is compartmentalized in a permissive region with the characteristics of r 1 and a restrictive region associated with r 5 ; a cask CR 2 is associated with a permissive region r 2 and another restrictive r 4 ; and finally, a cask CR 3 has only one region with the qualities of the intermediate region  (3) Taking into account the values and characteristics exposed in the previous points, the dimensions of the MILP-1 model for the Ascó # 01_6.0 instance are the following: (i) 6984 binary decision variables (x i,j ) corresponding to the assignment of fuel assemblies to virtual regions; (ii) 6 integer variables (y j ) associated with the number of fuel assemblies hosted by each virtual region; (iii) 4 representative integer variables (z k ) for the number of casks of each class; and (iv) 1195 constraints (1199 if the availability of casks is bounded).

Procedures
All compiled codes corresponding to the procedures used have been executed on a DELL Inspiron-13 (Intel(R) Core(TM) i7-7500U @ 2.70 GHz CPU 2.90 GHz, 16 GB de RAM, x64 Windows 10 Pro). The models and tools used to find the optimal solutions for this case study are the following: (1) Model MILP-1 [28] based on MILP that allows minimizing the cost (Γ * (T)) of the containers required (z * k ∀k) to transfer all the SNF from the NPP pool. MILP-1 also offers the number of SFAs assigned to each region (y * j ∀j), as well as an optimal assignment (x * i,j ∀i∀j) of each SFA to a virtual region. For the exploitation of MILP-1, we use the IBM ILOG CPLEX solver (Optimization Studio v.12.2, win-x86-64).
(2) Deterministic method A1[π(n)] that starts from the result of MILP-1 and offers an assignment of SFAs to the permissive and restrictive regions of the casks based on the sequence of assignment of SFAs to casks π(n). Here, 1000 consecutive allocation sequences are considered, which are generated by slightly mutating one of them to obtain the next, following the neighborhood transition rule P LS . MILP-1 plus A1[π(n)] (with lexicographic order of the fuel assembly codes) constitute the deterministic procedure DP[π(n)] proposed in [28]. (3) Local search heuristic A2[π(n)] that starts from a solution generated by A1[π(n)] and explores neighborhoods non-exhaustively until reaching a local optimum. The procedures MILP-1 plus A1[π(n)] plus A2[π(n)] (with the mutation rule P m applied to the π(n) sequence in it π∞ iterations) constitute the metaheuristic MS X 0 , π(n) proposed in this work.
In short, metaheuristic MS X 0 , π(n) works in this manner: from an optimal solution found by MILP-1, which is characterized by the triple of sets of values [(x * i,j ∀i∀j), (y * j ∀j), (z * k ∀k)], the algorithm A1[π(n)] is applied it π∞ times, transferring one by one and following it π∞ sequences π(n), the SFAs assigned to the virtual regions r j (∀j) toward the physical regions (permissive and restrictive) of the casks CR k (∀k), using the region-cask compatibility matrix B. Afterwards, the local optimization corresponding to heuristics A2[π(n)] is applied to all (it π∞ ) solutions generated by algorithm A1[π(n)].

Results
In addition to considering the solution from MILP-1 with the triplet → γ 0 = (1, 1, 1), which corresponds to minimizing the total number of casks, and in order to assess the impact of the adversity coefficients γ k (∀k) on the solutions offered by MILP-1, 24 other optimal solutions are also analyzed, fixing the list of arbitrary costs → γ = (γ 1 , γ 2 , γ 3 ) to the following values:  Using the deterministic procedure DP[π(n)] proposed in [28], after obtaining an optimal solution with MILP-1 [(x * i,j ∀i∀ j), (y * j ∀j), (z * k ∀k)], the algorithm A1[π(n)] is applied following the sequence π(n) corresponding to the lexicographic order of the fuel assembly identifiers. The global results of the 25 solutions are shown in Table 1. Table 1. Summary of results of the 25 optimal initial solutions (#0 to #24) provided by procedure DP[π(n)] proposed in [28] for the Ascó #01_06 instance with the following technical characteristics: (i) 1164 Spent Fuel Assemblies (SFAs), (ii) T = 6 years, (iii) 6 virtual regions (including the neutral region) defined by three values of the ratio ρ (1, 2, and 3), and (iv) 4 cask classes (including the neutral container). Specifically, in Table 1, the following data are collected for each solution (#0 to #24): (1) The cost → γ = (γ 1 , γ 2 , γ 3 ) imputed to each cask class. (2) The optimal number of casks of each class (z * k ∀k) and the total number (z * tot ) of casks required by each solution.  Table 1 shows 16 optimal solutions that require 37 MPC-32, which is a number that also coincides with the minimum number of casks, according to the lower bound ( 1164/32 ). However, these solutions show different values of the standard deviation of the thermal load, depending on the configuration of the 37 casks.
As an example, initial solution #9, one with the lowest thermal load dispersions between the casks: σ(Q) = 2.90, corresponds to the triplet of costs → γ 9 = (3, 2, 3), in which it is considered that the adversity to the use of casks of classes CR 1 and CR 3 is 1.5 times higher that of class CR 2 casks.
Based on solution #9, the FAs column in Figure 3 shows the number of SFAs assigned to each virtual region, i.e., y * 1 = 24, y * 2 = 274, y * 3 = 369, y * 4 = 457, y * 5 = 40. This is translated into a total requirement of 37 casks (see the MPC[z k , h j ] column in Figure 3), z * 1 = 2, z * 2 = 23 and z * 3 = 12 regionalized casks of classes CR 1 , CR 2 , and CR 3 , respectively.     In view of Figure 4, the following can be stated: • Cask n.1 of initial solution #9 is class with a permissive region ← of 12 positions and a restrictive region ← of 20. The 12 hottest assemblies in this cask are located in region , so locus-1 is occupied by the fuel assembly 751 with a decay heat of 1.3740 kW, while the fuel assembly 753 (the last in ) is located in locus-12 with 1.3234 kW. Complementarily, 20 "cold" elements are located in positions 13 to 32 of the restrictive region of MPC-32 n.1, assembly 23 in locus-13, and assembly 21 in locus-32, whose decay heats are 0.4189 kW, coinciding with the decay heat of all the elements located in this region. On the other hand, Figure 4 shows the detailed loads of the containers MPC-32 number n.1 (class CR 1 ), number n.3 (class CR 2 ) and number n.26 (class CR 3 ) of the initial solution #9. These loads are characterized by the identifiers and decay heat of the fuel assemblies located in each MPC-32.     In view of Figure 4, the following can be stated: • Cask n.1 of initial solution #9 is class with a permissive region ← of 12 positions and a restrictive region ← of 20. The 12 hottest assemblies in this cask are located in region , so locus-1 is occupied by the fuel assembly 751 with a decay heat of 1.3740 kW, while the fuel assembly 753 (the last in ) is located in locus-12 with 1.3234 kW. Complementarily, 20 "cold" elements are located in positions 13 to 32 of the restrictive region of MPC-32 n.1, assembly 23 in locus-13, and assembly 21 in locus-32, whose decay heats are 0.4189 kW, coinciding with the decay heat of all the elements located in this region. In view of Figure 4, the following can be stated:

•
Cask n.1 of initial solution #9 is class CR 1 with a permissive region R a ← r 1 of 12 positions and a restrictive region R b ← r 5 of 20. The 12 hottest assemblies in this cask are located in region R a , so locus-1 is occupied by the fuel assembly 751 with a decay heat of 1.3740 kW, while the fuel assembly 753 (the last in R a ) is located in locus-12 with 1.3234 kW. Complementarily, 20 "cold" elements are located in positions 13 to 32 of the restrictive region R b of MPC-32 n.1, assembly 23 in locus-13, and assembly 21 in locus-32, whose decay heats are 0.4189 kW, coinciding with the decay heat of all the elements located in this region. • Casks n.3 (initial solution #9) is class CR 2 , with its permissive region R a ← r 2 and its restrictive region R b ← r 4 . In region R a , the 12 hottest elements of this cask are located with element 255 occupying locus-1 and element 281 located in locus-12 with residual heats equal to 0.7813 kW and 0.6320 kW, respectively. In addition, all the SFAs located in the restrictive region R b are considered to have a decay heat of 0.4189 kW. On the other end, to illustrate the performance of the MS X 0 , π(n) metaheuristic in terms of improving the standard deviation of the thermal load in the set of casks, nine solutions have been selected from the 25 solutions presented in Table 1. The solution selection criterion meets the following conditions: (1) Have a total number of z * tot containers that is minimal; that is: z * tot = 37. (2) If two or more solutions with z * tot = 37, have the same configuration of containers by type (z * 1 , z * 2 , z * 3 ), the one that presents smaller standard deviation σ(Q) is chosen; for example, between solutions #2 and #9, we select solution #9.
The final solutions selected to apply Phase 3 to improve σ(Q) using metaheuristic MS X 0 , π(n) are shown in Table 2. The indicators selected to compare these solutions are as follows: (i) the standard deviations σ A1 (Q) and σ A2 (Q) corresponding to the initial solutions provided by A1[π(n)] and the improved solutions provided by A2[π(n)], (ii) the iteration it * corresponding to the best solution with MS X 0 , π(n) ; (iii) the computing time in seconds, CPU it * , corresponding to the iteration it * ; and (iv) the computing time in seconds, CPU 1000 , given to MS X 0 , π(n) to improve 1000 initial solutions provided by A1[π(n)] from a solution X 0 provided by MILP-1.
In view of Table 2, we can draw the following conclusions: (6) On average, the standard deviation σ(Q) is reduced approximately 20 times by applying the metaheuristic MS X 0 , π(n) versus procedure DP[π(n)] proposed in [28]. Table 2. Summary of results of the solutions selected to apply Phase 2 plus Phase 3 of the proposed procedure MS X 0 , π(n) . The table shows the values: σ A1 (Q), σ A2 (Q), it * , CPU it * , and CPU 100 . To illustrate the reduction in the dispersion of the thermal load in the set of casks, Figure 5 shows the detailed loads of the containers MPC-32 number n.1 (type CR 1 ), number n.3 (type CR 2 ), and number n.22 (type CR 3 ) from solution #10 provided by metaheuristic MS X 0 , π(n) .   Figures 6 and 7 show graphically the thermal loads corresponding to the solutions provided by procedure DP[π(n)] proposed in [28] (see Figure 6) versus metaheuristic MS X 0 , π(n) proposed in this work (see Figure 7) for each region (permissive and restrictive) and for each of the 37 casks. Note the significant differences that exist between the thermal loads of the 37 used MPC-32 in both figures, despite the fact that the containers present the same thermal limitation characteristics number by number (1 to 37). Figures 6 and 7 show graphically the thermal loads corresponding to the solutions provided by procedure [ ( )] proposed in [28] (see Figure 6) versus metaheuristic [ , ( )] proposed in this work (see Figure 7) for each region (permissive and restrictive) and for each of the 37 casks. Note the significant differences that exist between the thermal loads of the 37 used MPC-32 in both figures, despite the fact that the containers present the same thermal limitation characteristics number by number (1 to 37).

Conclusions
A hybrid multi-start procedure assisted by mixed integer linear programming and local search optimization algorithms, that we have called here metaheuristic [ , ( )], has proven to be competitive in solving the spent fuel cask loading problem with minimal storage cost and a minimum standard deviation of thermal load. This is confirmed through the experiments carried out with the   Figure 6. Solution #10 (procedure DP[π(n)] proposed in [28]): thermal loads for each of the 37 used MPC-32 as an addition of the ones in the permissive (Qa) and restrictive (Qb) regions.  [28] (see Figure 6) versus metaheuristic [ , ( )] proposed in this work (see Figure 7) for each region (permissive and restrictive) and for each of the 37 casks. Note the significant differences that exist between the thermal loads of the 37 used MPC-32 in both figures, despite the fact that the containers present the same thermal limitation characteristics number by number (1 to 37).

Conclusions
A hybrid multi-start procedure assisted by mixed integer linear programming and local search optimization algorithms, that we have called here metaheuristic [ , ( )], has proven to be competitive in solving the spent fuel cask loading problem with minimal storage cost and a minimum standard deviation of thermal load. This is confirmed through the experiments carried out with the  3  4  5  6  7  8  9  10  11  12  13  14  15  16  17  18  19  20  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  Thermal  In addition, Appendix B contains the detailed configurations of the 37 containers corresponding to solution #10 provided by metaheuristic MS X 0 , π(n) .

Conclusions
A hybrid multi-start procedure assisted by mixed integer linear programming and local search optimization algorithms, that we have called here metaheuristic MS X 0 , π(n) , has proven to be competitive in solving the spent fuel cask loading problem with minimal storage cost and a minimum standard deviation of thermal load. This is confirmed through the experiments carried out with the database of Ascó NPP. It allows generating instances with industrial dimensions (e.g., 1200 fuel assemblies and from 1 to 30 classes of regionalized casks, which translates into a number of decay heat virtual regions between 2 and 60).
The method proposed here, MS X 0 , π(n) , to solve the single-stage cask loading of spent nuclear fuel is computationally competitive in its three phases: MILP-1 mathematical program, deterministic algorithm A1[π(n)], and local search heuristic algorithm A2[π(n)].
(i) The MILP-1 model is used to minimize the cost of casks required to host the nuclear fuel; it offers, as a solution, an optimal assignment of fuel assemblies to virtual thermal regions, in less than 0.5 CPU seconds, using instances with 1200 fuel assemblies, 6 types of regions, and 4 classes of casks.
(ii) Algorithm A1[π(n)], which starts from an MILP-1 solution, can assign the SFAs to the real regions of the casks in less than 0.25 s for the tested instances. (iii) The local search algorithm A2[π(n)], based on 1000 initial solutions generated by A1[π(n)] by slight mutations of the π(n) sequences, is capable of reducing the initial values of the standard deviation σ(Q) 20 times compared to the solutions provided by procedure DP[π(n)] proposed in [28], in an average CPU time equal to 4.24 min.
In other words, metaheuristic MS X 0 , π(n) (made up MILP-1 plus A1[π(n)] plus A2[π(n)]) is capable of obtaining, on average, a solution with a Pearson's coefficient of variation lower than 0.75% in less than 260s CPU (1000 iterations). In the slowest approximation, MS X 0 , π(n) obtains a solution with a Pearson's coefficient of variation lower than 2.4%, starting from around 15.5% in little more than 12 min CPU (1000 iterations).
Future working lines taking advantage of these results include the following: (1) the formulation and exploitation of other models based on MILP for the regionalized cask loading according to several optimization criteria, and (2) the extension of the above models and procedures to the multi-stage load case. Acknowledgments: This work has been subsidized by the Ministry of Science, Innovation and Universities of the Government of Spain with the OPTHEUS project (ref. PGC2018-095080-B-I00), including Funds for European regional development, and the UPC project "Tool development for the simulation and optimization of the emptying of spent fuel pools in nuclear power plants" financed by ENDESA Generation. The authors also thank the AsociaciÓn Nuclear AscÓ-VandellÓs II (ANAV) the support provided and, especially, to Jordi Estrampes, head of Reactor Engineering and Nuclear Safeguards at NPP Ascó, for his tips and clarifications.