1. Introduction
The system of linear, constant real coefficient ordinary differential equations, referred to as the linear dynamic system in this article, appears in various disciplines. The linear dynamic systems can describe reactor depletion and radioactive decay dynamics and many related literature results concerning the radioactive decay linear dynamic system are available spanning over a century [1-8].
Many nuclear engineering and science applications require numerous iterations of radioactive decay calculations. An example requiring numerous iterative radioactive decay calculations is molten salt reactors (MSR) [9,10] where the reactor discharges a part of the molten salt inventory to the used fuel processor (UFP) for chemical treatment. The UFP treats the discharged used fuel to separate the fission products and feeds the processed fuel back to the reactor along with the exogenous fuel feed to make up the discharged salt from the reactor. The performance of the MSR system requiring online used fuel processing is tightly coupled with the performance of the UFP as the daily make-up fuel feed comes from the used fuel process operation. The evaluation of the MSR system performance requires numerous iterations of reactor depletion and used fuel process simulations. Another example is a pyroprocessing facility where radioactive materials are accumulated, divided, moved, and decayed daily [11-14]. The evaluation of the pyroprocessing facility via modeling and simulation also involves tens of thousands of radioactive material decay calculations. In these applications, computationally efficient radioactive decay calculations are crucial to avoid a prolonged overall calculation time. The radioactive decay calculation involves matrix exponential that can be addressed with various numerical and mathematical techniques. A review on this topic with extensive references is available in [15]. The foci have been given to find accurate, robust, and efficient algorithms for matrix exponential.
The intent of this article is not in contributing to the prescribed aspects. Rather, specifically having the radioactive decay application in mind, the focus is on identifying the structural information of the graph induced from the radioactive decay dynamics in order to fragment the problem into the set of smaller size subproblems. The computational gain depends on the degree of the problem fragmentation and the radioactive decay application exhibits a high degree of fragmentation allowing a multiple order of magnitude computational improvement. To do so, section 2 gives relevant preliminaries regarding the linear dynamic systems and graph theories. Section 3 presents two graph-based algorithms for solving the linear dynamic system. An application to the radioactive decay calculation is detailed in Section 4. Section 5 gives concluding remarks.
2. Preliminaries
Consider the system of linear, constant real coefficient ordinary differential equations given below.
Denote the set of real numbers as R. With matrix and vector notations, a compact representation of the above equation set can be written as below.
where
The constraint is that the solution X(t) satisfies the initial state X0, i.e., X(0) = X0. Then, the solution is expressed with matrix multiplication of the initial state X0 and the matrix exponential eAt as below.
In this article, we refer matrix A as the rate matrix and eAt as the transition matrix.
A directed graph G is an ordered pair (V, E) where V is the set of vertices and E⊆V×V is the set of directed edges. A directed graph G = (V, E) is induced from rate matrix A where V is the set of vertices representing the index set of equation (1) and E is the set of edges reflecting the influence structure of rate matrix A. That is,
The set {Vi} is called the partition of a set V if
A path of length k in a directed graph G = (V, E) is the sequence of vertex {v0, v1, …, vk} such that (vi–1, vi)∈E, for all i = 1, …, k. With this, we call that u∈V is reachable from v∈V if there exists a positive integer k such that {v0, v1, …, vk} forms a path of length k in G = (V, E) with v = v0, and u = vk. Let’s define the reverse edge set of E as below.
We call that a directed graph G = (V, E) is weakly connected, for every pair of vertices u, v∈V, if u is reachable from v in G' = (V, E∪E–1). A graph Gs = (Vs, Es) is called a subgraph of G = (V, E) if Vs⊆V and Es⊆E. A graph Gs = (Vs, Es) is called a strict subgraph of G = (V, E) if Vs⊂V or Es⊂E. A subgraph Ga = (Va, Ea) of a directed graph G = (V, E) is called a weakly connected component of a directed graph G if Ga is weakly connected and there is no subgraph Gb of G such that Gb is weakly connected and Ga is a strict subgraph of Gb. With this, we define the weakly connected components of a directed graph G = (V, E) denoted by {Gc,i := (Vc,i, Ec,i)} such that for all i, Gc,i is a weakly connected component and {Vc,i} and {Ec,i} are the partition of V and E, respectively.
Below gives an example illustrating the definitions presented above.
Example 1. Consider the set of linear differential equations below.
Rate matrix A in (3) can be constructed from the linear dynamic system in (5) as below.
Digraph G = (V, E) induced from rate matrix A is given in (6); Fig. 1 gives its graphical representation.
There are two weakly connected components: Gc,1 = ({1,2,3}, {(2,1),(3,2),(3,3)}) and Gc,2 =({4}, Ø).
The induced digraph facilitates an intuitive examination of the influence structure of the involved variables. Typically, the Depth First Search (DFS) over the undirected graph is used for finding the weakly connected components with O(|V|+|E|) complexity. A textbook treatise of an algorithm for finding the weakly connected components is available in [16].
3. Subdynamics Fragmentation
Denote the initial state imposed on the system given in Equation (1) by
Consider the following index set collecting the index of the nonzero elements in Equation (7). We call this the initially activated vertices.
The reach of i∈Î in the induced graph G = (V, E) can be defined as below.
The union of Ri gives the union of the reachable vertices from the initially activated vertices, denoted by RÎ. That is,
Then, we can construct the subgraph induced with RÎ, , where
The subgraph of the induced graph G = (V, E) can be constructed with an application of DFS with O(|V|+|E|) complexity.
The set of unreachable vertices from Î, denoted by , is the complement of RÎ. That is, . By grouping the reachable and unreachable vertices, Equation (2) can be rewritten as below.
where is the collection of the states corresponding to the reach of the initially activated vertices. On the other hand, represents the remaining vertices that are not in the reach of the initially activated vertices. Note that
and
As is not reachable from RÎ, we have . Then, Equation (9) becomes
Thus, the unreachable part of the linear dynamic system becomes
With , Equation (8) is further simplified to
The solution of Equation (10) is
Given the subgraph , building the rate matrix takes O(|RÎ|2). The matrix exponential can be calculated with O(|RÎ|3) [15]. The multiplication, , can be done in O(|RÎ|2.3737) with the algorithm given in [17]. Thus, overall complexity is governed by the matrix exponential computation step and O(|RÎ|3).
3.1 Algorithm with Weakly Connected Component (WCC)
Evaluating the activated part of the system with graph reachability analysis before deriving the solution provides an opportunity to reduce the size of the problem as in Equation (10). In this subsection, for the problem formulated in Equation (10), we attempt to divide the problem further into smaller subproblems with the application of the weakly connected components graph analysis described in Section 2, which can be performed with complexity. In a rough sense, two vertices of a directed graph, G, are weakly connected if a path can be found from one to the other in the undirected graph created by ignoring the direction of the edges in the directed graph, G.
First, denote the ith weakly connected component of by Gc,i := (Vc,i, Ei) where Vc,i⊆RÎ and Ec,i is the edges defined over Vc,i. That is,
Suppose that there are M weakly connected components of GÎ. By rearranging the index, one can get the following block diagonal rate matrix where rate matrix Ac,i is constructed from the connected component Gc,i := (Vc,i, Ei). We note that the value of off-diagonals elements representing the connection between the weakly connected components is zero. That is,
With the above rate matrix, the corresponding linear dynamic system can be described as below.
Then, for i = 1, …, M,
Rather than solving the linear system with Equations (4) or (11) in a monolithic manner, the solution can be derived for each diagonal block induced from the weakly connected components.
Note that synthesizing rate matrix Ac,i takes O(|Gi|2) and building all diagonal block matrices takes . As , one can see that
Equation (13) implies that the overall complexity of fragmenting the whole dynamics into the weakly connected component based subdynamics is not worse than building the rate matrix from Also, the overall complexity of building matrix exponentials of all subdynamics is , which is not worse than O(|RÎ|3), the complexity of computing matrix exponential . That is,
Other miscellaneous computation steps such as matrix multiplications and combining subdynamics solution can be performed with complexities less than O(|RÎ|3). Thus, the proposed solution strategy with the weakly connected component algorithm takes O(|RÎ|3). For brevity, we shall refer to the algorithm described in this section as WCC.
3.2 Algorithm with Minimum Starting Vertices (MSV)
The section assumes that is acyclic. The minimum starting vertex set of is the collection of vertices that do not have incoming edges. That is,
As is acyclic, all vertices in RÎ should be reachable from the minimum starting vertex set. The reachable vertices from the initially activated vertices Î can be expressed with as below.
For , one can construct the disjoint set of the reachable vertices as below.
Then, becomes a partition of RÎ. That is,
For A⊆{1, …, N}, the projection map PA : RN → RN, is defined as:
where
Denote the complement of RÎ by . As the non-zero initial condition index Î and are disjoint, we have
Note that, with any partition {V1, …, VM} of RÎ, becomes the partition of the index set {1, …, N}. The initial condition can then be expressed with the summation of the projected initial values with the partition of RÎ. That is,
Note that X(t) can be expressed as the sum of respective state evolution with the projected initial state X0 induced from any partition {Vi} of the index set V. That is,
As is a partition of RÎ, and become the partition of the index set V. Then, we have
As , we have
Let us focus on the following dynamics below.
The solution of the dynamics given in Equation (15) is
The set of unreachable vertices of G from , denoted by , is the complement of RÎ, that is, . Note that . Thus, we have and all components of X(0) belonging to become zero with the projection map . Grouping the reachable and unreachable vertices and rewriting Equation (15) result in
With , one can see that
Then, the solution of the unreachable part of the dynamics becomes for all t ≥ 0,
The dynamics of the reachable part is
The solution of the reachable part dynamics becomes
The above derivation implies that the state evolution activated with the initial state X0 can be expressed by combining subdynamics solutions described in Equation (16) that are governed by the reachable vertices from a given starting vertex and an appropriately partitioned initial value. For brevity, we shall refer to the algorithm described in this section as MSV. Below gives an example that illustrates the partition concepts described in this section.
Example 2. Consider the digraph G = (V, E) in Fig. 2. The shaded vertices in the figure are the starting vertices of the graph G. The reachable vertices from each starting vertex are R1 and R5. The partition of the vertex set V is {{1,2,3,4,6},{5}}. The corresponding initial value partition are and .
4. Application to Decay Chain Graph
The set of differential equations representing radioactive decay with no external neutron flux is shown in Equation (17). For k = 1, …, N
In Equation (17), xk(t) represents the amount of the kth nuclide at decayed to time t and ẋk is its time derivative. λk is the decay constants for the kth nuclide and fk,m is the moles of the kth nuclide created from the decay of one mole of the mth nuclide, i.e., branching fraction. Then, the rate matrix can be constructed as below:
Rate matrix A is constructed with Idaho National Laboratory Mass Tracking System Database built from ENDF/ B-VII database [18]. Decay digraph G = (V, E) that is induced from rate matrix A has the following properties:
-
|V| = N = 3631 vertices
-
|E| = 2367 transitions
-
2063 weakly connected components
-
Largest weakly connected component: 1299 vertices (mostly actinides)
-
Smallest component size: 1
-
Average component size: 1.76
-
Acyclic
Fig. 3 below shows decay digraph G = (V, E) clustering weakly connected components. Fig. 4 gives the number of reachable vertices from each vertex. Contrary to the size of the largest weakly connected component (1299 vertices), the size of the largest reachable vertex set is substantially small (36 vertices). The transition structure of the graph allures that an algorithm building the subdynamics with the minimum starting vertices is likely to be efficient. Fig. 5 shows the largest decay chain directed graph with 36 vertices (isotopes). The diagram includes neutron as a vertex as well.
MATLAB calculations were conducted to compare the performance of the algorithms described in this article.
The monolithic calculation is performed with expm MATLAB function where a scaling and squaring algorithm with a Pade approximation [15,19,20] is used. WCC and MSV algorithms also use the same expm function for their subsystem transition matrix computation.
For the above numerical calculations, the following resources are used:
Table 1 and Fig. 6 give the average computational time of the calculation methods described in this article for 100 trials. Given the number of nonzero initial isotopes to be simulated, selection of isotopes was randomized and given the value of 1. The values in the table are relative to the monolithic calculation time with 10 random nonzero initial isotopes. One can observe that the monolithic calculation ignoring the reachable vertices shows a consistent computation time. Once the reachable part of graph is considered, a significant reduction of computational time is observed for all algorithms. All three methods accounting the reachable part of the graph show a gradual increase in computational time as the number of reachable vertices increases. MSV algorithm shows the least steep increase of computational time followed by WCC algorithm and expm with reachability consideration. The superiority of expm function with reachability consideration over WCC and MSV algorithms for the smaller size initially activated vertexes may come from the implementation inefficiency of WCC and MSV algorithms where multiple graph searches were done due to the adaptation of MATLAB generic functions, which can be consolidated with a customized graph search algorithm.
5. Concluding Remarks
An attempt to view the system of linear, constant real coefficient ordinary differential equations via a graph structure allows various graph algorithm treatments to enable efficient computation of the solution. Connected component analysis, reachability analysis, and the synthesis of the minimal starting vertex set were used to break the problem into the smaller size subproblems. From a theoretical perspective, the worst-case computational complexity does not exhibit an improvement. However, empirical results demonstrated with radioactive decay application show a significant saving of computational efforts. Furthermore, we note that the subproblems are independent and allow parallel computations.