Design of experiments applied to lithium-ion batteries: A literature review

The statistical design of experiments methodology (DoE) has been a valuable tool since its conception for the understanding of the relationship between factors and responses. Although it has been employed successfully in different research fields and industries for years, its application to the evaluation of lithium-ion batteries (LIBs) is just getting recognition. LIBs are one of the most promising technologies for a complete transition to sustainable energies, are the main technology behind electric vehicles and are fundamental for the continual development of portable electronic devices. This paper presents a critical literature review of the available DoE works applied to the manufacturing and characterisation of LIBs. An overview of DoE and the most important available designs are first presented, followed by a general introduction of the statistical analysis required for the interpretation of the results including regression models. Several aspects of the LIBs such as ageing, capacity, electrode formulation, active material synthesis, thermal design, charging and parameterisation are discussed based on the main objective of the respective DoE studies found in the literature. A case study is presented to visualise the practical application of DoE to the LIBs field. Perspectives and future outlook are given to highlight opportunities and potential areas of research in the application of traditional and modern designs to the LIB’s field. This critical review contributes to a better understanding of the DoE methodology with a focus on LIBs or LIBs related aspects which will lead to faster developments in the field.


Introduction
The search for more sustainable alternatives to fossil fuel energy resources and the information era has led to the development of lithiumion batteries (LIBs). LIBs are helping in the substitution of oil powered cars by electric vehicles (EVs) [1], and are aiding in the transition to renewable energy sources by serving as energy storage devices [2]. Additionally, the demand for portable electronic devices such as mobile phones, laptops, tablets, digital cameras, wearables and drones, has resulted in a considerable growth in LIBs production [3]. However, battery performance and cost still need to be improved to facilitate a complete transition towards electrification.
LIBs performance is mainly described by their energy density, power density, capacity, lifetime and charging time, characteristics that are a function of electrochemical properties (e.g. current and ion density, coulombic efficiency, ohmic resistance, polarization and reaction losses) and physical electrode properties (e.g. thickness, mass loading, porosity, adhesion, conductivity) [1]. The electrochemical and physical properties, and finally, the performance of the battery is the result of the interaction of battery components (anode, cathode, binder, separator and electrolyte solution) [1] and the electrodes structures created in the several manufacturing process steps [4]. In order to optimise the performance, quality and cost of the manufactured LIBs it is first necessary to understand how the different components and manufacturing processes interact and their impact on battery characteristics.
The understanding of the relationship between the manipulated variables and the output variables is essentially done by experimentation, frequently by trial-and-error through a best-guess approach, or in a few cases by a one-factor-at-a-time approach. In the context of LIBs manufacture and performance, the manipulated variables are, for instance, the types and amounts of materials in the electrode formulation, the thickness and porosity of the electrode, charge and discharge temperatures, and equipment operating parameters, whereas the output variables can include the battery's energy density, discharge capacity and temperature rise.
In recent years, the combination of experiments and modelling has shown to be a promising alternative to only experimental work [5]. Some researchers have focused on reducing the number of experiments required to understand the relationship between battery performance and the manufacturing process by using models at different scales [6,7]. Schmidt et al. [4], for instance, combined a process chain simulation consisting of the coating, drying and calendering processes, and a homogenized cell model to identify the most relevant process parameters and their impact on battery performance.
Nevertheless, accurate experimental data underpins all modelling work for either its development, testing, validation or parameterisation. To reduce the time and efforts associated with experimentation, the statistical Design of Experiments (DoE), is a valuable tool to obtain the maximum amount of relevant data at the minimum economic cost and time. DoE has successfully been used in different areas and industries, such as the pharmaceutical [8,9], agricultural [10,11], energy and bioenergy [12], fuel cells [13], microencapsulation [14], analytical chemistry [15] and chemical and biochemical processes [12,16]. In the battery field, however, DoE seems to be just gaining recognition as shown in the number of published papers through time (Fig. 1).
This paper presents a review of DoE works applied to LIBs within the academic literature. The paper is organised as follows: first, an introduction to DoE is provided, covering the aspects of design methodology, available experimental designs and the statistical tools used for the analysis of the collected data. The article then presents and discusses a collection of works where DoE has been applied specifically to the LIBs field. The works have been classified in the following categories depending on the main objective of the DoE: battery ageing, energy capacity, formulation, active material synthesis, electrode and cell production, thermal design, charging, other applications, optimisation studies and model parameterisation. A DoE case study is presented to illustrate a practical example of the application of DoE and how the experimental effort is reduced compared when traditional methods are employed. Perspectives and future outlook are given highlighting potential areas of research. General conclusions are presented at the end. The paper contributes to a better understanding of the DoE methodology applied to the LIBs field and clarifies a few of its misconceptions.

DoE methodology
DoE can be defined as the branch of statistics involved in the planning, the collection and analysis of experimental data to ensure valid and objective engineering conclusions are attained [17]. DoE deals with the understanding of the effect of independent or input variables (factors) on the dependent or output variable (response) [11]. DoE is built on concepts of randomisation, blocking, replication, factorial approach and analysis of variance first introduced by Fisher [10,18]. The study areas in which DoE is often applied in engineering are comparative, screening, modelling, optimisation, robust design and formulation ( Fig. 2) [11,19].
DoE involves a series of methodical steps (Fig. 3) [20]. First, the problem as well as the objectives of the experimental study are stated. The problem can be classified in one of the study areas shown in Fig. 2. The response(s), as well as the factors and their suitable levels (settings of the process parameters) are chosen based on the objective of the experimental study. Several factors (>4) at only 2 levels are normally considered for screening studies, whereas 2 or 3 factors with 3 or more levels are common for optimisation or robust parameter design. Next, an experimental design is selected, which is the category comprising the series of experiments that are undertaken considering the number of factors, levels, replicates, blocks, randomisation and the consideration of an (empirical) model. The next step is performing the experiments as dictated by the experimental design. In most cases, particularly in screening, it is useful to run a few preliminary experiments to ensure that a suitable range of the factors has been chosen.
Once the data has been collected, it is processed using statistical methods. The statistical analysis commonly involves the analysis of variance (ANOVA), as well as the use of graphical methods (e.g., Pareto charts, histograms, normal or half normal probability plots and mean plots). In this step, a model, usually empirical, can be used to interpret the results and to represent the relationship between the factors and the response (Section 4.1). Finally, the evidence given by the statistical analysis provides a basis to arrive to objective conclusions.

Experimental designs
It is beyond the scope of this paper to discuss in detail the different experimental design methods. Further details are provided in a number of publications and educational texts, for example [11,21,22]. However, for completeness, a summary of the main experimental designs is provided below.

Factorial
This type of designs involves either the testing of all possible combinations of two or more factors and their levels (full factorial) or the testing of only a subset of the full factorial (fractional factorial design). The number of experiments (N) for a full factorial is given by N = L k , where L is the number of levels and k the number of factors. For a fractional factorial it is N = L k− p , where p generates the confounding pattern of the design [17].
Full factorials are rarely used when the number of factors or levels is relatively large (L > 3, k > 5). Common factorial designs are the 2 k and the 3 k . The former is mainly employed in screening (see Section 3.2), whereas the latter when a quadratic function between response and factors is suspected. Nevertheless, 3 k are common only when the number of factors is relatively low, and if quadratic or higher interactions are expected, response surface methodology (RSM) provides with better designs.
An example of a factorial DoE applied to LIBs is given by Rangappa and Rajoo [23] (see Table 7) who used a full 3 3 design for the identification of the main factors in battery temperature rise.

Screening
As previously mentioned, 2 k designs are primarily used for screening purposes. A special type of 2 k fractional factorials are the saturated designs in which up to k = N − 1 factors can be studied in N runs. An example of a saturated design used for screening is the Plackett-Burman in which N is a multiple of 4 [24]. Despite their effectiveness in finding main factors, the use of Plackett-Burman designs has not been reported in the LIBs field.
Screening designs have been mentioned for the study of vehicle-togrid (V2G) and grid-to-vehicle (G2V) strategies on battery degradation [25], and for the analysis of impurities on lithium extraction (2 4 design) [26]. A limited collection of papers reported the use of screening designs for the identification of main factors as a preliminary research step. Some works identified the main factors by orthogonal arrays (OAs) with more than two levels (e.g., [27][28][29][30]) which may not be economical.

Response surface methodology
After the work of Fisher [10], a wider application of DoE in industry came with the RSM proposed by Box and Wilson [18,31,32]. The methodology focuses primarily on process optimisation but is also used for robust product and process design. RSM acquires its name from the graphical representation of the surface created from the values of the factors and the resulting value of the response. For the case of two variables, the surfaces are a 3D plot with its associated contour plot (see for instance Fig. 6). RSM studies are recommended after screening designs in order to perform targeted experiments only on the main factors.
The relationship between factors and response is expressed by a mathematical model, normally a second-order polynomial (Section 4). If the model is deemed to represent the relationship between factors and levels, it can then be used to determine the settings to optimise the response.
Although there are different RSM designs (central composite, Box-Behnken, the small composite design, equiradial designs, hybrid design) [11], the most popular and widely used in industry are the Box-Wilson central composite design (CCD) and the Box-Behnken design (BBD) [31,33].
In CCD, a factorial constitutes the basis of the design to which axial runs and centre runs are added resulting in factors with five levels. The number of centre runs (c), are required to provide stability in the variance of the predicted response. The total number of experiments for a CCD is N = k 2 + 2k + c.
The BBD, conversely, consists of three evenly spaced levels and the number of experiments is N = 2k(k − 1) + c, for k > 2. The designs are formed by the combination of 2 k factorials with incomplete block designs [32].
A detailed application of RSM to electrode formulation is discussed by Lv et al. [34] (see Table 4).

Taguchi
The popularization of DoE in industry came with Taguchi's work on what was termed robust parameter design (RPD). RPD is a methodology and philosophy to determine the best settings of the controllable factors to make the process insensitive to noise factors in order to reduce process or product variation [11,32]. The product or process is, in this sense, optimised. Noise factors are defined as the variables that cannot be controlled in the actual process but that can be controlled or simulated in the experiments. Taguchi promoted the use of fractional factorial designs in the form of OAs (a class of orthogonal main effects designs) consisting of a N × k matrix [35]. The nomenclature used to define the Taguchi designs is L N ( L k ) , or simply L N . In Taguchi designs an OA of the control variables (factors) is crossed with an OA containing the noise variables. The signal-to-noise ratio (S/N) is used as a measure of the variability of the response with respect to the target value due to the noise factors. Three common S/N relationships can be specified depending on the objective of the RPD experiment: nominal-the-best, larger-the-better and smaller-the-better (Table 1).
Nevertheless, it has been shown that Taguchi designs can lead to a large number of experiments, are in most cases unable to identify parameter interactions, and moreover, the use of the S/N as the response does not allow to distinguish between the factors that affect the mean and the factors that affect the variance [32,36].
Most of the studies found in this review, discussed further in Section 5, are of the Taguchi type.

Mixture
Mixture designs are a special type of designs in which the factors (x) are the compounds in a mixture and the levels their proportions, related according to Eqs. (1) and (2), where q is the number of components. [37]: In mixture designs it is commonly assumed that the response is only a function of the proportions and not of the amount of the total mixture. Graphically, for two components, the factor space is a line, for three components is a triangle and for four components it is a tetrahedron. The common types of mixture designs are the simplex lattice, simplex centroid, simplex axial and extreme vertices design [21,38].
In the simplex lattice, there are m +1 equally spaced proportions (Eq. (3)) in the range from 0 to 1, and all mixtures with these proportions for each compound are evaluated (see for instance Fig. 4a). These designs  Table 1 Typical signal-to-noise ratio (S/N) relationships and their applicability. a Signal-to-noise ratio type Objective Equation

Nominal-thebest
To achieve a specific value of the response, e.g., electrode thickness Larger-thebetter Maximise the response, e.g., coating adhesion
In the simplex centroid design, the experimental points are in the centroid points, the vertices and the overall centroid (centre point of the triangle, for q = 3). The total number of experimental points is 2 q − 1 (Fig. 4b).
In simplex axial designs, the experimental points are on the axial of the components. The axial is defined as the line extending from the base point x i = 0, x j = 1 q− 1 , to the vertex x i = 1, x j = 0, for all j ∕ = i [21]. The total number of points in a simplex axial is given by 3q + 1, corresponding to q vertex points, q centroid of constraint planes, q axial points and the overall centroid.
In the extreme vertices design, a reduced design space of the simplex is created due to linear constraints or upper and lower bounds on the proportions. The number of runs depends on whether axial points and overall centroid point are added.
Recently, Rynne et al. [39] optimised a specific electrode formulation by a constrained mixture designed [21,40]. The research consisted in testing the proportions of active material (AM), conductive additives (CA) and polymer binder (B), and their effect on discharge/charge capacities. The analysis showed that the AM should be maximised while keeping a small fraction of CAs and minimum B content.

Optimal designs
The previous designs are also called traditional or classic designs [17,19]. Optimum designs, on the other hand, are an umbrella term for designs created based on specific objective optimality criteria. Because of their complexity, computer algorithms are often needed for the implementation of optimum designs, and are therefore also termed as computer-generated designs [41,42]. There are various types of optimal designs, but frequent ones are D-optimal, G-optimal, I-optimal, Aoptimal and V-optimal, each with different features (Fig. 5) [43,44].
Optimal designs are primarily used when dealing with non-standard models, a restricted number of experiments or restricted experimental  regions [32,45]. Few papers reported optimal designs, and these were mainly D-optimal. Only Rynne et al. [39] reported I-optimal in their formulation study.

Regression model
One of the goals of DoE is to obtain a mathematical model that represents the relationship between the response and the factors. For a process represented by a k number of x factors, the response y is given by Eq. (4) where ε is the random error.
If the expected response, E(y), is a function of the x k factors (Eq. (5)), then the surface represented by Eq. (6) is the response surface [11].
In DoE, an empirical model in the form of a polynomial is used to fit the experimental data. The simplest model to represent the response is the linear function (Eq. (7)), also called the main effects model since it only contains the effects of the x i factors. β 0 is the average value of the response.
More complex relationships can be studied by including an interaction term to the main effects model. For instance, a model with two-way interactions (x i x j(i∕ =j) ) can be expressed as in Eq. (8).
If curvature is shown or expected in the system, then higher order polynomials may be used. The higher order-order polynomial can also include higher interaction products (e.g., three-way) but for most DoE studies, normally up to the second-order model is considered (Eq. (9)).
The coefficients β i , β ij and β ii in Eqs. (7)-(9) are the first order, the interaction and the quadratic coefficients, respectively. The β parameters are estimated from fitting of the experimental data by the method of least-squares.
Eqs. (7)-(9) can be expressed in matrix notation according to: A special matrix resulting from a given experimental design and model is the hat matrix (H) defined as: The fitted regression model is given by: where ŷ is the least-squares estimate of E(y), or: θ in Eq. (13) is the least-square estimator of β defined by Eq. (14) according to the least-square analysis [32,45].
The difference between the observed values (y i ) and the fitted values (ŷ i ) is the residual, e i . The vector of residuals is therefore: The covariance matrix of θ is defined as: where X ' X is called the design information matrix. It is noteworthy that for mixture designs, the polynomials representing the response surface are different due to the constraint given by Eq. (2). The second-order canonical mixture (or second-order Scheffé model) and the special cubic models are given by Eqs. (17)-(18) [21,44]. Other models can be found in the review of Cornell [21].

Model diagnostics
The first step in the model diagnostics is to look at residual plots in the form of residual vs predictions, residual histogram, normal probability plot of residuals, among others [17,46]. Numerical indicators for  [34]; (b) a generic rising ridge system. model testing basically involve ANOVA, that is an analysis of the total variability in the response. The analysis starts with the computation of the regression sum of squares (SS R ), the residual sum of squares (SS E ) and their contribution to the total sum of squares (SS T ) according to Eqs. ((19)- (21)).
The regression (MS R ) and residual (MS E ) mean squares are then computed according to Eqs. (22) and (23), respectively, which are in turn used to determine the F-value (F-statistics test, Eq. (24)).
If F-value > F ∝,k,n− k− 1 , then the null hypothesis (H 0 : β 1 = β 2 = … = β k = 0) is rejected and the model equals noise. Values of F ∝,k,n− k− 1 can be obtained from tables in statistics books, where ∝ is the level of statistical significance.
It is also customary to use P-value criteria to reject the null hypothesis when the P-value < ∝. Calculation methods of the P-value are not straightforward, and it is normally done by statistical software (e.g., SAS, Design Expert, Minitab, R, MATLAB and JMP). ∝ is generally set between 0.05 and 0.1; the latter value is preferred in screening designs.
The coefficient of determination (R 2 ) computed from Eq. (25) is commonly used to explain the data variation that the regressed model can explain is. The closer R 2 to 1 the better the model is in representing the data. However, because the addition of terms in the model increases R 2 regardless of their statistical significance, the adjusted R 2 (R 2 adj ) (Eq. (26)), is a preferred representation of the data variation explained by the model. Unnecessary terms will decrease the value of R 2 adj .
A proposed estimate of the predictive capability of a model is given by R 2 pred , according to Eqs. (27) and (28).
h ii in Eq. (28) are the diagonal elements of the hat matrix (Eq. (11)). The method behind Eq. (27) is to remove one of the i data points, fit the regression model to the remaining n − 1 observations and use that model to predict the removed y i point. The process is repeated for all data points to compute PRESS.

DoE applied to batteries
The following sections present the application of DoE by subject research area depending on the main objective of the study.

Battery ageing
Ageing within a LIB is classified as cycle ageing and calendar ageing; the former is related to degradation processes linked to mechanisms to the charge-discharge cycle, whereas the latter to independent mechanisms of such cycles. Calendar ageing has been mainly attributed to the build-up of the solid-electrolyte-interphase (SEI) [47]. Cycling ageing is attributed to SEI growth, lithium platting in the anode, volume changes and material degradation in both electrodes [48]. These mechanisms are influenced by several factors like cell chemistry, cell temperature, state of charge (SoC), current magnitude, depth of discharge (DoD), frequency of the charge-discharge process [49,50], voltage (charging voltage, end-of-discharge voltage and end-of-charge voltage) [51], pulse duration [52] and superimposed AC current [53,54].
Ageing is mainly defined as capacity fade and power fade (increase in the internal resistance or impedance) [50]. The understanding of the relationship between the different ageing factors on these two performance indicators (responses) is important for both, the design of more durable LIBs and for the management of existing batteries through battery management systems (BMS). The relationship between ageing factors and the response is expressed mathematically by either a physico-chemical model or an empirical or semi-empirical model. Nevertheless, the several factors involved in the aging process make the physico-chemical modelling a challenging task [55], and the simple empirical and semi-empirical models are sometimes preferred. The calibration of these models (parameterisation) is done by fitting of experimental data.
Performing the minimum number of experiments is important in ageing tests since these can take from months to years, especially when dealing with new chemistries or in the development of new battery designs [48,56]. DoE is useful not only to determine the minimum and most valuable number of ageing tests conditions but also to obtain the parameters to underpin empirical model development. A detailed review of literature in this domain is summarised in Table 2.
Although different chemistries have been involved in the ageing studies (Table 2), generally, temperature (T) and SoC have been identified as the two main important factors in capacity fading. The effect of SoC is mainly due to interaction with temperature. While the effect of temperature and SoC was already determined by studies not using DoE [47], the usefulness of the DoE methodology is in arriving to such conclusions faster, as well as to obtain an empirical model to explain the observations.
Using a D-optimal design, Mathieu et al. [48] obtained cycling degradation data that was later coupled with independently tested calendar data to develop an integrated empirical model for both ageing types. The resulting model was a double quadratic expression with interaction (Eq. (29)). The response surfaces at constant charge current (I C ) and discharge current (I D ) are saddle points that can be used as an indication of the behaviour for capacity fade rate under the studied factors ranges. The model revealed the capacity fades plateaus corresponding to the graphite potential plateaus when SEI growth is the predominant ageing mechanism as pointed out by Keil et al. [47].

Energy capacity
Energy capacity is one of the LIB's key performance indicators and an active area of research. The required capacity of a LIB depends on its final application (e.g. portable electronic, EVs, storage unit) [63]. As shown in Table 3, only a few papers have used DoE to study the effect of electrode physical properties (e.g. thickness, volume fraction, porosity) (continued on next page) on energy density and power density. Gitzendanner et al. [64] employed a fractional factorial design for the identification of the main factors affecting the discharge capacity in the design of a 20 Ah prismatic cell for military and aerospace applications. The only statistically significant factor (∝ < 0.05) was the cathode/anode weight loading ratio. From the statistical analysis and engineering design considerations, suitable values of the factors were chosen resulting in cells exceeding the design goals, for instance higher cycle life. The work of Hosseinzadeh et al. [65] used DoE to assess the effect of individual parameters and parameter interactions on specific energy and specific power. A 1D electrochemical-thermal model previously developed [66] was first implemented in COMSOL Multiphysics and validated against literature experimental data. The model was then used to determine the values of the specific energy and specific power at different parameter settings according to a full factorial design to identify and ranked the main parameters based on a first order polynomial with interactions. The empirical model was used for optimisation of the responses based on the main factors and significant interactions. The contour plots, however, do not reveal optima for any of the responses, exhibiting instead a rising ridge for specific energy and a twisted plane for specific power, due to the interaction terms. Although 81 computer experiments (simulations) were done, corresponding to a 3 4 design, the same information could have been obtained in 27 simulations through a BBD, or in as few as 13 with a Definitive Screening design (DSD) [67,68].

Formulation
Electrode formulation comprises a mixture of components, mainly AM, CA and B. At the cell level, an electrolyte solution creates a medium for the movement of Li-ions between the electrodes. The proportions of these components have a direct effect on battery energy density and power density and will vary depending on the final application [70,71].
Factorial and RSM have been employed in the screening of several factors, including cathode and anode materials, solvent type and doping materials (Table 4). This is a valid approach if the response is not studied as a function of the relative proportions of the ingredients but rather as the total amount of the components. For the cases where the response varies as a function of the relative proportions, a mixture DoE should be used instead [21]. Lv et al. [34] used a CCD to determine the amounts of transition metals as co-dopers in LFP cathodes to optimise discharge capacities and capacity retention. The results are in this case a function of the amounts of the metals as well as the total amount of mixture. A mixture design could have been used to fix the total mixture amount varying only the proportions of the metals, this approach would have led to the determination of an optimum composition.
The first applications of mixture designs to the LIBs field seems to be the work of Park et al. [82] ( Table 4). The authors developed a tolerance design method to deal with mixture errors by using a crossed design of a fractional factorial with an extreme vertices design. The proposed method is applicable to any robust parameter design mixture problem. As an example, Park et al. [82] applied it to the tolerance problem in preparing a three-compound carbonate solvent mixture for electrolyte formulation. More recently, Rynne et al. [39] determined the optimum formulation in terms of charge capacities for varying compositions of active material (LFP and LTO), additive (carbon black and nanofibers) and binder (PVDF, TPE and HNBR). As pointed out by the authors, effects due to solvent type and rheology may have a bigger influenced on the response and these factors need to be carefully controlled in this type of studies. Studies combining a formulation problem (mixture design) with process factors have not been found in the literature. These type of studies are important, for instance, to determine the best electrode formulation as well as the optimum coating and drying conditions, to maximise energy density.

Active material synthesis
DoE can help in faster determinations of the influential parameters and optimum conditions for electrode material synthesis. However, only a few publications were identified employing DoE to this area (Table 5). Only one paper truly identified optimum conditions by RSM [83]. The rest of the papers determined the optimum conditions by Taguchi orthogonal designs. Taguchi designs, however, cannot identify the second order terms necessary to determine curvature, and can only identify interactions in a few cases [32,36]. The OA, as employed in the works in Table 5, helped only to the identification of the main factors affecting the synthesis process. The conditions reported as optimum are only the factor settings that resulted in the highest or lowest value of the response, but these are not necessarily the optimum.

Electrode and cell production
The quality, and therefore performance of the battery not only depends on the components of the electrodes, but also of the production processes at the electrode and cell level. Studies on process uncertainties and their effect on battery performance are scarce due to the considerable amount of time required in manufacturing the cells and the electrochemical measurements [4].
Several processes are involved in the electrode and cell production including mixing, coating, drying, calendering, welding, winding (or stacking or Z-folding), cutting, contacting, housing and electrolyte filling [1,91]. Of the many processes, however, only a few (welding, winding, coating and calendering) have been studied using DoE methods ( Table 6).
Contributions of DoE applied to cell manufacturing are presented by Meyer et al. [92] and Billot et al. [93]. The former optimised the web tensions for the electrode and the separator in the jelly roll winding process for plug-in hybrids EVs. A CCD identified the optimum settings for the web tensions for six different responses. Since no general optimum was found for all responses, desirability indices determined the optimal levels of the factors. Of the six responses studied, only one (spring back after 25 h) resulted in R 2 > 0.9, probably due to the omission of other influential factors in the model [92].
Billot et al. [93] presented a comprehensive study of the adhesion strength (a parameter linked to LIB's capacity) from four different  adhesion theories: mechanical, thermodynamic, chemical and micromechanical. Influencing key factors were chosen from each theory (mechanical: roughness of the foil; thermodynamic: surface tension of the foil and of the slurry; chemical: binder amount; micromechanical: coating thickness) through a series of DoE plans. The main factors identified were the amount of binder and the coating thickness. In a different DoE, the compression rate, porosity and roll temperature were analysed. All main factors were identified as highly statistically significant, as well as the porosity-temperature interaction. Contrasting with the findings in the literature where observations are normally done at room temperature, making the adhesion strength less noticeable, the results showed an increase in adhesion strength at low porosity and high temperature [93]. In another DoE, the effect of roll diameter, layer thickness, roll temperature, number of passes, porosity and web speed were studied. Results showed that smaller roll diameters result in lower adhesion strengths. Unfortunately, not many details of the DoE plans were presented, and only the use of a D-optimal design is reported in the study of the roll diameters.

Thermal design
Temperature is an important aspect of battery performance and safety, manifested in the forms of operating temperature, storage temperature and battery temperature distribution. In the mildest case, temperature contributes to the calendar and cycling mechanisms, and its effects are seen as battery degradation [63] (Section 5.1). In extreme cases, a poor battery temperature distribution and localised hot or cold spots could lead to accelerated battery degradation and in thermal runaway [96].
DoE plays an important role in battery design to study the effects of potential design parameters (materials, dimensions, coolant flowrates, added features, etc.) that guarantee an even battery temperature distribution, minimise temperature rise and released energy. These aspects in turn reduce the possibility of unexpected thermal runaway, or in the worst case, minimise its catastrophic effects.
The DoE types in thermal design found in the literature varied from fractional factorial designs to BBD and even full factorials ( Table 7).
Most of the studies have determined a linear or quadratic relationship between factors and responses. A comprehensive DoE study by Walker et al. [97] have shown that inclusion of a bottom vent in the battery reduces the severity of a thermal runaway and improves the prediction of the event. Design features studied also included battery capacity, manufacturer separator material, cell casing thickness and an internal short-circuiting device. on the responses were energy distribution, total energy released and mass distribution. A bespoke calorimeter was developed to examine the effect of thermal runaway. A lognormal distribution model explained the relationship between the factors and the response.

Battery charging
An interesting application of DoE in LIBs is in finding the optimum charging pattern, or optimum settings of charging parameters, that will result in faster charging times, higher charging capacities or increased cycle life. Two charging methods that have been proposed as alternatives to the simple constant current (CC), constant voltage (CV) and the more commercially deployed CC-CV methods, are the pulse charging (PS) and the multi-stage CC (MS-CC). In PS, the intermittent voltage or current used to charge the battery reduces dendrite formation by avoiding accumulation of the Li + ions on the electrode surface, resulting in larger charging currents and lower charging times at the fixed cut-off voltage [98]. The optimisable pulse charging parameters are frequency, peak charge current amplitude and duty cycle [99]. In MS-CC, the total charging time is divided into stages, each set to a predetermined charging current. When the cut-off voltage is reached during charging at the present current, charging moves to a new stage under a new current until the cut-off voltage is reached again. The process is repeated until the final stage number is reached (normally between 4 and 5) [100]. The challenge in MS-CC is to find the optimum current value for each stage. Although the optimum can be found through optimization algorithms (including genetic algorithms, particle swarm optimisation, ant colony system), these rely on accurate battery models which are difficult to obtain due to the varied and complex electrochemical LIBs' properties [101,102]. DoE has been applied as a substitute tool for finding the optimum settings of the PS method and for finding the optimum current for each of the MS-CC, avoiding the reliance on a battery model and on the timeconsuming testing of all the current candidates. In particular, OAs as promoted in the Taguchi method have been the preferred DoE type in these kind of studies (Table 8).
Liu and Luo [103], and in more detailed Liu et al. [101] proposed the use of Taguchi designs to determine optimum charging patterns for the MS-CC method. A few variations of the proposed methodology have been developed (Table 8), but the central idea is that through an iterative process, the best charging current combination for each stage is obtained by finding the best cost function from the charging conditions determined by an initial OA and an initial set of current candidates. If the termination criteria are satisfied at this iteration, the optimum charging pattern is that given by the best cost function; if not, a new iteration with updated current candidates is tested by creating a new OA. The process is repeated until the termination criteria is satisfied.
Despite the effectiveness in the use of OA for finding the best conditions in MS-CC methods, the number of elements in the testing matrix (design matrix) can be reduced by first creating the factorial design and then applying optimal design criteria taking advantage that no factor interactions are needed, i.e., only the main effects of each charging stage are needed. For instance, a design based on D-optimality would comprise 11 experiments instead of 16 as in the L 16 array, for a main factors model. Optimum pulse charging parameters are also better obtained by RSM instead of OA, since the latter do not allow in most of the cases for the study of two-way interactions, or second order terms as mentioned in Section 3.4.

Other applications
DoE is also valuable in studying aspects surrounding LIBs, which when integrated to the whole production chain, can result in lower battery development costs, lower energy consumption or improved performance. The primary examples are presented in Table 9, and include research on improved synthesis conditions for electrode materials, a study on the effects of impurities in lithium extraction and a study on corrosion processes. Notably, all the studies involve full or fractional factorial designs.  [106] a N/A: Not available.

Optimisation studies
One of the primary applications of DoE is for process optimisation. This step is normally carried out focused on two or three main factors identified through a screening design. When the relationship between factors and response can be approximated by a second order or higher polynomial (Eq. (9)), it is possible to find the stationary points of the equation corresponding to minimum, maximum or saddle points. The contour plots for response surfaces exhibiting maximum or minimum points are concentric ellipses or circles (Fig. 6a). For surfaces with a saddle point, the contour plots correspond to hyperbolas. When the response surface is a rising or falling ridge (Fig. 6b), it is an indication of a system where the optimum is outside the studied experimental region [11]. In this case, appropriate adjustments in the factor ranges are required to shift the region towards the optimum.
Of the sixteen papers found discussing optimisation (Table 10), only six used RSM as the proper method for optimisation (one is a mixture problem), the rest employed OAs or full factorial (one case). The works presented as optimisation by OAs were used as part of the Taguchi methodology. Perhaps most of the confusion in these works arise from the misinterpretation of the S/N relationships where larger-the-better, smaller-the-better or nominal-the-best are confused with optimum values. Taguchi methods were developed to make processes robust to noise factors and primarily to the identification of main effects, and it is applied to find optimum settings in terms of best values for minimum variation. Most of the optimisation results reported under the Taguchi method are therefore the identification of main factor effects rather than optimum values since no stationary point is being reported.
The OA and the analysis of means can be used for the identification of main factors and the optimisation carried out by other methods. For instance, Lee et al. [69] identified the four main factors from an initial pool of eight, which were then considered design variables for the specific energy optimisation based on Newman's electrochemical model and the micro genetic global optimisation algorithm.
Only a few papers presented response surface plots or contour plots, and where mainly ridge systems or bended planes. The surface reported by Lv et al. [34], on the other hand, clearly shows an optimum in the study of doping quantities of vanadium and titanium in LFP cathodes (Fig. 6a). The authors also reported the statistical analysis of the results (R 2 , ANOVA, F-value, P-value), something not generally reported in the reviewed works despite being the statistical analysis a core principle of DoE.

Model parameterisation
The design and efficient management of LIBs relies on accurate models to describe their behaviour at different operating conditions and for different chemistries. Such models are used, for instance, in BMS controlling voltage, current, SoC, state of health (SoH) and temperature [121].
Available models can be classified into three main categories: empirical models, electrochemical models (EMs) and the equivalent circuit model (ECM). The empirical models are the simplest of all but also have the lowest accuracy [51]. EMs involve the modelling of physico-chemical phenomena, are more complex, but are also more accurate [48]. ECMs are a trade-off between the empirical and the EMs, are formed of basic electrical components (voltage sources, resistances, capacitances) and in some instances of nonlinear elements [122][123][124]. A detailed EM is the Doyle-Fuller-Newman (DFN) model, also known as the pseudo-two-dimensional (P2D) model [123,125]. However, their complexity has hindered its widespread application and simpler EMs have been proposed, like the single particle model (SPM) or SPM with electrolyte dynamics (SPMe) [123,126]. The effects of temperature on battery behaviour are also important for accurate modelling, and require coupling of electrical-thermal models [127].
Regardless of the type of model, parameter estimation is an inherent Charging time and discharge capacity ratio [107] Implementation of an SOC-based four-stage constant current charger for Li-ion batteries Battery cycle life, series resistance and charge transfer resistance [109] Search for optimal pulse charging parameters for li-ion polymer batteries using Taguchi orthogonal arrays Charging time and normalized discharge capacity [112] New charging strategy for lithium-ion batteries based on the integration of Taguchi method and state of charge estimation Energy efficiency charging time and T variation [113] Obtaining optimal membership functions using fuzzy-based Taguchi method N/A L18(2 1 × 3 7 ) 5 charging stages (3 for each) Charging time, charge capacity and T [114] Optimization of a fuzzy-logic-control-based five-stage battery charger using a fuzzy-based Taguchi method N/A L18(2 1 × 3 7 ) 5 charging stages (3 for each) Charging time, charging efficiency and T variation. [115] The implementation of consecutive orthogonal array method on searching optimal five step charging pattern for lithium-ion batteries LMO L18(2 1 × 3 7 ) 5 charging stages (3 for each) Charge capacity [116] Search for an optimal five-step charging pattern for li-ion batteries using consecutive orthogonal arrays LMO L18(2 1 × 3 7 ) 5 charging stages (3 for each) Charge capacity [101] Search for an optimal rapid-charging pattern for li-ion batteries using the Taguchi approach LMO L18(2 1 × 3 7 ) 5 charging stages (3 for each) Charge capacity and charging time [103] a N/A: Not available.
part of the modelling process. While some parameters can be measured or predefined, others are obtained from experimental data. The methods and complexity of the parameter estimation from experimental data will vary depending on the model use but are in general more complicated for the ECM and EM.
Here, we highlight the use of DoE for parameter identification based on two approaches. The first is the use of DoE to reduce the number of experiments needed to construct the response surface from which the parameters are obtained. The second, DoE as a tool for parameter estimation and identifiability, that is, determining whether these can be obtained from experimental data [128]. An example of the former is presented by Mathew et al. [127] (Table 11), where a CCD constructed the response surface to get the parameters of an electro-thermal model. A comparison of the use of DoE with Non-DoE methods and look-up tables revealed that although the same accuracy was obtained from the three approaches, the DoE reduced the experimental time by around 75%.
The number of experiments can also be reduced by specifying the number of experiments based on D-optimal design. This approach will also maximise the parameter estimation information, since as mentioned in Section 3.6, D-optimal designs maximise the determinant of the information matrix (Eq. (16)). In this way, Forman et al. [129] proposed a framework for optimal experimental design applied to battery experimental and modelling development. The framework was first demonstrated on an empirical cubic model from 14 experiments, and later in the identification of 88 DFN model parameters [130].
Since the unbiased least-square estimator satisfies the Crámer-Rao inequality [129], the best possible covariance in estimating θ is given by the inverse of Eq. (30), called the Fisher information matrix. Parameter estimation based on Eq. (30), are referred as Fisher-based optimal experimental designs. Several battery models have been studied and their parameterization obtained by the Fisher-based approach ( Table 11). The relationship between estimator performance and model structure (ill-conditioning) was also studied by López C et al. [131]. The use of the Fisher-based optimal experimental design has revealed that some of the parameters of the DFN model are unidentifiable. Taguchi optimization of the carbon anode for Li-ion battery from natural precursors

Carbon
To identify the best natural precursors and operating settings to obtain the best carbon from a pyrolysis process and acid and alkali treatment.

L9(3 4 ) Precursor, T, Pretreatment and gas. (3 for each). 2 repeats
Charge capacity and S/N S/N N/A N/A [120] a N/A: Not available.
6. DoE application case study

Experimental design
To illustrate the benefits of applying the DoE methodology to the LIBs field, the experimental data reported by Hamed et al. [138] on LiNi 1/3 Mn 1/3 Co 1/3 O 2 (NMC) electrodes is studied in the present work through a DoE approach. The authors were investigating the correlation between the electrodes' rate capability and their microstructural properties through a series of experiments and the P2D model. The research involved the study of three factors: dry thickness (δ = 38 µm, 76 µm and 108 µm), NMC loading (m NMC = 87 wt%, 92 wt% and 95 wt%) and binder/carbon weight ratio (B/C = 1, 1.25, 1.67 and 2.5), at 3, 3 and 4 levels, respectively. The examined responses were the tortuosity (τ), the rate capability coefficient (γ) and the fractional long-range energy loss (ω). A series of 36 experiments were performed which are the product of testing all the factors and level combinations, i.e., performing the full factorial experimental design. The study successfully demonstrated, by a series of 2D plots, the influence of the selected factors on tortuosity and the performance limitations of the electrodes.
As mentioned in Section 3.1, full factorials are rarely used for L > 3, since the number of experimental runs increases considerably. If the 3 factors would have been studied at 3 levels each, the required number of experiments for the full factorial would have resulted in, perhaps, a more manageable number of 27 experiments. Nevertheless, the selection of a 4th level for B/C may have been deemed necessary, or a requirement, for development purposes. An alternative to the full factorial design used in the original publication is presented in the current work for the same number of factors and levels. The proposed design is a RSM I-optimal design aimed at developing a predictive empirical model containing the statistically significant terms and to visualise the effect of the factors on the responses. The optimal design was created in Design-Expert [139] treating the factors as discrete variables to be able to use the reported data in the original publication. The minimum number of experiments to fit a quadratic model (Eq. (10)) for the given combination of factors is 10, but a total of 15 experiments were considered to obtain a robust model. Even with 15 experiments, the I-optimal design chosen reduces considerably the experimental effort by more than 50%. Failure to include enough experiments for parameter estimation would have required augmenting the design with extra runs. The design matrix for the RSM I-optimal design is presented in Table 12. The corresponding values of the responses were taken from the reported data by Hamed et al. [138] and are also presented in Table 12.

Tortuosity
Linear regression analysis revealed the quadratic model, including all model terms, to be suitable to represent the tortuosity as a function of δ, m NMC and B/C (R 2 = 0.89). However, the calculated R 2 pred was negative implying that the model is not useful for prediction, mainly because of the inclusion of non-statistically significant terms at a 90% level (∝ > 0.1). The model was then reduced by forward elimination and P-values criterion (∝ < 0.1). Furthermore, the response data was transformed to the inverse function (y transf = 1/y to obtain a normal distribution. The model, after removing non-statistically significant terms (Eq. (31)), shows a strong relationship of tortuosity with the three main factors and the quadratic contribution of the NMC content (P-values < 0.0067), with high correlation and predictive capabilities (R 2 = 0.95 and R 2 pred = 0.86) and adequate precision (signal-to-noise ratio) of 23.1. Although Hamed et al. [138] had already pointed out the significant sensitivity of tortuosity to the slurry formulation and the electrode thickness, the advantage of the DoE approach is that a response surface is obtained whose visual examination can guide further experimental efforts. The response surface shown in Fig. 7a is a 3D representation of Fig. 1 (a)-(b) in the original paper [138]. Fig. 7a also shows that a maximum tortuosity is expected around 92% NMC content, electrodes thicker than 108 µm and higher B/C ratios than 2.5. Further experiments beyond those conditions would be required to determine such a maximum point. Additionally, more experimental points in that region are required to improve the model prediction, since a large deviation between the experimental tortuosity (13)τ = and the modelling value (23)τ = is observed at those conditions.

Rate-capability coefficient
A rate-capability coefficient (γ) was used by Hamed et al. [138] to quantify the reduction of discharge energy by the increase in C-rate. γ was explained mainly as a function of the NMC loading, and a critical value of 92% NMC was reported. The RSM analysis revealed, however, that such critical value strongly depends on the B/C ratio as given by the significantm NMC × B/Cinteraction term (P-value = 0.05) in the model equation (Eq. (32)). At the lowest B/C ratio (1), the critical NMC content is around 89%, whereas at the highest B/C ratio (2.5), the critical NMC is around the reported value of 92%. The real critical NMC value is, nevertheless, outside the studied bounds, as can be seen in the ridge system in Fig. 7b obtained from the plot of Eq. (32). The curvature of the system is mainly provided by the squared term of the NMC content. Although the R 2 is relatively high (0.83) and the adequate precision is 9.9 (values>4 are desirable), the R 2 pred is only 0.48, discouraging the use of the model for prediction aims. Validation experiments would be needed to confirm the validity of the obtained model.

Fractional discharge energy loss
To understand the relative contributions of the short-and long-range transport in the rate-capability limitations, a dimensionless index quantifying the fractional discharge energy loss (ω) was introduced by the authors in the original paper [138]. The ANOVA shows a significant quadratic model for this response (P-value = 0.0022), with thickness and NMC content as the only significant factors (Eq. (33)); i.e., B/C has no effect on ω. The response surface (Fig. 7c) is a 3D representation of the 2D plot in the original paper (Fig. 2 in [138]) and helps to clearly visualise the effect that the factors have on ω. Maximum fractional It is important to point out that although the DoE approach aids in determining the relationship between the different variables, an indepth analysis is always required to truly establish a cause-effect relationship between the factors and the response, as it was done in the work of Hamed et al. [138].

Perspectives and future outlook
Although DoE has been applied to several aspects of the LIB research, as shown in the present review, its application to the different manufacturing processes (mixing, coating, drying, calendering, cutting) has not been widely reported in the academic literature. The understanding of operating parameters on electrode quality and performance is important for process optimisation and ultimate LIB quality and life. The DoE methodology can be applied to the study of manufacturing processes to first determine the main factors influencing the responses, followed by a robust parameter design or an optimisation study. For instance, in the coating process, an initial list of factors including slot die head type, web speed, coating ratio, bump roll gap, web tension, drying temperature and drying air speed can be reduced to the two or three main factors for subsequent analysis.
The obtention of manufacturing data, especially at pilot-plant and industrial scales is, nonetheless, challenging due to the several operating variables involved and the usual time and budget constraints. Fast and efficient designs are therefore particularly important to generate large enough valuable data with the minimum experimental cost. In this regard, designs such as Plackett-Burman can be used to identify the main influencing factors when interactions or higher order terms can be neglected. The more recently proposed Definitive Screening designs [67,140] are also promising in the application of DoE to manufacturing processes since main factors as well as curvature can be detected with the minimum number of experiments.
Hart-to-change factors (HTC) and the need for blocking could be expected in the application of DoE to the manufacturing processes. Electrode drying temperature is an example of an HTC factor whose analysis would involve the application of split-plot designs, not yet reported in the LIBs field. If optimisation is the objective of the DoE and involves an HTC factor, then a RSM with split-plot design would be required. The available literature on this kind of designs and in particular designs that include blocking is, however, rare highlighting the need for further research in this area [141].
To deal with competitive quality criteria normally found in the design of LIBs, for instance, maximum gravimetric capacity but at a minimum electrode thickness, maximum mechanical stability but with minimum insulating polymer binder, the use of desirability functions is recommended since these can be applied to the simultaneous optimisation of multiple responses [142]. Nevertheless, only the work of Meyer et al. [92] has reported the use of desirability functions in the LIBs' field. In particular the authors studied the winding process as a function of the electrodes and separator web tensions.
An alternative to CCD and BBD are the Doehlert designs [143] which present the advantages of mixed level factors and the construction of a new experimental region from previous adjacent points [144]. Doehlert  designs have not been reported yet to the LIB's field, so further studies should be explored to quantify their value and their experimental economy compared with traditional designs and optimal designs. D-optimal designs and RSM are more efficient alternatives to the orthogonal arrays used in finding the optimum pulse charging parameters, largely due to a lower number of experiments and the possibility to detect parameter interactions and second order terms.
Even though DoE approaches have already been used to study ageing in LIBs, the studies so far have only focused on the cell and battery level. Application of the DoE methodology to study degradation at the pack level where several batteries are connected in series and parallel are also important since these are more related to real battery applications [145]. Factors such as chemistry, formulation, cell sizes, temperature, Crate, cooling system, cooling medium, SoC, SoH, storage time and cycling can be studied by optimal designs to identify the key factors, factor interactions, as well as to obtain an empirical model. The experimental data obtained can also be used for parameterisation of existing theoretical models.
Excluding the work of Park et al. [82] on mixture tolerances, there are no DoE papers reported on RPD. RPD studies, specially through RSM, would be a very valuable contribution for the LIBs' manufacturing processes at the different scales (electrode, cell, battery and pack).
The examples presented above highlight future areas of research in which the potential advantage of DoE application has yet to be fully explored and quantified within the field of LIB research. Within the context of LIB manufacturing, the authors propose the study of the individual electrode manufacturing processes to, first, identify the main operating variables and output electrode physical properties at each stage, followed by a series of DoEs to determine optimum settings of the key variables and electrode physical characteristics that will ultimately maximise LIB performance.

Conclusions
The DoE methodology was briefly introduced in this paper and its application to LIBs has been presented by a critical literature review of the available papers on the subject. The review shows how the DoE methodology has already been applied to several aspects of the LIBs field, from active material synthesis and electrode formulation, to the study of influencing electrochemical and physical properties on LIBs' key performance indicators such as capacity and ageing. The review shows how DoE has also helped in the mechanical design of safer batteries by determining suitable operating conditions or to evaluate the addition of safety features. Better charging strategies have also been evaluated by DoE avoiding complex or specialised optimisation algorithms. The best set of parameters have been identified for empirical, semiempirical and theoretical models taking advantage of the maximum information characteristic inherent in DoE.
There are still, however, major opportunities for the application of DoE to LIBs, specially related to manufacturing processes where cost and time are the usual constraints. Efficient designs such as Definitive Screening designs are good candidates for the identification of main factors and curvature in the system. If the number of affordable experiments is fixed, optimal designs could help in determining the best set of experiments.
Ageing studies, which can take several years to complete, can be established using DoE approaches to guarantee that the maximum of relevant information is obtained from the minimum number of experiments. The application of DoE for ageing studies at the pack level are especially needed due to the several factors involved.
Although Taguchi designs are still the most widely used in DoE, even for studies different from RPD, other designs are generally more efficient and should be the preferred option depending on the objective of the study. Optimisation studies, for instance, should be done by RSM. In general, the use of optimal designs is recommended. As exemplified by the case study presented in the review, an I-optimal design can help to understand qualitatively and quantitatively the effect of the input variables on the responses, and to identify optimum regions and settings.
A better understanding of the DoE principles and knowledge of the available designs will result in a wider use of the methodology and ultimately in faster developments in the LIBs' field.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.