Optimal operation guidelines for direct recovery of high-purity precursor from spent lithium-ion batteries: hybrid operation model of population balance equation and data-driven classifier

The direct resynthesis of precursor from spent lithium-ion batteries (LIBs) via co-precipitation is a crucial step in closed-loop cathode recycling systems. However, design and operation strategies for producing high-purity precursors have not been comprehensively explored or optimized. Herein, we propose the optimization of co-precipitation during the recovery of spent LIBs to achieve impurity-free precursor resynthesis. By incorporating the thermodynamic equilibrium model of the leaching solution of spent LIBs into a population balance equation (PBE) model, we identified the operating ranges that prevented the formation of impurities. Bayesian optimization was employed within the screened operating ranges to determine the optimal operating conditions for minimizing both operation time and maximum particle size. This optimization was performed for both unseeded batch and semi-batch systems. The results demonstrate that the selection of an optimal semi-batch operation can reduce the operation time by 23.33% and increase the particle size by 54.75%, owing to the high nucleation and particle growth rate during the initial time step. By employing an optimization approach based on the PBE model, this study provides detailed operational guidelines for batch and semi-batch co-precipitation, enabling the production of high-purity precursor materials from spent LIBs, while minimizing both operating time and maximum particle size.

1. Introduction

The ongoing energy transition towards renewable power is driving the expansion of the rechargeable battery market. Specifically, lithium-ion batteries (LIBs) have experienced remarkable expansion in the global battery market, surging from 12.9 billion USD in 2012 to 52 billion USD by 2021, with projections towards a further increase to 77.42 billion USD by 2024 (Asif & Singh, 2017). Among LIBs, LiNi_xMn_yCo_1−x−yO₂ (NCM) batteries, which currently occupy 69% of the global production market, are expected to surpass 80% (Curry, 2017). This dominance is attributed to the superior power rating and energy density compared with LiFePO₄ batteries. With the dramatic expansion of NCM production, the accumulation of spent LIB is estimated to range from 120000 to 170000 tons, solely within the electric vehicle sector (Nshizirungu et al., 2021). Given the limited reserves of NCM sources, recycling of the main metals, such as Li, Co, Mn and Ni, has emerged as the sole solution for achieving a closed-loop battery manufacturing cycle.

The spent NCM battery recycling industry currently relies on hydrometallurgy and pyrometallurgy processes with annual recycling capacity ranging from 1000 to 100000 tons of spent LIB (Baum et al., 2022). In pyrometallurgy, a battery undergoes thermal decomposition and reductive roasting to convert it into a metal alloy. Owing to the flexibility of feed and availability of existing facilities, pyrometallurgical processes are commonly employed in large-scale recycling of over 20000 tons of spent LIB per year (Baum et al., 2022). Owing to the formation of alloys, the recovery of valuable metals, including cobalt, nickel and lithium, of high purity is impossible (Makuza et al., 2021; Jeon, Kim, Eun et al., 2022). In contrast, valuable metals, including cobalt, nickel, manganese and lithium, can be selectively extracted with high purities via hydrometallurgy (Yao et al., 2018). After chemical leaching of cathodes, pH adjustment or liquid–liquid extraction are conducted to selectively extract metals as salts (Vasilyev et al., 2019; Fan et al., 2021; Kim, Kim et al., 2023; Liang et al., 2024; Li et al., 2021). Due to the mild operating conditions (80 °C, 1 bar) and high recovery efficiency, most research and industrial efforts have been directed towards advancing the hydrometallurgy processes, by optimizing leaching agents and chemical precipitants (Jung et al., 2021; Jeon, Kim, Kim et al., 2022). Despite the advantages of the selective extraction of valuable metals, secondary conversion processes are essential for regenerating cathodes for new batteries (Wan et al., 2022). This limitation hinders the direct application of the process in closed-loop battery manufacturing cycles.

To address this issue, recent experimental studies have proposed the direct regeneration of Ni_xCo_yMn_z(OH)₂ precursors from leaching solutions (Fang et al., 2022; Lee et al., 2023; Yang et al., 2018; Yu et al., 2023). Typically, ammonia is injected into leaching solutions as a complexing agent to obtain metal–ammonia complexes. Under basic conditions (pH > 10), these complexes prevent undesired precipitation reactions (e.g. metal hydroxides) (van Bommel & Dahn, 2009) and release metal cations to form spherical hydroxide particles:

Fang et al. (2022) added NH₃·H₂O (0.4 mol L⁻¹) to the reactor after adjusting the pH to 8.0, to trigger co-precipitation. Following a 12 h batch operation, a spherical Ni_0.5Co_0.2Mn_0.3(OH)₂ precursor with a diameter of 10–15 µm was successfully regenerated from spent LIBs. Park et al. (2019) investigated the impact of Fe as an impurity in regenerated precursors to test the feasibility of the direct resynthesis of precursors from spent LIBs. Co-precipitation was performed in a simulated leaching liquor with varying Fe compositions ranging from 0 to 1.0%. The gradual increase in the Fe concentration of the leaching solution resulted in the deterioration of structural perfection, ultimately leading to a degradation in the capacity and Coulombic efficiency of regenerated cathodes. Consequently, for cathodes derived from a 1% Fe leaching solution, the discharge capacity and Coulombic efficiency were 3.2 and 2.5% lower than those of cathodes obtained from pure leaching solutions, respectively.

Previous experimental studies have successfully proved the feasibility of the direct resynthesis of precursors from spent LIBs. However, industrial application of the resynthesis concept still faces challenges in process operation. First, impurities, including Al, Fe and Cu ions, should be filtered to ensure the yield and quality of precursors. Impurity ions are commonly precipitated by increasing the pH, but this also results in the precipitation of the target metals, such as Ni, Co and Mn, owing to the overlapping pH range of metal hydroxide precipitation (Song & Zhao, 2018). The loss of the target metals in this step increases the operating cost because addition of pure metal sulfates (NiSO₄, CoSO₄ and MnSO₄) would be required in the co-precipitation process for the crystallization of the target precursors, including Ni_0.5Co_0.2Mn_0.3(OH)₂ and Ni_0.3Co_0.3Mn_0.3(OH)₂ (Feng et al., 2018). Additionally, there is still the chance of undesired precipitation from the other ions in the leaching solutions, including SO₄²⁻, Na⁺ and NH₄⁺, even after impurity removal. Second, the dynamics of crystallization (e.g. nuclei formation and secondary particle growth) in the co-precipitation steps are strongly affected by the operating conditions and complex equilibrium of the other ions in the leaching solution. The ammonia injection, temperature and NaOH addition (for pH control) should be carefully controlled to minimize the particle size under minimum operation time, while constraining undesired side reactions.

To address the aforementioned challenges, reliable modeling is essential to capture the dynamics of particle size evolution and reduce the need for extensive experiment. Once the model is established, feasible operating conditions can be identified, allowing for operation optimization. However, unlike conventional precursor synthesis, resynthesis from the spent LIB introduces modeling complexities due to the presence of impurity ions. During crystallization, complex equilibria between metal ions, the chelating agent and impurities affect the kinetics of nucleation and particle growth. The intricate interaction between kinetics and thermodynamic equilibria makes it difficult to track particle size evolution and complicates operation optimization while preventing impurity formation.

To overcome these issues, we developed a population balance equation (PBE) model of crystallization, integrated with the thermodynamic equilibrium model of complex equilibria, salts and ion dissociation within the leaching solutions. Using the developed model, we first focused on determining the feasible operating conditions that constrain the precipitation of salts other than the precursors. Under these feasible conditions, the optimization aimed to minimize both operation time and particle size distribution. Smaller and narrower particle size distributions are typically more desirable, as they reduce Li-ion transport distances and increase the cathode's specific surface area, leading to improved charge/discharge rates and overall battery performance. Thus, the operating trajectories were optimized to achieve shorter operation times and optimal particle size distribution.

The high nonlinearity of the PBE model makes it impossible to represent the feasible domain of operating conditions using linear hyperplanes. This complexity significantly increases the search time for feasible operating conditions and leads to higher computational costs during the optimization procedure. To address this limitation, this study employed surrogate modeling using a deep neural network (DNN). After the DNN has been trained on a dataset of operating conditions and impurity formation, the pre-trained model predicts the probability of a given operating condition belonging to the feasible domain. This allows the optimizer to pre-screen operating conditions without requiring full PBE simulations, thus reducing both computation time and the search space for feasible solutions.

Additionally, in terms of operation strategies, both batch and semi-batch modes offer distinct advantages and limitations. Batch operation is easy to control but can lead to longer operational times due to decreasing supersaturation, especially with fixed feed concentrations. Semi-batch operation allows high supersaturation to be sustained by adjusting chelating agents or temperature but increases the risk of impurity formation and unfavorable increase of particle size, which can affect NCM cathode performance. As this is the first study to model and optimize precursor resynthesis, both unseeded batch and semi-batch modes were explored, with a performance comparison conducted to highlight the strengths of each mode. The novelty and major contributions of this study are as follows:

(i) A thermodynamic equilibrium model of the leaching solution is used to screen the feasible operation range of the resynthesis process. The formation of numerous impurities helps determine the reliable operating conditions for obtaining high-purity precursors from a leaching solution. By integrating the equilibrium model with the PBE model, simulations can be conducted to track the evolution of the particle size distribution (PSD) while considering the intricate equilibria of ions during leaching.

(ii) Optimization can provide the operating trajectories of chelating agent concentration, flow rate and temperature, for minimizing the operating time and particle size. The optimal trajectory driven by this model-based approach provides detailed guidelines regarding the operation of the resynthesized precursor. Additionally, by pre-screening the infeasible operational domains, the optimal trajectory can ensure the production of a theoretically pure precursor in the presence of impurity ions.

(iii) Two different operation options (batch and semi-batch) were optimized within the feasible operating domains. A comparison of the optimal PBE dynamics and performances of different options can help in selecting the appropriate operation to achieve the target performance (e.g. reducing the operation time and minimizing the particle size).

2. Mathematical modeling

Fig. 1 illustrates the process of resynthesizing the NCM precursor from the cathode of a spent LIB. Initially, the spent NCM cathode undergoes sulfuric leaching for conversion into a metal slurry. Prior to precursor resynthesis, metal impurities in the slurry, such as Al³⁺ and Cu²⁺, are filtered out as metal hydroxides at pH values above 9. Subsequently, the purified slurry is fed into a reactor with aqueous NH₄OH (aq) for co-precipitation. Controlling the NH₄OH (aq) feed facilitates the formation of metal–ammonium complexes, while preventing the formation of metal hydroxides. These metal complexes co-precipitate above the equilibrium concentration, forming NCM precursor particles. Unlike conventional precursor synthesis using pure metal sulfates, co-precipitation from leaching solutions requires consideration of the equilibrium shifts among impurity ions during PSD evolution.

Figure 1 Schematic diagram of the resynthesis of NCM precursors from spent NCM cathodes.

Thus, in this section, a PBE model of co-precipitation, incorporating the thermodynamic equilibria of metal ions and impurities (e.g. Na⁺ and SO₄²⁻) in leaching solutions, is developed. As detailed in Section 2.1, the PBE model governs the PSD evolution through the nucleation and 1D particle growth rates. By integrating the PBE model with the thermodynamic equilibria of the leaching solutions, the mass and energy balances in the reactor were iteratively solved in discretized time steps, as detailed in Section 2.2. The kinetic parameters of the PBE model were estimated using experimental data, as explained in Section 2.3. Additionally, the operating conditions affect the formation of impurities during co-precipitation. Hence, a feasible operating domain to ensure impurity-free precursor production is identified in Section 2.4. Furthermore, an additional data-driven classifier model was trained to selectively screen feasible operating conditions. To optimize the operating trajectory, the trained model was utilized to filter out feasible operating conditions. To simplify the problem, the molar composition of the leaching solution per kilogram of the NCM cathode material is utilized in Section 2.4, and further optimization is presented in Section 3. Details regarding the operating conditions of the leaching process and composition of the leaching solution are provided in Section S1 of the supporting information.

2.1. Population balance model

Under the assumption of semi-batch and unseeded operations, the PBE describing the evolution of the population density under nucleation and particle growth is given by equation (3) (Rawlings et al., 1993):

where V is the slurry volume, f is the population density of the particles, L is the particle size, G is the particle growth rate, B is the nucleation rate and δ(L) is the Dirac delta function. The Dirac delta function is used for establishing spatial boundary conditions, signifying the occurrence of nucleation at the minimum nucleus size, L₀.

The particle growth rate depends on the temperature, particle size and supersaturation ratio, and it can be expressed by equation (4):

where k₀, E_g, k₁, k₂ and α are parameters estimated from experimental data. S is the supersaturation ratio, which is the driving force of co-precipitation given by (Mersmann, 2001)

where [M(NH₃)_n]²⁺ is the total concentration of the metal–ammonia complex (M is Ni, Mn or Co) and K_sp is the solubility product of the metal hydroxides.

The nucleation rate can be expressed as follows:

where k_b and β are parameters estimated from experimental data.

Under the batch and semi-batch operations, the mass balance of the solute is defined as follows (Hu et al., 2004):

where c_i is the solute concentration in the slurry, Q_in and c_in are the flow rate and concentration, respectively, of the input stream, and ρ_c and k_v are the precursor density and volumetric shape factor, respectively.

The energy balance of the slurry is defined as follows:

where ρ and c_p are the density and specific heat of the slurry, respectively; Q_k and H_k are the flow rate and enthalpy of the input stream, respectively; is the heat of co-precipitation; and H_ext is the heat removed by the cooling system.

2.2. Model solution with thermodynamic equilibrium

The evolution of the PSD should be trackable for model validation and operation optimization. Thus, instead of using moment methods (Luo et al., 2017; Hulburt & Katz, 1964; Brock & Oates, 1987), the solution technique reported by Hu et al. (2004) was employed to simulate the PSD evolution in each time step. In the discretized PBE [equation (3)], the population balance and particle growth can be converted into a set of algebraic equations, simplifying the process compared with solving complex partial differential equations. Formulation was performed only in batch operations with a constant slurry volume. However, the experimental data were collected during a semi-batch operation (Feng et al., 2018). Thus, in this study, a set of algebraic solutions was reformulated as a semi-batch operation for parameter fitting. Details of the reformulation can be found in Section S2 of the supporting information.

Consequently, the population density of the precursor is expressed as follows:

where the first index j denotes the time interval, the second index i indicates the series number of particle size points, and L_j,i and V_j are the particle size points and slurry volume at time interval j, respectively.

L_j+1,i is given by the following equation:

To update the PSD according to equations (9) and (10), the supersaturation ratio [equation (5)] should first be estimated in the equilibrium state of the leaching solution system. In contrast to conventional precursor synthesis systems based on metal sulfates, complex equilibria involving impurities present challenges for this process. The concentration of the metal–ammonia complex [M(NH₃)_n]²⁺ [equation (5)] can be affected by the side reactions of metal ions with impurities, which ultimately leads to the depletion of metal resources for precursor production. Additionally, under various operating conditions, there is a possibility of undesired generation of impurities. To address these complexities in the thermodynamic equilibrium modeling of the system, a Pitzer–Debye–Huckel model (Chang & Lin, 2019) based equilibrium solver within the commercial process simulator Aspen Plus V12 (Al-Malah, 2022) was used. This solver considers the metal–ammonia complex equilibria and all potential side reactions involving ionic equilibria and salt formation. Further details on the reactions are provided in the supporting information (Section S3).

The overall simulation algorithm for the PBE model is shown in Fig. 2. The operating conditions including the temperature, concentration and flow rate of NH₄OH (aq) in the thermodynamic model were updated according to the molar composition of the leaching solution. A process simulator computed the mass and energy balance of the thermodynamic equilibrium in the leaching solution. To ensure convergence, a maximum of 500 iterations were considered with a relative error tolerance of 0.01 for computing the mass and energy balance in the simulator. If the equilibrium simulation failed to converge, the PBE simulation was terminated. Otherwise, the supersaturation ratio (S) and pH at equilibrium were calculated. To maintain the pH at the target value (pH_target), the input volume of NaOH (aq) is iteratively updated using parameter α until the current pH meets the error criteria of pH_target. Here, for faster convergence, an α value of 5 is used, while allowing for an error_pH of 0.05. Once the pH condition is satisfied, the PSD is updated using a discretized PBE model [equations (9) and (10)] at the jth time step with an interval of Δt. In this study, a time step of 1 min is chosen. On the basis of the updated PSD, the process simulator solved the mass and energy balances of the system [equations (7) and (8)] at thermodynamic equilibrium. Finally, depending on the updated molar concentration of the system, the preceding calculation steps were iteratively performed either until the end of the simulation time or until the S value was less than 1, while ensuring the mass and energy balances of the process simulator.

Figure 2 Simulation framework of the PBE model integrated with a thermodynamic model.

2.3. Parameter estimation

To estimate the kinetic parameters in the PBE model, in situ data of the median diameter of the secondary particles of the NCM 111 precursor [Ni_1/3Mn_1/3Co_1/3(OH)₂] at a pH of 10.6 were used (Feng et al., 2018). Additionally, the maximum particle size of the secondary particles after co-precipitation and concentration of the precursor was used for parameter estimation. The following experimental setups were considered as the operating conditions of the simulation: 1 bar and 60 °C, with flow rates of 20 mL h⁻¹ of 2.0 M MSO₄ (aq) and 15 mL h⁻¹ of 5.0 M NH₃ (aq); additionally, a 4.0 M NaOH solution was added to maintain a pH of 10.6, and a semi-batch operation was conducted over a duration of 180 min. Here, seven unknown parameters, namely k₀, E_g, k₁, k₂, α, k_b and β, are estimated using three different types of experimental data (median diameter, precursor concentration and maximum particle size). The PBE model is highly nonlinear and has a high computational cost owing to its integration with a thermodynamic solver for the calculation of thermodynamic equilibrium. Thus, Bayesian optimization (BO) (Martinez-Cantin, 2014), which demonstrates superior performance within resource limitations, was chosen over an evolutionary algorithm. BO operates in a feedback loop of two main steps: surrogate modeling of targeted objective functions, and selecting the next searching point via an acquisition function. The surrogate modeling is commonly built upon a Gaussian process (GP) with posterior mean and variance based on previously searched points. The acquisition function then uses the GP to identify the next optimal point for exploration, balancing exploration of uncertainty and exploitation of searched regions. This reduces the number of function evaluations required to find an optimal solution. In contrast, evolutionary algorithms typically require a large number of evaluations to search the entire space, as they rely on population-based exploration methods without incorporating probabilistic knowledge of the search space. As a result, evolutionary algorithms tend to be less efficient, especially in high-dimensional or computationally expensive problems. BO's ability to quantify uncertainty makes it especially well suited for cases like this, where each function evaluation is computationally costly and the search space is highly nonlinear (Kim, Han et al., 2023).

Given the parameter ranges listed in Table 1, BO was implemented to minimize the error in the simulation results, and it is formulated as follows:

where θ and D are the parameter and its domain, respectively; the first term is the mean squared error (MSE) of the median particle size between experiment () and simulation (), the second term is the error of the precursor concentration after co-precipitation, and the third term is the maximum particle size error after co-precipitation. All these terms are combined into a single objective function using w₁, w₂ and w₃.

Table 1
Kinetics parameters and ranges for parameter estimation

Parameter	Range	Units
Eg	[10000–70000]	J mol^−1
k0	[0–2000]	m min^−1
k1	[0–1000]	m^−1
k2	[0–10]	Dimensionless
kb	[10000–200000]	Particle number (kg^−1 min^−1)
α	[0–2]	Dimensionless
β	[0–2]	Dimensionless

Owing to the different scales of the experimental data, w₁, w₂ and w₃ in equation (11) were set to 1, 1000 and 1, respectively. In this study, BO employed the upper confidence bound (UCB) acquisition function (Kaufmann et al., 2012), denoted as a(θ; λ) = μ(θ) + λσ(θ), where λ is set to 2.5. The choice of λ emphasizes the exploration of parameter regions in light of the highly nonlinear nature of the PBE model.

2.4. Feasible operating domain

In the presence of various ions, such as Cu²⁺, Al²⁺, Fe²⁺, Na⁺ and SO₄²⁻, along with the target metal ions, there are numerous chances for the formation of impurities, including metal hydroxides [e.g. Al(OH)₃, Cu(OH)₂ and Fe(OH)₂] and metal sulfates (e.g. CoSO₄, MnSO₄, NiSO₄, Al₂SO₄, Na₃SO₄OH and Na₂SO₄) in the leaching solution under basic conditions. By targeting the pH required for the precipitation of impurities into metal hydroxides, impurities such as Cu²⁺ and Al²⁺ can be screened out without the loss of target metal ions, particularly at a pH of approximately 10. However, during the co-precipitation of the metal slurry, impurity salts can persist owing to the presence of Na⁺ and SO₄²⁻ in conjunction with the chelating agent (NH₄⁺), as governed by the complex equilibria listed in the supporting information (Table S2). Thus, the selection of operating conditions, including the concentration and flow rate of the chelating agent and operating temperature of the NCM precursor resynthesis process, not only affects the dynamics of the PSD evolution but also influences the formation of impurity salts. Thus, in this study, prior to optimizing the operating trajectories, we identified the feasible domain of the operating conditions using the thermodynamic equilibrium model described in Section 2.2, to ensure the absence of impurities. To simplify this problem, the pH of the leaching solution was maintained at 10.6. The molar composition of the cathode material and operating conditions of the leaching process can be found in Section S1 of the supporting information.

To identify the feasible domains, the operating conditions are firstly sampled from the following ranges with 8000 samples. NH₄(OH) (aq) concentration: 0–10 M; NH₄OH (aq) flow rate ratio multiplied by leaching volume: 0–10; and co-precipitation temperature: 25–90 °C (Bommel & Andrew, 2009; Feng et al., 2018; Barai et al., 2019). Latin hypercube sampling (LHS) [equation (12)] (Loh, 1996), which divides each range into equivalent probabilities, was chosen for sampling to ensure a broad investigation of the equivalent operating domain. For a given leaching solution under the sampled conditions, the required volume of 4 M NaOH (aq) for a pH of 10.6 was calculated using the thermodynamic equilibrium models, and the model finally calculated the moles of solid impurities:

where represents the samples for dimension j, π_j is a random permutation of indices, u_j is an independent and identically distributed sample of the uniform distribution, and n is the number of samples.

The nonlinearity of thermodynamic equilibrium makes representing a feasible domain using linear hyperplanes impossible. This complexity poses a challenge in screening out the infeasible operating domain during the optimization of the operating trajectories of crystallization. To address this limitation, a DNN was employed to classify whether a given operating condition fell within the feasible domain. In the classifier, the inputs comprised the operating conditions, including the concentration and flow rate of NH₄(OH) (aq) and temperature. The output was the probability of belonging to a feasible domain. For training and validation, the 8000 samples were randomly divided into training (80%) and validation (20%) sets. The network architecture was constructed by stacking five fully connected layers with 128, 64, 32, 16 and four hidden nodes, employing leaky-ReLU activation. In the output layer, a sigmoid function was used as the activation function. Training was conducted to minimize the binary cross-entropy (BCE) loss []. For training, we used an adaptive moment-estimation optimizer with a learning rate of 0.0001. Training and validation of the DNN were performed using PyTorch version 2.1.1 (Paszke et al., 2019).

where y_i represents the label of the ith sampled operating condition (zero for infeasible and one for feasible), represents the predicted probability of belonging to a feasible domain and N represents the batch size of the training dataset.

3. Problem formulation

The co-precipitation of the leaching solution involves several factors, including PSD, metal recovery from the leaching solution and operation time. To minimize the operation time with a certain level of metal recovery, maintaining a high supersaturation ratio by operating at the upper limits of the feasible domain may seem reasonable owing to the high rates of nuclei formation and particle growth rate, represented in equations (4)–(6). However, owing to the rapid growth of particles and consistent formation of large nuclei, a narrow PSD is unachievable. A large particle size increases the operating time of the li­thia­tion process and degrades the electrode performance (Zhang et al., 2022). Consequently, relying solely on a high supersaturation ratio is impractical for controlling the properties of the PSD.

To effectively optimize both the operational time and PSD, the operating conditions should be adjusted in accordance with the evolution of the PSD. Optimal trajectories of operating conditions exist to minimize the operation time while satisfying the specific criteria for the required PSD. In this study, the simulation framework outlined in Section 2 was employed to optimize the operating conditions. Additionally, the selection of the operation mode can affect the operational time, metal recovery and PSD, whether using batch or semi-batch operation. Consequently, optimization was carried out separately for both batch and semi-batch operations.

3.1. Batch operation optimization

For optimization, three operating conditions were considered: temperature of co-precipitation, injection concentration and volume of the chelating agent, NH₄OH (aq). In batch operations, the chelating agent is initially fed into the leaching solution at a given concentration and volume, and co-precipitation proceeds while the given temperature is maintained. The optimization aims to minimize both the operation time required for the complete recovery of metal ions from the NCM precursor and the PSD of the final product. In the dynamics of the NCM precursor, it is noted that the maximum particle size is positively associated with the coefficient of variation in the PSD, as referenced in Section S4 of the supporting information. Thus, rather than adopting a multi-objective approach, a single objective function was formulated through the summation of operational time (min) and maximum particle size (µm). This ensured the convergence of the optimization in the presence of high-dimensional decision variables and the high nonlinearity of the PBE model. The operating conditions were constrained to the feasible domain, to facilitate the convergence of the PBE simulation. In this study, a penalty method was employed to address constraint violations (Freund, 2004).

Therefore, the operating condition optimization problem is formulated as follows:

where t is the co-precipitation time; d_max is the maximum particle size after co-precipitation; T, c and V are the temperature, concentration and volume flow rate of NH₄OH (aq), respectively; D_feasible is the feasible domain; error_PBE is the mass and energy balance error of the PBE model; and p_infeas and p_error are the penalty values.

BO was employed with the UCB acquisition function, as described for parameter estimation in Section 2.3. To ensure that simulations are conducted only within the feasible domain (D_feasible), the data-driven classifier model trained in Section 2.4 was utilized to filter out the operating conditions during optimization. Thus, the new operating conditions obtained via the optimizer were fed into the classifier. If the operating condition is determined to fall within the infeasible domain, then the penalty term (p_infeas) is returned to the optimizer without conducting a PBE simulation, prompting the optimizer to search for the next operating condition. This approach guarantees that the BO dynamically explores only the feasible operating conditions throughout the optimization loop. Additionally, to address the PBE simulation convergence issues, an additional penalty term was included in the objective function. If the PBE simulation terminates because of mass or energy balance errors in the process simulator (Fig. 2), then the penalty term p_error is added to the objective function. By imposing high positive values of p_infeas (1000) and p_error (100), the optimizer searches for the optimal operating conditions that ensure impurity-free production of the NCM precursor, while achieving PBE simulation convergence.

3.2. Semi-batch operation optimization

In a semi-batch operation, additional chelating agent injection during co-precipitation is permitted with temperature control. Beyond the selection of operating conditions, optimization of the operating trajectory can be effective in controlling the particle growth and nucleation rate by regulating the supersaturation ratio and temperature. Herein, the complexity of optimization was reduced by discretizing the operating trajectory into five points corresponding to metal recovery percentages ranging from 0% to 100% with intervals of 20%. The operating conditions considered were the temperature, concentration and volume flow rate of NH₄OH (aq). Consequently, with the three types of operating conditions, each comprising five points, 15 decision variables were considered for optimization. The objectives and constraints were the same as those of the batch operation optimization described in Section 3.1.

The operating trajectory optimization problem is formulated as follows:

where i denotes the discretized interval point of metal recovery; t_i is the co-precipitation time at the ith metal recovery point; d_max is the maximum particle size after co-precipitation; T_i, c_i and V_i are the temperature, concentration and volume flow rate, respectively, of NH₄OH (aq) at the ith metal recovery point; and p₁ and p₂ are the penalty functions, as described in equations (15) and (16), respectively.

For equivalent comparison with the batch operation optimization, the semi-batch operation optimization follows the same optimization procedure with identical penalty values for constraint handling, as described in Section 3.1.

4. Results and discussion

In Section 4.1, the results of the parameter estimation are presented through a comparison between the PBE simulation and experimental results, focusing on the in situ measurement data of the median diameter, maximum particle size and concentration of the particles after co-precipitation. Section 4.2 details the identification of the feasible domain and training results of the operating condition classifier. Finally, in Section 4.3, operating optimization is conducted using the PBE model integrated with the trained classifier. This section compares the results obtained from two different operating modes: batch and semi-batch.

4.1. Parameter estimation result

For parameter estimation, BO was implemented to minimize the summation of three different errors [equation (11)]: median particle size, maximum particle size and concentration of the precursor after co-precipitation. The parameters obtained after 1000 iterations of the BO step are listed in Table 2. With an MSE of 9.85, the median particle size evolution was in reasonable agreement with the experimental data, as depicted in Fig. 3(a). Furthermore, the precursor concentration and maximum particle size at the culmination of co-precipitation were calculated to be 91.8 mM and 63 µm, respectively, while the corresponding experimental data were measured to be 85 mM and 50 µm, respectively [Figs. 3(b), 3(c)]. The PSD after co-precipitation did not perfectly align with the experimental data, owing to the limitations of the 1D PBE model. To avoid the precise kinetics of the PSD evolution, high-fidelity model development supported by in situ experiments is necessary. However, this entails a significant increase in computational costs, rendering it impractical for subsequent operational optimization. Hence, the current results regarding the precursor concentration indicate that the fitted parameters offer a reasonable mass balance for co-precipitation, accompanied by a low MSE error during the evolution of the PSD. Operating optimization was conducted using the PBE model with the fitted parameters, and the corresponding results are described in Section 4.3.

Table 2
Parameter fitting results

Parameter	Value
Eg	61041.97
k0	588.51
k1	161.65
k2	1.26
kb	176404.20
α	0.38
β	1.36

Figure 3 PBE model simulation under optimal parameters: (a) median particle size evolution, (b) concentration of precursor and (c) PSD evolution.

4.2. Identification of feasible operating domain

Fig. 4 presents the calculated results of the thermodynamic equilibrium for the sampled operating conditions described in Section 2.4. The scattered points represent the operating conditions associated with the formation of the solid impurities. Within the infeasible domain, distinct types of impurity formation become apparent with the formation of the following solid impurities: CoSO₄, CoSO₄·6H₂O, Na₂SO₄ and (NH₄)₂SO₄. Due to the predominant concentration of Na⁺ and SO₄²⁻, the formation of Na₂SO₄ was widely distributed in the infeasible domain. Furthermore, the high equilibrium con­stants of CoSO₄ and CoSO₄·6H₂O contributed to unavoidable Co²⁺ losses throughout the entire infeasible domain. Increasing both the concentration and volume of aqueous NH₄(OH) led to the formation of (NH₄)₂SO₄. During co-precipitation, this solid formation not only degrades the purity of the NCM precursor but also hinders the formation of metal–ammonia complexes, limiting rapid nucleation and particle growth owing to the low supersaturation ratio.

Figure 4 Infeasible operating domain for resynthesis from the leaching solution with different types of impurity formations. Type 1: CoSO4·6H2O and Na2SO4; type 2: Na2SO4; type 3: (NH4)2SO4 and CoSO­4·6H2O; type 4: (NH4)2SO4; type 5: CoSO4 and Na2SO4; type 6: (NH4)2SO4 and CoSO4.

The DNN model was trained using 8000 samples to represent the nonlinearity of the identified infeasible domains. Instead of predicting the types and quantities of impurities produced, this study employed a model to classify whether the given operating conditions fall within an infeasible domain. Thus, the BCE loss for both training and validation sets, plotted in Fig. 5, shows sharp decreases from 0.73 and 0.75, respectively, to 0.01. The low dimensionality of the input and simplicity of the model contributed to this rapid reduction in loss. Fig. 6 shows the confusion matrix for the trained classifier, indicating a high accuracy of 99.75%. This level of accuracy underscores the reliability of the classifier in filtering the operating conditions for the co-processes, ensuring that only the NCM precursor is produced, while excluding the formation of solid impurities. The trained model was employed to filter out the feasible operating conditions during operating optimization, as described in Section 3.

Figure 5 BCE loss evolution of the solid impurity classifier.

Figure 6 Confusion matrix of the trained operating condition classifier.

4.3. Operating optimization

4.3.1. Batch operation

Fig. 7 illustrates the search trajectory of the operating conditions during BO. In the pre-trained classifier, penalization for selecting an infeasible region aids in reducing the number of searches in that region during optimization. Consequently, the search points converged to the optimal points, as listed in Table 3, after 196 functional evaluations. In the infeasible domain, most of the searched points were concentrated at NH₄OH (aq) concentrations exceeding 6 mol L⁻¹. Although a high concentration can boost metal conversion into the precursor, it is unsuitable for leaching solutions owing to the formation of impurities. Hence, the optimizer selected a concentration at least 51.4% lower than the infeasible point to ensure impurity-free production. To validate the optimal case, we selected 19 base cases with objective functions at most 20% higher than the optimal case. The corresponding operating conditions were simulated and compared with the optimal results in terms of PSD, recovery and purity.

Table 3
Optimization results of batch operation

Optimal operating condition	Value
c (mol L^−1)	2.92
V (L)	540.16
T (°C)	64.36

Optimization results	Value
t (min)	60
dmax (µm)	24.53
Metal recovery	0.78

Figure 7 Search space of operating conditions of batch mode via Bayesian optimization.

In the base cases, the operating conditions spanned the following range: c from 1.53 to 2.97 mol L⁻¹, V from 214.83 to 637.77 L and T from 56.36 to 75.73 °C. Fig. 8 illustrates the co-precipitation dynamics for both optimal and base cases. Regarding c, the optimal case deviated by only 1.68% from the upper limit of the base cases. A high c leads to an increased concentration of metal–ammonium complexes, which in turn results in a higher supersaturation ratio. However, the selection of a large V in the optimal case reduced the overall concentration of metal complexes in the leaching solution, resulting in an intermediate supersaturation ratio, as shown in Fig. 8(a). Particle growth and nucleation rates were determined using the supersaturation ratio and temperature, as described in equations (4) and (6), respectively. In the optimal case, the T value fell within the range of the base case. Consequently, the optimal operating conditions yielded intermediate nucleation and particle growth rates compared with those of the base cases, as shown in Figs. 8(b) and 8(c).

Figure 8 Comparison of the PBE dynamics of the optimal and base cases of batch operation: (a) supersaturation ratio, (b) nucleation rate and (c) maximum growth rate.

Fig. 9 illustrates the evolution of the particle size and recovery of the precursor from the leaching solution in both base and optimal cases. The base cases ranged from a minimum operation time of 58 min to a maximum of 81 min, with corresponding maximum particle size in the PSD of 35.57 and 22.26 µm, respectively, as depicted in Fig. 9(a). Conversely, in the optimal case, employing intermediate nucleation and growth rates effectively minimized both the maximum particle size and operation time, as detailed in Table 3. Moreover, compared with the base cases, the optimal case yielded competitive metal recovery from the leaching solution, achieving 78% recovery of precursors with 100% purity.

Figure 9 Comparison of batch performance of the optimal and base cases: (a) maximum particle size and (b) recovery of the NCM precursor.

4.3.2. Semi-batch operation

In contrast to the batch operation mode, semi-batch operation permits the additional injection of NH₄OH (aq) feed with temperature control during operation, resulting in the optimization of 15 operating conditions. The search trajectory for the average operating conditions is shown in Fig. 10. Compared with the batch mode, the semi-batch mode involves a higher dimension of decision variables, leading to convergence after 528 function evaluations despite the support of a pre-trained classifier. Notably, the infeasible points were mostly concentrated at NH₄OH (aq) concentrations higher than those of both optimal and base cases. The optimization results for the semi-batch operation are summarized in Table 4. The optimal case was validated by comparing its simulation results with those of 18 base cases with objective functions at most 50% higher than those of the optimal case.

Table 4
Optimization results of semi-batch operation

Optimal operating trajectory	Value
ci (mol L^−1)	[8.99, 0.34, 1.60, 0.52, 5.64]
Vi (L)	[49.95, 237.95, 217.81, 166.45, 247.82]
Ti (°C)	[75.70, 76.58, 57.73, 71.52, 51.08]

Optimization results	Value
t (min)	46
dmax (µm)	37.96
Metal recovery	0.80

Figure 10 Search space of the operating conditions of the semi-batch mode via Bayesian optimization.

As illustrated in Fig. 11, the operating trajectory of the 18 base cases is depicted as a bounded area, whereas the trajectory of the optimal case is enclosed within this area. This indicates that the optimization process effectively balances the operating conditions to minimize both the operation time and particle size within the boundaries of the base cases. For instance, during the initial stages of co-precipitation (metal recovery < 20%), NH₄OH (aq) was injected at high concentration and low volume, nearing the upper and lower bounds, respectively, as observed in Figs. 11(a) and 11(b). However, for metal recovery of 20–80%, the optimal case selects a high volume with concentrations near the lower bounds. Additionally, in terms of temperature, the optimal case maintained a high operating temperature until a metal recovery of 40%, then shifted to a lower value during the intermediate stage (metal recovery from 40% to 60%), as shown in Fig. 11(c).

Figure 11 Optimal operating trajectory of semi-batch mode: (a) NH4OH injection concentration, (b) NH4OH injection volume and (c) operating temperature.

Given the trajectories of the base and optimal cases, Fig. 12 shows the dynamics of the PBE in semi-batch operation mode. The supersaturation ratio in the optimal case followed the concentration trajectory, initially peaking at 27.85, before dramatically decreasing to 3.99 owing to the injection of a low-concentration NH₄OH (aq) solution at a high volume, as depicted in Fig. 12(a). Because the nucleation and growth rates are determined by the temperature and supersaturation ratio, these rates follow the trajectory of the supersaturation ratio, as shown in Figs. 12(b) and 12(c). Sustaining high nucleation and growth rates in the early stages facilitates the efficient recovery of the metal ions during leaching into the precursor, thereby reducing the operation time. However, maintaining a high growth rate throughout the operation is inefficient for particle size minimization. Thus, a significant shift to lower rates helps restrict particle growth during the remaining operation time.

Figure 12 Comparison of the PBE dynamics of the optimal and base cases of semi-batch operation: (a) supersaturation ratio, (b) nucleation rate and (c) maximum growth rate.

To support the effectiveness of the optimal trajectory, Fig. 13 presents the progress of the maximum particle size and metal recovery for both the optimal and base case trajectories. The base cases span from a minimum operation time of 54 min to a maximum of 88 min, with corresponding maximum particle sizes in the PSD of 43.61 and 28.23 µm, respectively, as depicted in Fig. 13(a). In contrast, the optimal case effectively balanced both operation time and particle size, resulting in a 14% reduction in operation time compared with the base cases. Although the maximum particle size was 34.47% higher than the minimum particle size among the base cases, the optimal case achieved a 47.78% reduction in operation time compared with the corresponding base case. Additionally, in comparison with the base cases, the optimal trajectory demonstrated competitive metal recovery from the leaching solution, achieving 80% precursor recovery with 100% purity.

Figure 13 Comparison of the semi-batch performance of the optimal and base cases: (a) maximum particle size and (b) recovery of the NCM precursor metal.

4.3.3. Comparison of batch and semi-batch operation

In the preceding section, the effectiveness of the optimal operating conditions was demonstrated through comparisons with base cases across batch and semi-batch operation modes. To find out which operation mode is the most feasible for the co-precipitation of the leaching solution in terms of minimizing both the particle size and operation time, this section compares the dynamics and performance of different operation modes under optimal operating conditions.

In the batch operations, the initial supersaturation and temperature values determine the dynamics of particle nucleation and growth, as shown in Fig. 14(a). For instance, selecting high initial supersaturation and temperature values can effectively decrease the operational time but increase the maximum particle size. Specifically, the growth rate [equation (4)] is proportional to k₂ (1.26) of the particle size term. This implies that a large particle size at the initial time leads to fast growth of particles owing to operating conditions, but it cannot guarantee minimum particle size after crystallization. Consequently, the optimizer selects intermediate points as the optimal conditions, as listed in Table 3, resulting in lower nucleation and growth rates than those in the semi-batch mode, as shown in Figs. 14(b) and 14(c). In contrast, the primary advantage of semi-batch operation is that it allows the control of the driving forces of nucleation and particle growth (supersaturation ratio and temperature) during co-precipitation. Thus, to optimize the semi-batch operation, a significant decrease in supersaturation from a high initial value to a low value was selected, as shown in Fig. 14(a). This enabled the attainment of a large number of particles during the early operation time and restricted particle growth after reaching a certain level of metal recovery. As a result, at operation times below 18 min, the nucleation and growth rates were higher than those of the batch mode but decreased thereafter. Fig. 15 shows the maximum particle size and metal recovery across different operation modes. Because of its high nucleation and growth rates, the semi-batch mode achieved 80% metal recovery within 46 min, as shown in Fig. 15(b). In contrast, the operation time in the batch mode was 30.43% higher than that in the semi-batch mode. The lower nucleation and growth rates employed in the batch operation achieved a maximum particle size of 24.53 µm, which is 35.36% lower than that of the semi-batch operation, as shown in Fig. 15(a). Although the current study was conducted with limited computational resources (1000 iterations of BO), the comparison results clearly demonstrate that the semi-batch operation offers advantages in terms of saving operation time. When the maximum particle size is larger than that of the base case, further optimization settings or algorithms can be considered to effectively reduce the maximum size. Potential approaches include reducing the discretized interval of the operating trajectories and employing multi-objective BO.

Figure 14 Comparison of the PBE dynamics of the batch and semi-batch operations: (a) supersaturation ratio, (b) nucleation rate and (c) maximum growth rate.

Figure 15 Performance comparison of the batch and semi-batch operation: (a) maximum particle size and (b) recovery of the NCM precursor metal.

5. Conclusions

This study proposes an operational optimization process for co-precipitation during the production of NCM precursor from spent LIB cathode materials. By employing a data-driven classifier, 1D PBE simulations were conducted within an impurity-free operation domain. A PBE model was integrated with a data-driven classifier to enable the optimization of operating trajectories and minimize both operation time and particle size while ensuring 100% purity of the produced NCM precursor.

Unlike conventional co-precipitation methods, utilizing leaching solutions from spent LIB requires consideration of the complex equilibrium between the source metal ions (Mn²⁺, Co²⁺ and Ni²⁺) and impurities (e.g. Na⁺, SO₄²⁻, Al³⁺ and Cu²⁺) during PSD evolution. To address this issue, a thermodynamic equilibrium model was integrated into the numerical solver of the 1D PBE model. This integrated model accounts for the potential side reactions in the calculation of the supersaturation ratio, which is the driving force for nucleation and particle growth. Consequently, the proposed model can effectively handle the impact of impurities during PSD evolution, which is crucial for determining the purity of the produced precursor. The infeasible operating conditions, which lead to impurities [CoSO₄·6H₂O, Na₂SO₄ and (NH₄)₂SO₄], were identified using the thermodynamic model of the leaching solution. A DNN model served as the operating condition classifier to screen the infeasible domains with a high accuracy of 99.75%. This pre-screening reduces the number of operating conditions to be searched for minimizing both the operation time and maximum particle size of the final PSD. Under the optimal operating conditions, simulation results demonstrated that 100% purity of precursor could be theoretically achieved.

Two operating modes, namely batch and semi-batch modes, were optimized. In the batch mode, a maximum particle size of 24.53 µm with a metal recovery of 78% was achieved after 60 min of operation. Comparatively, employing the semi-batch mode saved 23.33% of operation time with a metal recovery of 80% but resulted in a 54.75% increase in the maximum particle size of the PSD. A comparison of the co-precipitation performance in different operation modes not only provided detailed operating guidelines but also assisted in the selection of the most suitable operation mode based on different objectives, including the reduction of the operation time and particle size of the product.

The current exploration of NCM precursor regeneration from spent LIBs is confined to experimental studies that lack the optimization of operating conditions or detailed investigations of performance, including impurity formation, operation time and particle distribution evolution. Addressing this critical gap is imperative to implement recycling approaches on an industrial scale. This study addresses this need by offering insights into the theoretical performance of the recycling process through modeling and optimization under different operation modes.

Despite its major contribution, this study has limitations that should be considered in future studies to improve the reliability of modeling and optimization. First, regarding PBE modeling, several kinetic factors, including multimodal distribution, agglomeration and breakage kinetic terms, can ensure the reliability of the simulation results. Additionally, integrating the PBE model with the kinetic Monte Carlo simulation would allow for the capture of more detailed microscopic kinetics (Kim, Pahari et al., 2024; Kim, Shah et al., 2024), particularly in relation to the multi-dimensional growth and morphology of both primary and secondary particles. This integration would enhance the reliability of proposed operation optimization, enabling operating conditions to more effectively meet target microscopic properties, such as particle shape and morphology, which are key factors influencing the electrochemical performance of the final cathode product (Koshika et al., 2022). However, expanding to a high-fidelity model without experimental studies can lead to issues regarding model convergence, overfitting and computational cost. Thus, developing a high-fidelity model should be accompanied by experimental studies to screen the model components on the basis of the real-world behavior of the NCM regeneration process. Furthermore, implementing surrogate modeling can potentially address the computational cost associated with a high-fidelity model during operational optimization (Choi et al., 2023; Zheng, Zhao et al., 2022; Zheng, Wang et al., 2022; Wu et al., 2023). Specifically, incorporating physical information into the model, such as using the physics-informed neural network approach, can embed the fundamental physical laws of crystallization directly into the surrogate model (Wu et al., 2023). This method effectively reduces the risk of overfitting by ensuring that the model respects known physical behavior. Moreover, it enables more reliable predictions of PSD and other detailed particle properties, even when computational resources are limited. Second, regarding optimization, this study discretized the operating trajectory and formulated single objective functions owing to limited computation resources. Supporting profile generation algorithms (Choong & Smith, 2004) and employing multi-objective algorithms, such as expected hypervolume improvement (Li et al., 2018), the non-sorting genetic algorithm (Deb et al., 2002; Joo et al., 2024) and the particle swarm algorithm (Hong et al., 2022; Hong et al., 2023), can improve the optimization results, specifically in the semi-batch operation model. To ensure the reliability of the current study, the proposed model and optimization must be expanded in future studies.

6. Related literature

The following references are cited only in the supporting information for this article: Kazakov et al. (2012), Van't Hoff (1884) and Zhang et al. (2018).