Jean-Marc Auger
email: mauger@laposte.net

Sylvain Martin

Frédéric Gruy

Inverse populations balance problems: heuristics for aggregation/breakage of solid clusters in liquid metals

Keywords: populations balance, aggregation, breakage, reaction networks, inverse methods, kernels, heuristics

Populations balance models have been applied to the simultaneous aggregation and breakage of solid clusters in agitated liquid metal. Challenges for these systems lie in kinetics constants determination (kernels) from experimental clusters concentrations measurements. However, due to difficulties inherent to molten metals experimentation (severely limiting available data points and replication), and lack of literature information concerning the physical phenomena involved, more advanced inverse methods were needed. Taking into account physical characteristics, theoretical approaches using reactions networks proved the existence of at least one stable positive equilibrium state. This result allowed the construction of two distinct fitting algorithms, aimed at solving the corresponding inverse problem (kernels determination from experimental data). The first of these heuristic methods accurately identifies kernels parameters from perfect steady state data, while the second, based on transitory states, reliably leads to correct results with as few data points as 2, and measurements errors as high as 5%.

Introduction

Populations balances are a common class of models with a wide range of applications, from biology and population ecology, to chemical engineering [START_REF] Ramkrishna | Population Balances: Theory and Applications to Particulate Systems in Engineering[END_REF] . This paper will address specifically the use of populations balance models in the context of particulate clusters aggregation and breakage in agitated liquid media.

Generally speaking, populations balance problems are characterized by three main elements:

-a set of boundary conditions (initial state, eventual entrance or exit fluxes for clusters) -a reaction network (describing every possible interaction between clusters) -a set of kernels (linking reactions kinetics to physical data and parameters) -in our case, the aggregation kernel will describe the mechanisms by which clusters colliding through liquid drag merge into larger clusters, while the breakage kernel will address clusters separation under liquid turbulence While literature shows numerous examples of kernels models built upon theoretical considerations, experimental validations are notoriously less frequent and much more complicated to implement. Determination of experimental kernels characteristics (kernel model type, and the corresponding physical parameters) is thus a difficult but critical task, needed to advance comprehension in this domain. Usual methods involve optical in-situ measurements of clusters concentrations during aggregation experiments, with the resulting concentrations evolutions with time being matched to those provided by populations balance simulations. While this approach is suitable for ordinary transparent media (such as water and most common solvents), the opacity of other fluids, such as liquid metals, makes it highly impractical, when at all possible. This obstacle is compounded by the uneasy transposition of results obtained in simple cases to more difficult ones, owing to:

-considerable quantitative changes in fluid properties (viscosity -ranging from 2x10 -4 for Rb to 5.5x10 -3 Pa.s for Co-, volumetric mass…) -qualitative modifications of involved physics: wettability (and thus surface tension), for instance, has differing meaning and mechanisms in water and liquid steel (due to the different bonds -hydrogen or metallic -present in these two media) General investigations of the inverse problem linked to population balances (e.g. the determination of aggregation or breakage kernels from clusters concentrations data) revealed it as ill-posed, adding theoretical difficulties to the experimental ones. Following works from Narsimhan et al. [START_REF] Narsimhan | Analysis of drop size distributions in lean liquid-liquid dispersions[END_REF] , Sathyagal et al. [START_REF] Sathyagal | Solution of Inverse Problems in Population Balances-II. Particle Break-Up[END_REF] , and Wright et al. [START_REF] Wright | Solutions of Inverse Problems in Population Balances-I. Aggregation Kinetics[END_REF] , many studies managed to devise successful protocols for kernels determination under the hypothesis of scaling behavior. Other regularization approaches involved approximation of kernel expressions with parametric functions such as Hermite splines [START_REF] Dürr | Parameter Estimation and Method of Moments for Multi Dimensional Population Balance Equations with Application to Vaccine Production Processes[END_REF] , or by using the method of moments 6 , [START_REF] Singh | Forward and Inverse Problems in Population Balances[END_REF] . All these procedures typically require access to a large amount of experimental data, either to validate their hypotheses, or to carry out their calculations. Sensitivity to experimental noise has only been investigated in some rare studies, such as [START_REF] Vikhansky | Droplets Population Balance in a Rotating Disc Contactor: An Inverse Problem Approach[END_REF] or [START_REF] Chakraborty | Inverse Problems in Population Balances. Determination of Aggregation Kernel by Weighted Residuals[END_REF] (limited to elementary kernels forms: sum, product, and constant). The corresponding observations seem to imply that most solving methodologies for inverse population balances additionally suffer from ill-conditioning.

In this work, we will provide tools and algorithms to determine kernels characteristics (general kernels models, then orders of magnitude for their parameters), in the context of realistic laboratory model experiments, and under the most commonly encountered conditions in liquid metals:

-non repeatability of initial conditions; while some experimental conditions such as temperature or media hydrodynamics can easily be imposed, controlled introduction of particles in liquid metals has been shown by decades of research to be a currently impossible task ( 10 , [START_REF] Agarwala | Fabrication of Aluminium Base Composite by Foundry Technique[END_REF] , [START_REF] Pai | Fabrication of Aluminium-Alumina (Magnesia) Particulate Composites in Foundries Using Magnesium Additions to the Melts[END_REF] , [START_REF] Auger | Wettability-Enabling Coating on Oxide Particles through Controlled Milling[END_REF] , [START_REF] Pai | Role of Magnesium in Cast Aluminium Alloy Matrix Composites[END_REF] , [START_REF] Auger | Preparation of Alumina Particle Suspension in Liquid Tin Using a Pre-Coating Process[END_REF] , [START_REF] Schultz | Microstructure and Hardness of Al2O3 Nanoparticle Reinforced Al-Mg Composites Fabricated by Reactive Wetting and Stir Mixing[END_REF] providing examples of such attempts). The initial clusters concentrations distribution can thus be measured, but not set.

-limited number of experimental points; this restriction follows from the use of samples for measurements, typically through metallographic analysis. As sampling disturbs the media, it must thus be kept to a minimum (up to a maximum 3 or 4 samples per experiment).

-concentrations measurements errors: metallographic analyses are carried out on 2-dimensional cross-sections from sampled ingots, while clusters are 3-dimensional structures, leading to unavoidable errors. Furthermore, the method is destructive, forbidding measurements duplication to lessen uncertainty through statistical methods. X-ray tomography could provide similar data but seems more appropriate to study individual cluster 3D characteristics (as shown in [START_REF] Li | Morphology of Nonmetallic-Inclusion Clusters Observed in Molten Metal by X-Ray Micro-Computed Tomography (CT)[END_REF] ) than global concentration distributions.

The second part of this paper will detail the formalism and hypotheses used for populations balance problems under the aforementioned conditions, then discuss the existence and nature of potential steady states. The third part will concern itself with new heuristic algorithms to determine kernels from realistic experimental data, and their limitations.

Theoretical analysis of the population balance problem

Chemical kinetics formalism will be used for both aggregation and breakage reactions, even though these phenomena are physical in nature for the cases considered.

Model hypotheses:

H1 Quantization of clusters size Existence of a minimal cluster size is postulated, and all clusters are considered aggregates of these elementary particles. This hypothesis is both realistic (as monocrystalline oxide can be considered discrete, and ultimately is at atomic scale) and not (aggregation phenomena and properties for atoms differ widely from those of macroscopic clusters, leading to transitions in models long before reaching atomic scale). However, experimental use of nonreactive (thus unbreakable) elementary particles of calibrated size coincides well with this modeling.

H2 Continuity of clusters concentrations

While classically used, this approach is highly unrealistic, as single clusters can only react in full, and concentrations values should thus be discrete. However, literature generally postulates that long term stochastic variations are similar to the corresponding deterministic steady state, even if formal proofs have only been given under very restrictive conditions [START_REF] Anderson | Product-Form Stationary Distributions for Deficiency Zero Chemical Reaction Networks[END_REF] .

H3 Closed system / elementary particles conservation While highly unrealistic for most practical cases, this hypothesis can be justified in the context of model experiments where unitary particles are introduced once in the medium and not removed. Results from these simplified experiments could be transposed to more complex situations. One example for such adaptation would be to take clusters entrance into account via additional nucleation mechanisms. From this hypothesis, with V the total volume considered, and n the total number of elementary particles, conservation Eq. 1 can be deduced.

∑ i= 1 i=∞ i⋅[ i]= ∑ i= 1 i=N ≤n i⋅[ i]= n V (1)
Following from hypotheses H1, H3, and the necessary finite volume of experimental devices, one can deduce the existence of a maximum cluster size N=n (corresponding to the totality of elementary particles aggregated as a single cluster). An interesting consequence is that simulations involving finite size kernels are sufficient to completely describe the aggregation/breakage problem as formulated in Eq. 2.

∀ (i,j )∈ℕ 2 , i+j>N, { k i,j =0 k i,j =0 } (2)
Using the notations from symbols table, the most general expression for quantized aggregation/breakage population balance problems is thus given by Eq. 3.

[

i]= ∑ j=1 j= ⌊ i 2 ⌋ ϕ j,i-j -∑ j=1 j=N -i ( 1 +δ i,j ) ⋅ϕ i,j (3)
H4 Aggregation and breakage phenomena as binary reversible reactions (see Eq. 4).

Reversibility accounts for some unusual aspects of this work, as pure aggregation or breakage have been studied

separately in several instances ( 19 , 20 , 21 , or [START_REF] Murillo | A Mean Field Model for Brownian and Turbulent Coagulation of Polydispersed Aerosols[END_REF] , for example), but much more scarcely together. This hypothesis implies the existence of at least one reversible aggregation reaction in the system, but does not assume reversibility for each. R i,j : i + j ⇔(i+j ) , ∀ (i,j )∈ℕ 2 (4) H5 First order kinetics for both aggregation and breakage reactions (e.g. use of mass action kinetics).

This hypothesis is commonly used for diluted suspensions, and expressed as Eq. 5, with ki,j and k * i,j the kinetics coefficients corresponding to aggregation and breakage kernels respectively. H5 also implies that kinetics coefficients are independent from clusters concentrations.

ϕ i,j =k i,j ⋅[ i ] ⋅[ j ]-k i,j * ⋅[ i+j ] , ∀ (i,j ) , i+j N ⩽ (5)

Existence of equilibrium states

Theorems developed for chemical reactions networks can be used to evaluate the existence and uniqueness of equilibrium states for coupled chemical reactions systems. The peculiarity of this approach resides in its independence from kinetics coefficients and initial conditions, being based on the structure of reactions networks alone. While no general method was discovered to determine networks equilibrium properties, various theorems have been devised to conclude under more restrictive conditions (summarized for instance in [START_REF] Ji | Uniqueness of Equilibria for Complex Chemical Reaction Networks[END_REF] or [START_REF] Joshi | A Survey of Methods for Deciding Whether a Reaction Network Is Multistationary[END_REF] ). Examples include deficiency zero and deficiency one theorems [START_REF] Feinberg | Chemical Reaction Network Structure and the Stability of Complex Isothermal Reactors-I. The Deficiency Zero and Deficiency One Theorems[END_REF] , higher deficiency algorithm [START_REF] Ellison | The Advanced Deficiency Algorithm and Its Applications to Mechanism Discrimination[END_REF] , or injectivity [START_REF] Craciun | Multiple Equilibria in Complex Chemical Reaction Networks: I. The Injectivity Property[END_REF] and concordance [START_REF] Shinar | Concordant Chemical Reaction Networks[END_REF] tests.

In the context of chemical reactions networks, species refer to pure chemical species (in our case, a species includes all clusters of a given size), while complexes correspond to sets of species appearing on at least one side of a reaction (they are thus every single clusters sizes and every couples of 2 clusters sizes).

Reaction graphs (see [START_REF] Craciun | Multiple Equilibria in Complex Chemical Reaction Networks: II. The Species-Reaction Graph[END_REF] for other possible representations) corresponding to the population balance problem with the aforementioned hypotheses are summarized in Figure 1. General properties and concepts pertaining to these networks will be detailed in the following sections:

Reversibility

Weak reversibility defines the existence, for each linked complexes couple (C1, C2), of oriented paths linking C1 to C2, and C2 to C1. As population balance reaction graphs present no cycle, the network is either fully reversible (when each effective reaction is reversible), or not reversible at all (when at least one reaction is not reversible).

Weak reversibility is equivalent to complete reversibility for the population balance networks studied.

Number of linkage classes l

A linkage class is defined as a set of complexes that are linked together through reactions. In the case of population balance problems, there are (N-1) linkage classes (as expressed in Eq. 6), each one corresponding to all reversible reactions leading to the formation of each clusters size ranging from 1 to (N-1).

l ( N ) =N -1 (6)

Number of complexes m

As visible in Figure 1, the reaction network contains (N-1) single species complexes, and k binary complexes per linkage class linked to k-size single clusters complex, summed in Eq. 7.

m ( N )=( N -1)+ ∑ k= 2 k=N ⌊ k 2 ⌋ =( N -1)+ ⌊ N 2 ⌋ ⋅ ( N -⌊ N 2 ⌋) (7)

Stoichiometric sub-space dimension s

The stoichiometric sub-space is defined as the vector space spanned by all chemical reactions belonging to the network. In the case of population balance problems, the easiest way to determine its dimension is to exhibit a (mathematical) basis for this space.

It can be noticed that each reaction of the network defined in Eq. 4 can be decomposed in a sum of reactions involving the aggregation of single elementary particle to existing clusters, as expressed in Eq. 8. [START_REF] Vikhansky | Droplets Population Balance in a Rotating Disc Contactor: An Inverse Problem Approach[END_REF] These (N-1) reactions are independent, and thus form a basis for the stoichiometric sub-space, leading to the conclusion (Eq. 9) that its dimension equals (N-1):

R i,j = ∑ l=j k=i+j-1 R 1 ,l -∑ k=1 k=i R 1 ,k
s ( N ) =N -1 (9)

Stoichiometric compatibility class

A stoichiometric compatibility class is defined as the set of clusters concentrations distributions that can be obtained from a given initial state (initial concentrations distribution) through any combination of the network reactions. Due to the H3 hypothesis of total elementary particles conservation, the stoichiometric compatibility class associated to a starting concentrations distribution with a total concentration of unitary particles n/V, is a subset of concentrations distributions with (obviously) the same concentration of elementary particles of n/V.

From the following observations:

-the set of concentrations distributions with n/V elementary particles concentration is of dimension (N-1) (N clusters sizes concentrations which are linked by Eq. 1)

-any stoichiometric compatibility class of the population balance problem shows the same dimension as the stoichiometric sub-space, which was determined in the preceding paragraph as (N-1)

It can be concluded about the identity between the stoichiometric compatibility class associated with any clusters concentrations distribution totaling n/V elementary particles per volume unit, and clusters concentrations distributions with n/V total elementary particles concentration.

The stoichiometric compatibility class of n/V elementary particles concentration can be represented, in the vector space of clusters concentrations distributions, by a simplex bounded by N vertices (each corresponding to a concentrations distribution with a unique non-zero clusters size concentration).

Deficiency

Deficiency is an important function of reaction networks, as its values can help determine the existence and uniqueness of equilibrium states. Its expression (using definitions and notations detailed previously) is given as Eq.

10.

δ= mls (10) This allows the determination of deficiency as a function of maximum cluster size N (Eq. 11). [START_REF] Agarwala | Fabrication of Aluminium Base Composite by Foundry Technique[END_REF] Deficiency will be evaluated for populations balance networks as a function of maximum cluster size N, and its consequences in terms of equilibrium states will be detailed, using deficiency zero and one theorems (as reminded below, from [START_REF] Feinberg | Chemical Reaction Network Structure and the Stability of Complex Isothermal Reactors-I. The Deficiency Zero and Deficiency One Theorems[END_REF] ) for simple cases, and properties of embedded networks when N value increases:

δ ( N )= ⌊ N 2 ⌋ ⋅ ( N -⌊ N 2 ⌋) -N+1

Deficiency Zero Theorem

If a mass action system is applied to a weakly reversible and deficiency zero reaction network, then, independently from kinetics constants, the system has exactly one equilibrium concentration in each positive stoichiometric compatibility class.

Deficiency One Theorem

If a weakly reversible reaction network shows the three following properties:

-deficiency for each of its linkage classes is less than 2 -deficiency of the whole network equals the sum of its linkage classes deficiencies -each of its linkage class contains only one terminal strongly linked component Then every mass action system permitted by the network admits a positive equilibrium.

N=2 case

Deficiency equals zero. If the corresponding network is weakly reversible (e.g. the only reaction in the network is reversible -all kinetics constants are strictly positive), according to zero deficiency theorem, the reaction network for these population balance problems has then exactly one stable positive equilibrium state.

N=3 case

Deficiency equals zero. If the corresponding network is weakly reversible (e.g. the 2 reactions in the network are reversible -all kinetics constants are strictly positive), according to zero deficiency theorem, the reaction network for these population balance problems has then exactly one stable positive equilibrium state.

N=4 case

Deficiency equals one. According to deficiency one theorem [START_REF] Feinberg | Chemical Reaction Network Structure and the Stability of Complex Isothermal Reactors-I. The Deficiency Zero and Deficiency One Theorems[END_REF] , deficiency of each linkage class in the network should be evaluated in order to prove existence and/or uniqueness of equilibrium states. From the data summarized in Table 1, each linkage class has deficiency equal to zero, while the sum of all deficiencies is one, which precludes the use of deficiency one theorem. However, the deficiency one algorithm [START_REF] Feinberg | Chemical Reaction Network Structure and the Stability of Complex Isothermal Reactors-II. Multiple Steady States for Networks of Deficiency One[END_REF] , used through the CRNToolbox implementation (developed by Martin Feinberg Group for Chemical Reaction Network Theory at The Ohio State University), can successfully determine that such networks admit either exactly one stable positive equilibrium state, or none, depending on the reversibility of individual reactions, as summarized in Table 2.

N>4 case

For N greater than 4, the corresponding deficiency is greater than 1 and increasing with N. Simple theorems and classical approaches (such as 31 , 32 , 33 , 34 , 35 , 36 , 23 , 37 , 38 , or [START_REF] Craciun | Multiple Equilibria in Complex Chemical Reaction Networks: II. The Species-Reaction Graph[END_REF] ) cannot be applied anymore. According to [START_REF] Joshi | Atoms of Multistationarity in Chemical Reaction Networks[END_REF] , it is however possible to conclude on the existence of steady states using the following property: a reaction network admits at least the same number of stable positive equilibrium states as one of its embedded networks with which it shares the exact same stoichiometric sub-space. In the context of reactions networks, an embedded network is constructed from an initial network by the removal of any number of reactions and/or chemical species.

An example of such an embedded network applied to population balance is presented in Figure 2, and obtained through the removal of every reaction not corresponding to the aggregation of single elementary particles to clusters.

As presented before, dependency of characteristic functions for the embedded network can be expressed in Eq. 12,

Eq. 13, and Eq. 14 with N as a variable.

m ( N )=2⋅( N -1) (12) l ( N ) =N -1 (13) s ( N ) =N -1 (14)
Deficiency can then be expressed in Eq. 15 as a function of the maximum cluster size.

δ ( N ) =m ( N )-l ( N )-s ( N )=0 (15)
By construction, the embedded network stoichiometric sub-space is included into the original network stoichiometric sub-space. Furthermore, these sub-spaces share the same dimension (Eq. 9 and Eq. 14), and are subsequently one and the same. Consequently, under the condition of embedded network weak reversibility (equivalent to the reversibility of the aggregation reactions of single elementary particles to clusters, with this choice of embedded network), zero deficiency theorem applies to the embedded network, which then admits exactly one stable positive equilibrium state. From 33 , the reaction network for population balance problems has then at least one stable positive equilibrium state.

It should be noted that such reasoning can be applied to any embedded network comprised of exactly one (effectively) reversible reaction per linkage class of the original network.

As a conclusion for N>4 population balance problems, the existence of at least one stable positive equilibrium state is guaranteed, provided (sufficient condition) there is at least one effectively reversible aggregation/breakage reaction leading to the formation of each clusters size ranging from 2 to N.

Multiplicity of equilibrium states

An example will show the potential for multiple equilibrium states in a given population balance problem. Kernel data are summarized in Table 3, and the corresponding boundary conditions are following: maximum cluster size of N=5, total elementary particles concentration of 1 particle per volume unit. In this case, two equilibrium states shown in Table 4 are allowed, depending on initial conditions (e.g. starting clusters concentrations distribution) of the problem.

In order to assess the dependency of equilibrium states multiplicity with boundary conditions (and especially total elementary particles concentration), population balance simulations were carried out (until steady state) using this same kernel and multiple starting points inside various stoichiometric compatibility classes. As demonstrated before, in the case of population balance problems verifying H3 hypothesis, a compatibility class corresponds to all clusters concentrations distributions with a given elementary particles total concentration.

As mentioned in 1.2.5, each stoichiometric compatibility class can be seen as a simplex in the clusters concentrations distributions vector space. As a way to maximize the regularity of simulation starting points, they were chosen as evenly spaced points on each edge of the simplex.

Edges are defined as linking two distinct vertices of the simplex, which in turn correspond to clusters concentrations distributions with only one non-zero concentration. Concentrations distributions belonging to an edge can thus be expressed as Eq. 16.

(i,j ) ∈(1

: N= 5) 2 ,i<j, { ∀ l ∈[ 0,1 ] , ( [1]=0 , .. . , [i ]= n V⋅i ⋅l, . . . , [ j ]= n V⋅ j ⋅(1-l ) ,. . . , [ N ]=0
)} [START_REF] Schultz | Microstructure and Hardness of Al2O3 Nanoparticle Reinforced Al-Mg Composites Fabricated by Reactive Wetting and Stir Mixing[END_REF] Five evenly spaced concentrations distributions on each edge for every compatibility class have been used as starting points. For N=5, the simplex has 10 edges, for a total of 35 starting concentrations distributions, defined by Eq. 17.

(i,j ) ∈(1

: N= 5) 2 ,i<j, { ∀ l ∈ { 0 , 1 4 , 1 2 , 3 4 ,1 } , ( [ 1 ]=0 , .. . , [i ]= n V⋅i ⋅l, .. . , [ j ]= n V⋅ j ⋅(1-l ) ,. . . ,[ N ]=0 )} (17)
From these initial distributions, steady states were computed through population balance simulations, using tools detailed in part 3: "Kernels determination through experimental data fitting". Three consecutive time step sizes (2.10 -3 s; 2.10 -6 s; 2.10 -9 s), have been applied 10 4 times each. Final instant concentrations variations were verified to be less than 10 -10 particles per volume unit and second, confirming that steady states were actually reached.

Comparison between equilibrium states derived from starting concentrations distributions corresponding to the same stoichiometric class has been carried out through Bhattacharyya distance calculations. The upper bound of computed distances for each class (equivalent to a given total elementary particles concentration value) is represented in Figure 3. As visible on the figure, low and high values of particles concentration are characterized by a single equilibrium state, while intermediary values allow the existence of multiple (in fact two) balance states. Since the kernel used is the same for all concentrations, this last parameter alone cannot be used to determine the multiplicity of solutions for population balance problems with clusters sizes higher than 4. As a conclusion, while the reactions network structure allows the possibility of multiple equilibrium states for N greater than 4, boundary conditions (total particles number, maximum clusters size) impose the actual number of steady states.

Kernel models

Apart from the 3 classical elementary kernel forms (namely, constant, sum, and product), which have the advantage of providing analytical solutions [START_REF] Dubovskii | Mathematical Theory of Coagulation[END_REF] , but do not correspond to any realistic experimental setup, several kernels formalisms have been developed for various applications (coagulation of solid suspensions in liquids, polymerization kinetics, liquid drops coagulation in gas...). Generally speaking, aggregation kernels were far more investigated than breakage ones, hence explaining their higher mathematical complexity and number of variants.

In the limited scope of solid clusters coagulation and breakage in liquids, 6 general kernel classes (3 for aggregation, and 3 for breakage) can be derived from literature reviews ( 41 , for example). They are detailed in Table 5. Each of these mathematical expressions corresponds to several physical models (for example, Kruis-Kuster and Saffman-Turner models are both variants of the turbulence type 1 general model), the physical meanings of coefficients varying wildly with the hypotheses used in their construction.

As easily visible from the expressions in Table 5, aggregation (respectively breakage) kernel coefficients ki,j (respectively ki,j * ) can locally (e. g. for an (i,j) clusters couple) take a value of zero, only if the corresponding kernel is a constant and null. From this observation, it can be deduced that either each aggregation reaction is reversible, or none is. In combination with H4 hypothesis, it can be concluded that classical kernel expressions lead to reversible reactions exclusively, which in turn (using demonstrations from the preceding paragraphs) ensures the existence at least one stable positive equilibrium state, independently from maximum clusters size.

Nature of equilibrium states

The most general form of equilibrium is the complex balanced state, which is characterized by the lack of evolution of complexes (and consequently species) concentrations with time. A more restrictive case is the so-called detailed balanced state, in which each elementary aggregation/breakage reaction is balanced. As such, detailed balanced states are much simpler to investigate theoretically, since they can be studied through the classical tools of single reaction kinetics. We will identify some necessary conditions for the existence of detailed balanced state in the context of population balance problems.

If an equilibrium state is detailed balanced, using the chemical kinetics formalism, and under the hypothesis of reversibility for each reaction, a general expression for equilibrium constants is given by the Eq. 18 system. [START_REF] Anderson | Product-Form Stationary Distributions for Deficiency Zero Chemical Reaction Networks[END_REF] Easily verified recursive reasoning using Eq. 18 gives relationships between generic equilibrium constants, and constants corresponding to the reactions involving the aggregation of elementary particles. The resulting expressions in Eq. 19 constitute a necessary condition for detailed balanced states.

{ K ij = [(i+j )] eq [ i ] eq ⋅[ j ] eq = k ij k ij * } ( i= 1 . . . N -1 ,j= 1 . . .i,i+j ≤ N )
{ K ij = ∏ l=max ( i,j ) i+j-1 K 1 l ∏ k=1 min (i,j ) -1 K 1 k } ( i,j ) ∈(ℝ ∖ {1}) 2 (19)
A simpler incremental formulation, as shown in Eq. 20, can be used as a criterion to rule out the possibility of detailed balanced states. If this criterion is different from 1 for at least one (i,j) couple in a given kernel, balanced states can only be of the complex variety for this kernel. [START_REF] Cheng | Kinetics of Fragmentation[END_REF] Formal expressions for this criterion have been calculated for each of the nine combinations of physical kernel models detailed in the preceding paragraphs. As easily visible in Table 6, criterion equality with 1 can only be satisfied when the maximum clusters size N equals 2. This configuration is, by definition, that of the equilibrium of a single reversible reaction, which can be considered as both detailed and complex balanced. For i or j strictly higher than 1, and independently from aggregation and breakage kernels models parameters, criterion values differ from 1.

{ K 1 ,i ⋅K i+1 ,j K 1 ,i+j ⋅K i,j =1 } ( i,j ) ,i+j<N
It can thus be concluded that realistic kernels models restrict equilibrium states for population balance problems, under the aforementioned hypotheses, to complex balanced states.

As a general conclusion for this part, chemical reactions networks considerations proved the existence of at least one positive equilibrium state for each population balance problem, independently from kinetics models (kernels) involved. Balance states can thus be used as a basis for inverse methods (e.g. determination of kernels models and parameters from experimental data) in all cases, which is simpler than using transitory states. Furthermore, general kernels models types constructed from physical consideration have been identified in literature, reducing the scope of kernels identification to 9 aggregation/breakage models combinations. Finally, a criterion has been devised to identify the nature of population balances equilibrium states. Its application demonstrated that physically realistic kernels can only attain complex balanced states.

From these results, inverse methods applied to populations balance problems convey two main challenges. The first is the determination of best fitting constants for each of the 9 kernels types combinations based on a set of experimental data. The other resides in the ranking of these models against each other, in order to identify which one is the best suited to describe the aggregation/breakage phenomena studied. These problems will be addressed in the following part.

Kernels determination through experimental data fitting

Data fitting can be carried out along two distinct strategies, involving either the use of equilibrium state concentrations distributions (which existence was demonstrated in part 2), or the analysis of multiple transitory concentrations distributions. Each approach shows marked and distinctive advantages:

-Equilibrium data is fundamentally easier to obtain, as each populations balance experiment will finally reach this state. Time and duration measurements are also unnecessary since equilibrium is by definition a steady state.

-Dynamic concentrations distributions are more difficult to assess (in-situ experimental measurements are needed), and more sensitive to errors (for both concentrations and time). However, the potential to enforce any experimental time frame (duration between measurements) can be especially useful in the analysis of complex experiments involving multiple aggregation/breakage mechanisms that occur at different time scales. By choosing the right experimental time intervals, some phenomena can be made negligible, thus simplifying results interpretation.

The following paragraphs will first provide some information about simulation tools and their use to generate realistic pseudo-experimental data. Algorithms constructions corresponding to each of the aforementioned strategies will then be detailed and validated using simulated data.

Simulation tools

Between the lack of generalized theoretical tools, and minimal number of population balance problems allowing analytical solutions, numerical simulations currently constitute the only reliable method when studying both transitory and steady states. To this end, the C++ framework for population balances (PopulationBalance: [START_REF] Auger | PopulationBalance C++ Library for Population Balance Problems Simulations[END_REF] ) has been devised to provide simulation and analysis tools, with performance in terms of CPU time and memory usage as the main objective.

Tools in this framework were devised with the hypotheses mentioned in the first part of this paper in mind and are thus optimized for similar cases. Adaptation to other situations is nonetheless possible by using clever boundary conditions (Neumann and Dirichlet types), at the expense of CPU time. Hypotheses H1, H4 and H5, however, are too ingrained in the internal computing algorithms to be easily bypassed.

This framework allows the manipulation of two different types of objects: kernels and clusters concentrations distributions.

Three different kernel types can be generated: simple mathematical expressions (such as sum or product), general physical models as detailed in the previous part, and specific models based on physical parameters (Saffman-Turner, or Levitch kernels for example) were implemented.

Basic operations between kernels (corresponding to coefficient to coefficient operations) allow the construction of hybrid or even discontinuous kernels models.

Manipulation of clusters concentrations distributions follows the rules classically applied to probability distributions, with one exception: distances computations. Probability distributions distances can only be calculated under the assumption that the two distributions are normalized, a condition that is incompatible with clusters conservation.

Definitions for Total Variation Distance, Hellinger, Bhattacharyya, and Manhattan distances are thus modified by first normalizing the concerned clusters distributions, before applying the standard formulas reproduced in Eq. 21.

d TVD (( [i ] 1 ) i≤ N , ( [i] 2 ) i≤ N ) =max i≤ N ( |[ i ] 1 -[i ] 2 | ) d Manhattan (( [i ] 1 ) i ≤ N , ( [i ] 2 ) i≤N ) = ∑ i≤ N |[i ] 1 -[ i] 2 | d Hellinger (( [ i ] 1 ) i≤ N , ( [i ] 2 ) i≤N ) = √ 1-∑ i≤ N √ [i] 1 ⋅[ i] 2 d Bhattacharyya (( [i ] 1 ) i≤ N , ( [ i] 2 ) i≤ N ) =-ln ( ∑ i≤ N √ [i] 1 ⋅[i ] 2 ) (21)
Interactions between clusters concentrations distributions and kernels correspond to population balance advancement. Such time steps are applied using either an explicit Euler method, with the objective to minimize CPU time, or a second order Runge-Kutta method. Time steps used for simulations of experiments ("experimental simulations") are deliberately very small (10 -8 s) in comparison to those (10 -6 s) of simulations used in the various fitting algorithms ("analysis simulations"), to better reflect that experimental results correspond to infinitely small time increments.

Steady states are generated from a starting clusters concentrations distribution (unless noted otherwise, this initial state correspond to a maximal number of clusters, e.g. all clusters are of elementary size). Time steps (10 -8 s) are then successively applied until the Bhattacharyya distance between clusters distributions at two consecutive steps is less than a limit value. This limit is set to 10 -18 , and by analogy to Cauchy criterion, ensures that simulations are stopped close enough to the equilibrium state.

Transitory states are similarly obtained from a starting clusters concentrations distribution corresponding to a maximal number of clusters. The difference is that time steps of 10 -8 s are only applied until a fixed total duration has elapsed, to obtain an out of balance clusters concentrations distribution. Such distributions are typically generated in pairs, the first corresponding to a simulated duration from the initial state, and the second to a simulated time interval following the first state. This time interval is considered a part of transitory states data, while duration from initial distribution is not (duration between experimental measurements can easily be chosen, but precise control of the initial state is far more difficult).

Concentrations distributions thus generated simulate perfect experimental data, which is highly unrealistic, as measurements are always subject to noise. Therefore, a good fitting algorithm must not only identify the right kernel models (kernels expressions and constants values both) from perfect data, but also demonstrate robustness when applied to noisy information. Considering that concentrations are measured independently for each cluster size, three hypotheses have been used to accurately model the noise added to simulated concentrations distributions:

-noise is applied independently for each clusters size concentration -noise follows a centered normal distribution -noise standard deviation is proportional to the concentration it is applied to, with the same factor for all clusters concentrations in a given distribution From now on, a 4% error (or 4% noise) on concentrations distribution will describe a distribution for which the standard deviation of the error applied to each cluster concentration equals 4% of the undisturbed concentration value.

Considering that proper kernels models are currently unidentified for liquid metals, computation parameters (durations, constants, clusters concentrations…) used thereafter in this paper are voluntarily unrealistic, as the aim of this part is to provide generalized examples and validate methods.

Equilibrium states fitting

Using data corresponding to steady states presents several advantages, one being that equilibrium is practically easily reached (prolonged experiments always lead to equilibrium by definition), and another laying in the limitation to only one necessary measurement point. Most methods aimed at fitting theoretical models to experimental results rely on the use of some objective function, whose parameters can be divided between the experimental data set to fit and model parameters (constrained multivariate optimization). In the case of inverse populations balance problems, the data set is a measured equilibrium clusters concentrations distribution ([i])i≤N, and the model parameters consist of aggregation and breakage kernels K and K * (kernels mathematical expressions and their corresponding constants).

From Eq. 3 and 5, a possible definition for an objective function C would be given by Eq. 22.

C ( ( [

i] ) i ,K,K* ) = ∑ i= 1 i=N | [ i ] |= ∑ i= 1 i=N | ∑ j=1 j=i-1 ( k j,i -j ⋅[ j] ⋅[ i-j]-k j,j -i * ⋅[ i] ) -∑ j= 1 j=N -i ( 1 +δ i,j ) ⋅ ( k i,j ⋅[ i ] ⋅[ j ]-k i,j * ⋅[ i+j ] ) | (22)
By construction, criterion C has the following properties:

1-positive values 2-zero is reached only when the clusters concentrations distribution ([i])i≤N corresponds to an equilibrium state for kernels K and K * (distribution and kernels will be called compatible)

3-if kernels K, K * and distribution ([i])i≤N are compatible, then for any positive real number a, kernels a.K, a.K * (a.K defined as the kernel whose coefficients equals those of K multiplied by a) and distribution ([i])i≤N are also compatible (direct consequence of complex balanced state, which was demonstrated as the only equilibrium type possible in Part 2)

For a chosen equilibrium distribution ([i]eq)i≤N, criterion C can thus be used as an objective function whose minimization with regards to K and K * allows the identification of compatible kernels. Due to the general multiplicative coefficients common to all kernel models (α coefficients in Table 5), property 3 implies that the set of parameters allowing compatibility for an aggregation/breakage kernel models couple is either infinite or empty.

Consequently, minimizing algorithms cannot be applied directly to criterion C with the chosen kernel models. In order to avoid this effect, the multiplicative coefficient of aggregation kernel expressions can be locked to a single value (1 will be used), reducing the number of parameters to determine.

The algorithm used for C criterion minimization (a non-linear problem) is the Nelder-Mead heuristic method (detailed in [START_REF] Nelder | A Simplex Method for Function Minimization[END_REF] ), which is based on simplexes properties. Implementation simplicity and efficiency when applied to few parameters (2 or 3 depending on the kernels considered) make this method particularly apt to the fitting of kernels models to experimental data.

Examples of application for kernels determination: perfect experimental data

Pseudo-experimental data were generated using the following basis: Kernels:

-Turbulence type 1 aggregation, α=1 -Exponential law 2 breakage, n=0.5, α=0.25, β=3

Step size: 10 Using the previously described fitting method, optimal parameters for each of the 9 kernel models combinations have been determined, and the results summarized in Table 7. Two interesting constatations can be made:

-first, kernel constants determined by fitting the right models (Turbulence type 1 and Exponential law 2) match perfectly (less than 0.2% error) the values used for pseudo-experimental data generation.

-second, the correct kernel models couple is easily identified, as it corresponds to the lowest value for C.

Furthermore, the C criterion value reached at the end of the minimization procedure shows a marked discriminating effect, being six orders of magnitude smaller for the right models than for non-matching kernels.

Generality of these results was easily confirmed, by applying the same fitting method to pseudo-experimental data generated using various kernels models and parameters.

Examples of application for kernels determination: noisy experimental data

In order to assess the effect of noise on the fitting procedure, a slightly disturbed pseudo-experimental data set was generated from: Kernels:

-Turbulence type 1 aggregation, α=1 -Exponential law 2 breakage, n=0.5, α=0.25, β=3

Step size: 10 As shown in Table 8, the turbulence type 1 and exponential law 1 kernels combination is falsely identified as best fitting the pseudo-experimental data. Furthermore, even using the right theoretical models leads to highly incorrect constants (the minimal error is around two orders of magnitude). This effect could be linked to two different phenomena: an important sensitivity of the algorithm to entered data (ill-conditioning), or a systematic error following the use of noisy data (bias). In order to discriminate between these causes, 200 pseudo-experimental data sets were generated using the same kernels and parameters as previously mentioned, with the fitting results corresponding to the correct kernels models being summarized in table Table 9 for 25 of them. Errors on kernels parameters determined by the algorithm do not follow a systematic pattern, and seem dispersed between several orders of magnitude. This effect hints at a low condition number (ratio between relative error of computed parameters and relative error of entered data) for this approach.

As a general conclusion, the inverse method devised for kernels determination from equilibrium states is highly accurate for both kernel models discrimination and constants determination, provided that equilibrium data is noiseless. However, even small magnitude noise lead to incorrect kernels identification, and huge errors in constants determination. This sensitivity to data accuracy is a telling symptom of an ill-conditioned problem-solving algorithm.

Transitory states fitting

As mentioned before, experimental or pseudo-experimental data used to define transitory states consist in two consecutive clusters concentrations distributions, together with the time interval separating these 2 measurements.

The issue is thus to identify an aggregation/breakage kernels couple, such that there exist a population balance "trajectory" linking the two aforementioned distributions with the correct time interval. This kind of problem is generally the province of least squares methods. However, a population balance trajectory is defined not only by two kernels models and the corresponding parameters (3 or 4 in total), but also by a starting concentrations distribution (100 or more additional parameters). The considerable number of parameters involved, together with the non-linear nature of the problem, renders classical multivariate fitting methods inapplicable. The following paragraphs detail a variant version of the least squares method devised to circumvent this limitation.

Using the first measured clusters concentrations distribution ([i]initial)i as an initial state, a populations balance simulation (PBS) is carried out with kernels candidates K and K * for a time interval Δt corresponding to the duration between experimental (or pseudo-experimental) measurements (formally expressed as Eq. 23). The resulting simulated distribution ([i]simul)i is then compared to the second (pseudo-)experimental distribution ([i]final)i through their Bhattacharyya distance, which constitute the objective function T (Eq. 24) to minimize with respect to K and K * .

PBS ( K,K*,Δt ) : ( [ i ] initial ) i≤ N → ( [i ] simul ) i≤ N (23) T ( K ,K *,Δt, ( [ i ] in itia l ) i ≤ N , ( [i ] fin a l ) i ≤ N ) = d B a tta ch a ryya (( [i ] fina l ) i ≤ N , ( [ i ] sim u l ) i≤ N ) (24)
Similarly to equilibrium states fitting, the method used for T criterion minimization (also a non-linear problem) is the Nelder-Mead heuristic approach. Due to the more complex definition of objective function (which includes a complete population balance simulation), this algorithm displays longer CPU times.

This fitting method provides results like those obtained using equilibrium states (some reproducibility examples are given in Table 10 for lightly disturbed transitory distributions):

-accurate kernel models determination and parameters identification are possible, when based on exact (pseudo-)experimental data -an inability to identify the correct kernel models using even slightly noisy (1% error) (pseudo-)experimental data.

Determination of kernel parameters in these conditions suffers from a reproducibility problem, with instability of computed values ranging several orders of magnitude.

In order to understand the reasons for this demonstrated sensitivity of the fitting algorithm to imprecise data and determine potential workarounds, statistics of errors in clusters concentrations distributions had to be further investigated.

Measurements errors analysis

Models for errors in clusters concentrations distributions were selected for their agreement with experimental measurement procedures. Each clusters size concentration is measured separately, leading to independent errors that are often considered the result of multiple random phenomena. From central limit theorem, individual concentration errors can be realistically thought-out as following a normal distribution, centered on an undisturbed concentration value. Hypotheses pertaining to pseudo-experimental data generation can thus be summarized as Eq. 25.

∀ i ∈(1 ;N ) , [ i ]=[ i] +X i with X i =Ν (0 ,k⋅[ i ]) ( 25
)
Variable k is the measurement error as a proportionality factor to clusters concentration.

In order to check an eventual bias linked to concentrations errors, knowledge of the probability distribution of noisy concentrations distributions is needed. As literature does not provide any example of theoretical study for this problem, it will be qualitatively evaluated through numerical simulations. Starting from a base concentrations distribution (maximum clusters size N=100), 100000 disturbed distributions were generated by adding 5% noise to individual concentrations. Bhattacharyya distances between base and disturbed distributions were calculated, and the corresponding frequencies determined using ten evenly spaced classes, as presented in Figure 4. It should be noted that the use of different distance definitions, such as Hellinger, Manhattan, Total Variation Distance or Mahalanobis, leads to the same characteristics, whose consequences will be discussed in detail:

-the probability distribution does not follow a normal law (no symmetry around the maximum probability distance).

Statistical properties of errors in clusters concentrations distributions can be likened to those of multivariate normal distributions.

-highest probability is not reached for a distance of zero. This observation reveals the existence of a statistical bias in disturbed experimental data when considering clusters distributions distances. The actual clusters distribution is most probably situated on a sphere at some distance around the measured distribution.

-probability around a distance of zero is very low. This effect is of extreme import to the accuracy of fitting algorithms. It is highly improbable that any measured distribution closely matches the real distribution. Since the transitory states fitting method is based on distance minimization between experimental (thus noisy) points and simulated points, the algorithm as presented will ensure (statistically) a high error in kernels parameters determination.

This kind of bias is classically dealt with through the multiplication of experiments under similar conditions.

Provided enough data points, the corrected clusters concentrations distribution will be at the center of the noisy points cloud. In the case of clusters concentrations, such correction is obtained by a separate averaging of the concentration for each cluster size. This approach, however, is limited by experimental repeatability. Sufficient data accumulation is impossible in the context of liquid metal tests, due to uncontrollable initial conditions and destructive measurement methods (metallography). As a replacement, another corrective method based on statistical maximum likelihood can be used instead, as described in the following paragraph.

Maximum likelihood statistical correction

The objective of this part will be to construct a function (the likelihood function) of an experimental (noisy) data point (e.g. clusters concentrations distribution), which is statistically minimal when simulated points approach the corresponding undisturbed concentrations distribution.

Under the hypotheses already detailed in the previous part in Eq. 25, normalizing and centering of the normal law leads to the expression of Eq. 26.

[ i ]-[i ]

k⋅ [ i ] ≃ [ i ]-[i ] k⋅[i ] =Ν (0,1) (26)
By definition of Χ 2 law, summation of Eq. 26 for all values of i gives Eq. 27.

∑ i= 1 i=N ( [ i ]-[ i ] k⋅ [ i ] ) 2 Χ ≃ 2 ( N ) (27)
If N is sufficiently large (N>100), Χ 2 law can be accurately approximated by a normal law leading to Eq. 28.

∑ i= 1 i=N ( [ i ]-[i] k⋅ [ i ] ) 2 ≈ Ν ( N ,2⋅N ) (28)
Centering this normal law (Eq. 29), the probability of the expression becomes maximal at zero.

( ∑ i= 1 i=N ( [ i ]-[i] k⋅ [ i] ) 2 ) ⋅ 1 N -1≈ Ν (0,2 ) (29)
Finally leading to the maximum likelihood function expression f (Eq. 30), which follows the half-normal distribution.

f (( [ i ]) i , ([ i ]) i ,k,N )= | ( ∑ i= 1 i=N ( [ i ]-[ i ] k⋅ [ i ] ) 2 ) ⋅ 1 2⋅N -0 . 5
| ≈|Ν (0,1)| [START_REF] Feinberg | Chemical Reaction Network Structure and the Stability of Complex Isothermal Reactors-II. Multiple Steady States for Networks of Deficiency One[END_REF] f has the required properties to be used in the transitory states algorithm in place of Bhattacharyya distance (confirmed by 100000 simulations, N=100, k=5% -see Figure 5):

-f is positive, and consequently f minimal value is zero.

-f follows the half-normal law, ensuring a maximal probability when f is minimal.

Conclusion:

For a given disturbed concentrations distribution (maximum clusters size N, k% error) ( [ i] ) i , the clusters concentrations distribution ([i])i such as the function f defined by Eq. 31 is minimal, offers the most probable correspondence with the undisturbed concentrations distribution.

f (( [ i ]) i , ([ i ]) i ,k,N )= | ( ∑ i= 1 i=N ( [ i ]-[ i ] k⋅ [ i ] ) 2 ) ⋅ 1 2⋅N -0.5 | (31)

Transitory states fitting with maximum likelihood correction

Maximum likelihood correction can be used to improve kernels fitting methods without relying on measurements multiplication. The algorithm previously described for transitory states fitting is altered by substituting the objective function T with another, TC, constructed from the maximum likelihood function f (as defined in the previous paragraph). This modified algorithm is described by Eq. 32 and TC function in Eq. 33.

-Turbulence type 2 aggregation, α=1

-Exponential law 2 breakage, n=0.5, α=1, β=3

Step size: 10 Fitting results samples using correct kernels models are detailed by tests N°1 to 5 in Table 12. While still much higher (around one order of magnitude) than the original data errors (2% only), the variability of computed kernels parameters was dramatically diminished using maximum likelihood corrected algorithm. It should be noted that anomalous parameters values are still computed in some cases (α value for breakage kernel model in N°4 test, for example), but the frequency of these happenings is decreased compared to that shown by using T criterion (see Table 10). Generally speaking, multiplicative factor α tends to be more precise than exponential factor β, which in turn is more precise than exponent n.

As a conclusion, maximum likelihood correction greatly improves the accuracy of kernels parameters determination. However, as a statistics-based approach, this method is still dependent on actual measured data (two clusters concentrations distributions), that can occasionally take values too far removed from the most probable ones to be effectively corrected. In order to establish the limits of this improved algorithm (and thus determine the best experimental setup to optimize its use), sensitivity to two easily controlled experimental conditions has been investigated: the time interval between transitory states measurements, and how far removed they are from equilibrium state.

Sensitivity to time duration between (pseudo-)experimental points:

Cases N°11 to 15 in Table 12 were generated using the following characteristics (corresponding to a duration decrease between the pseudo-experimental transitory states in contrast to cases N°1 to 5): Kernels:

-Turbulence type 2 aggregation, α=1

-Exponential law 2 breakage, n=0.5, α=1, β=3

Step size: 10 between measurements for transitory states diminishes. This effect is like the increasing slope error observed if base points are closer from each other, when performing linear regression algorithms.

Sensitivity to distance from equilibrium state:

Cases N°6 to 10 in Table 12 were generated using the following characteristics (corresponding to an increase in duration from initial simulation state compared to cases N°1 to 5): Kernels:

-Turbulence type 2 aggregation, α=1

-Exponential law 2 breakage, n=0.5, α=1, β=3

Step size: 10

Conclusions

Part 2 investigated the specifics of population balance problems, notably using reaction networks theory, with a special interest toward eventual equilibrium states and their attributes.

-A set of hypotheses that closely matches model experimental setups (clusters quantization, particles conservation, concentrations continuity, binary reversible aggregations, and mass action kinetics) allowed the use of embedded reaction networks for population balance problems with maximum cluster size N. The existence of at least one stable positive equilibrium state was demonstrated, provided (sufficient condition) there is at least one effectively reversible aggregation/breakage reaction leading to the formation of each clusters size ranging from 2 to N.

-Boundary conditions (total elementary particles concentration) and kernels mathematical formulations have been shown to determine the multiplicity of equilibrium states. However, no simulation involving physical kernels expressions has been observed to generate multiple equilibrium states during this study.

-Possibility for detailed balanced equilibrium state can be assessed through a simple criterion based on kernels formulae. This criterion was constructed under the hypothesis of reversible elementary aggregation reactions.

-Using the 6 general kernels expressions found in literature for aggregation and breakage of solid clusters in liquid media, equilibrium states for realistic population balance problems were thus demonstrated to be of the complex balanced variety.

Part 3 concerned itself with the population balances inverse problem (e.g. determination of aggregation and breakage kernels from experimental data). Making use of the specifically developed PopulationBalance C++ framework for implementation and testing, two different fitting heuristics have been devised to solve this non-linear problem:

-The first algorithm, based on equilibrium states, is a simple and highly efficient method to both identify kernels types and compute their parameters, when applied to undisturbed (perfect) experimental data. However, extreme sensitivity to measurement errors makes it unusable in any realistic experimental context and exposes the inverse problem as an ill-conditioned one.

-Similar observations were made on a second algorithm constructed from least squares methods applied to transitory states. Sensitivity to measurements uncertainties evidenced a systematic bias that could be traced back to statistical errors spread in clusters concentrations measurements.

A statistically unbiased objective function (TC) was devised as a substitute to distributions distances in the algorithm. This improvement reduced error ranges for computed parameters to a single order of magnitude, and allowed accurate kernels identification using minimal experimental data (only two points) with 5% or less errors.

Study of the corrected algorithm sensitivity to experimental measurement parameters led to the following (and opposite) recommendations, in order to optimize kernels determination:

-measurement points must be separated by the greatest time interval possible -measurement points must be as far removed from equilibrium as possible

- 8 s

 8 Total elementary particles: 5/volume unit Maximum clusters size: N=200 Initial state: [1]= 5 clusters/volume unit Noise: 0% (no noise)

- 8 s

 8 Total elementary particles: 5/volume unit Maximum clusters size: N=200 Initial state: [1]= 5 clusters/volume unit Noise: 1%

- 8 s

 8 Simulated time from initial state: 0.01 s Simulated duration between transitory states: 0.025 s Total elementary particles: 5/volume unit Maximum clusters size: N=150 Initial state: [1]= 5 clusters/volume unit Noise: 2%

- 8 s

 8 Simulated time from initial state: 0.01 s Simulated duration between transitory states: 0.005 s Total elementary particles: 5/volume unit Maximum clusters size: N=150 Initial state: [1]= 5 clusters/volume unit Noise: 2% Comparison of results N°1 to 5 and N°11 to 15 shows of marked decrease in fitting accuracy when the time interval

FIGURESFigure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .

 12345 FIGURES

  The latter case is classically more stable, but configurations involving at least one concentration reaching zero can lead to unavoidable negative concentrations values on the following time step, and must thus be carefully monitored. Typical durations for complete populations balance simulations with clusters size reaching several thousand elementary particles can be estimated (depending on precise simulation

	Sources	for	this	framework,	documentation,	and	examples	of	use	are	provided	at
	https://gitlab.com/AugerJM/populationbalance.
	3.2 Pseudo-experimental data generation
	As a means to validate inverse methods algorithms, data obtained through controlled experimental parameters
	(kernel types, kernel constants, initial state, experiment duration, measurement errors…) would be needed. Due to
	the practical impossibility to carry out such experiments (experimental kernels are by definition unknown), clusters
	concentrations distributions were generated using populations balance simulations.

conditions) as a few minutes to a few hours. A time step adaptation mechanism was also added to the Euler algorithm, in order to forbid clusters concentrations from reaching negative values.

  Accuracy of kernels parameters determined by the algorithm is visibly degraded with data points nearer to the equilibrium state (comparison of cases N°1-5 to 6-10), in terms of both aberrant values frequency and amplitude (for example, computed β coefficient values range between 24 orders of magnitude in only 5 tests).

	-8 s
	Simulated time from initial state: 0.05 s
	Simulated duration between transitory states: 0.025 s
	Total elementary particles: 5/volume unit
	Maximum clusters size: N=150
	Initial state: [1]= 5 clusters/volume unit
	Noise: 2%

Table 1 .

 1 Deficiency data for each linkage class in N=4 reaction networks

	ABBREVIATIONS AND NOTATIONS
	(i)	Size i cluster (i elementary particles aggregated)
	(i+j)	Size (i+j) cluster
	(i)i≤N	Cluster sizes set
	[i]	Volumetric concentration for size i clusters
	[i]eq	Volumetric concentration for size i clusters at equilibrium
	[ i]	Evolution rate for volumetric concentration for size i clusters
	[ i]	Volumetric concentration for size i clusters measurement (random variable)
	([i])i≤N	Volumetric concentration distribution of clusters by sizes
	( [ i ] ) i ≤ N	Evolution rate for volumetric concentration distribution of clusters by sizes
	ki,j	Kinetic constant for i and j clusters aggregation into (i+j) clusters
	k * i,j	Kinetic constant for (i+j) clusters breakage into i and j clusters
	Ki,j	Equilibrium constant for i and j clusters aggregation into (i+j) clusters reaction
	Φi,j	Reaction rate for aggregation/breakage of i and j clusters
	N	Maximum clusters size (as a number of elementary particles)
	n	Total number of identical elementary (unbreakable) particles
	V	Total volume considered
	δi,j	Kronecker delta
	δ	Deficiency of population balance network
	l	Number of linkage classes of population balance network
	m	Number of complexes in population balance network
	s	Dimension of stoichiometric sub-space of population balance network
	K	Aggregation kernel

Table 2 .

 2 Equilibrium states conditions in N=4 reaction networks

	Clusters reactions reversibility
				Stable positive equilibrium states
	1+1→2	1+2→3	1+3→4	2+2→4

Table 3 .

 3 Example of potentially multi-stable aggregation/breakage kernel

	Maximum cluster size N = 5
	Kinetic constants
	k1,1=1	k * 1,1=592.48963
	k1,2=95.437698	k * 1,2=201.38279
	k1,3=419.0933	k * 1,3=2074.8396
	k1,4=1421.4501	k * 1,4=368.28533
	k2,2=34076.512	k * 2,2=753.42078
	k2,3=507535.96	k * 2,3=11190.471

Table 4 .

 4 Equilibrium states corresponding to the same population balance problem (kernels described in Table3, total elementary clusters concentration of 1 per volume unit)

	Clusters size	Equilibrium	state	1	Equilibrium	state	2
		(clusters/volume unit)		(clusters/volume unit)
	[1]	0.93783354			0.84858688
	[2]	9.8266x10 -3			1.6201x10 -2
	[3]	4.3674x10 -3			1.0742x10 -2
	[4]	4.3674x10 -3			1.0742x10 -2
	[5]	2.3881x10 -3			8.7629x10 -3

Table 9 .

 9 Examples of kernel parameters fitting based on disturbed equilibrium state data (1% error) -reproducibility tests using correct theoretical kernel models

		Best fitting aggregation kernel constants	Best fitting breakage kernel constants	Fitting criterion
		α	α	β	n	C
	Theoretical values	1	0.25	3	0.5	0
	N°1 test	1	0.1546	2400	3.23446	0.00051
	N°2 test	1	0.135	-1.6x10 -14	46.68	0.00299
	N°3 test	1	0.1278	1.34x10 89	114.103	0.0027
	N°4 test	1	0.134876	13.34	2.18957	0.00593
	N°5 test	1	0.134956	87733	5.16	0.0065858
	N°6 test	1	0.14008	-2.2369x10 15	50	0.014257
	N°7 test	1	0.166198	86.2013	2.06	0.021
	N°8 test	1	1.42x10 160	372	0.00169	0.016
	N°9 test	1	0.13748	7x10 235	211	0.0284
	N°10 test	1	0.1436	-286	6.63	0.0579
	N°11 test	1	0.132886	-100	5.49	0.03299
	N°12 test	1	0.113926	-5.80392	1.418	0.0565
	N°13 test	1	0.123744	-5.56x10 11	37	0.0893
	N°14 test	1	2.93	10.17	0.3929	0.0615
	N°15 test	1	0.1236	9.4x10 307	284	0.0357
	N°16 test	1	0.131805	2.4x10 149	156	0.147
	N°17 test	1	4.49x10 80	189	0.0028	0.181

Table 10 .

 10 Examples of kernel parameters fitting based on disturbed equilibrium state data (2% error)reproducibility tests using correct theoretical kernel models -T criterion

		Best fitting aggregation kernel constants	Best fitting breakage kernel constants
		α	α	β	n
	Theoretical values	1	1	3	1.5
	N°1 test	1099.2	6.47607x10 -27	28.8895	250.435
	N°2 test	1094.63	5.54897x10 -28	7.5477	470.113
	N°3 test	1098.46	2.52139x10 -27	28.8702	250.264
	N°4 test	49.9474	4.40021x10 -5	6.73752x10 -13	0.000330121
	N°5 test	30.8188	1.98245x10 -13	1.79745x10 -13	1.22582x10 -5

Table 11 .

 11 Example of kernel parameters fitting based on disturbed transitory states data (5% error)

	Kernel models	Best fitting aggregation kernel constants	Best fitting breakage kernel constants	Fitting criterion
	Aggregation	Breakage	α	α	β	n	T
	Turbulence 1	Exponential 1	4.05888	0.842453	1.45318	0.815055	3.91454
	Turbulence 1	Exponential 2	3.99292	7.77817	3.4312	0.214202	3.92658
	Turbulence 1	Power law	4.07621	0.68399	0.0399478	-	3.94365
	Turbulence 2	Exponential 1	1.32957	3.10842	17.3232	1.67806	4.48134
	Turbulence 2	Power law	1.572811	3.76433	1.318x10 -14	-	5.31992
	Turbulence 2	Exponential 2	4.14059	28.1234	12.1649	38.7169	6.7654
	Brownian	Exponential 2	27.2112	0.299807	3.9293x10 -12	1.37684	102.421
	Brownian	Exponential 1	25.7857	3.1725x10 -14	5.7908	0.185404	102.438
	Brownian	Power law	25.7857	8.27797x10 -15	0.00251735	-	102.438

Table 12 .

 12 Examples of kernel parameters fitting based on disturbed transitory states data (2% error-reproducibility tests using correct theoretical kernel models -TC criterion minimization) : effects of proximity to equilibrium state and time interval between measurements

		Best fitting aggregation kernel constants	Best fitting breakage kernel constants
		α	α	β	n
	Theoretical values	1	1	3	1.5
		Time from initial state 0.01s -Time interval 0.025s
	N°1 test	0.997673	1.18519	2.4185	3.19936
	N°2 test	0.961397	0.854704	1.86165	4.0272
	N°3 test	0.982993	1.24039	3.3828	2.21797
	N°4 test	0.94833	0.071937	5.11014	3.54243
	N°5 test	1.03542	1.15425	3.97135	2.38127
		Time from initial state 0.05s -Time interval 0.025s
	N°6 test	0.95	0.857545	6.44983	0.981753
	N°7 test	0.983844	0.934813	1.75802x10 12	26.2046
	N°8 test	1.01492	1.12928	25.8302	2.14533
	N°9 test	0.970671	0.914566	5.63656x10 -12	1.67278
	N°10 test	0.981762	0.922952	1.95157	3.49263
		Time from initial state 0.01s -Time interval 0.005s
	N°11 test	1.0203	3.54121	14.1176	0.703108
	N°12 test	0.978134	1.49193x10 -13	2.05187	4.75196
	N°13 test	1.06494	9.41688	477.607	3.29802
	N°14 test	1.00915	2.24547	12.0581	0.968625

ACKNOWLEDGMENTS

This research work was carried out in the context of French/German FLOTINC project (ref. ANR-15-CE08-0040), funded by ANR (Agence Nationale de la Recherche) and DFG (Deutsche Forschungsgemeinschaft).

PBS ( K,K*,Δt ) : ( [ i ] initial ) i≤ N → ( [i ] simul ) i≤ N [START_REF] Angeli | A Tutorial on Monotone Systems-With an Application to Chemical Reaction Networks[END_REF] TC ( K,K*,Δt,N,k, ( [ i ] initial ) i≤ N , ( [i]

Minimization of TC objective function with respect to K and K* is similarly carried out using Nelder-Mead heuristic approach.

In order to assess the effect of noise on the maximum likelihood corrected fitting procedure, a disturbed pseudoexperimental data set was generated from: Kernels:

-Turbulence type 1 aggregation, α=4

-Exponential law 1 breakage, n=0.5, α=1, β=3

Step size: 10

Turbulence type 1

Breakage

Kernel type Expression as functions of clusters sizes Constants Specific examples

Power law k i,j * =α⋅(i+ j ) ( n -1 ) n ≥ 0, α ≥ 0 Valentas [START_REF] Valentas | Analysis of Breakage in Dispersed Phase Systems[END_REF] Ramkrishna [START_REF] Ramkrishna | Drop-Breakage in Agitated Liquid-Liquid Dispersions[END_REF] Exponential law 1

Exponential law 2

Selomulya 53

Exponential law 1

)

Exponential law 2

)

Power law

Turbulence type 1 Exponential law 1

)

)

Power law

Turbulence type 2

Exponential law 1

Exponential law 2

⋅ ( j

⋅ ( (i + j )