Dmitry Gromov

Ingo Bulla

Ethan O Romero-Severson

Silvia Oana

Serea

Numerical optimal control for HIV prevention with dynamic budget allocation $

Keywords: HIV, AIDS, Highly active antiretroviral therapy, Treatment as prevention, Pre-exposure prophylaxis, Optimal control, Resource allocation

This paper is about numerical control of HIV propagation. The contribution of the paper is threefold: first, a novel model of HIV propagation is proposed; second, the methods from numerical optimal control are successfully applied to the developed model to compute optimal control profiles; finally, the computed results are applied to the real problem yielding important and practically relevant results.

Introduction

Since the outbreak of the global HIV/AIDS pandemic in the early 1980s about 35 million people died from AIDS-related illnesses [START_REF] Serea | Global Health Observatory (GHO) data, HIV/AIDS[END_REF]. Although the annual number of new cases of HIV (incidence) has been decreasing globally, among men who have sex with men (MSM) in developed countries incidence is increasing [START_REF] Beyrer | The global response to HIV in men who have sex with men[END_REF][START_REF] Hall | Etimation of HIV incidence in the United States[END_REF]. The current public health challenge in these populations is effective triage of limited prevention resources. Two recent advancements in HIV prevention science are relevant to this issue. First, new forms of highly active antiretroviral therapy (HAART) are so effective that viral load in infected $ The work of O. S. Serea was partially supported by The French National Research Agency, Project ANR-10-BLAN 0112 * Corresponding author Email addresses: dv.gromov@gmail.com (Dmitry Gromov), ingobulla@gmail.com (Ingo Bulla), eoromero@lanl.gov (Ethan O. Romero-Severson), oana-silvia.serea@univ-perp.fr (Oana Silvia Serea) patients become undetectable and that patients on effective treatment are essentially non-infectious [START_REF] Donnell | Heterosexual HIV-1 transmission after initiation of antiretroviral therapy: a prospective cohort analysis[END_REF]; we will refer to this as the Treatment as Prevention (TaP) modality. Second, new treatments targeted at uninfected persons have been shown to reduce the rate of infection to nearly zero when treatment is adhered to [START_REF] Baeten | Antiretroviral Prophylaxis for HIV Prevention in Heterosexual Men and Women[END_REF], which is referred to as Pre-Exposure Prophylaxis (PrEP). We have therapies that are effective at both blocking transmissions from infected persons and protecting uninfected persons from becoming infected; however, there is little consensus on how these technologies should be deployed [START_REF] Cohen | HIV treatment as prevention: Debate and commentary. will early infection compromise treatmentas-prevention strategies[END_REF][START_REF] Jeffrey W Eaton | HIV treatment as prevention: systematic comparison of mathematical models of the potential impact of antiretroviral therapy on HIV incidence in south africa[END_REF]. The main obstacle is that while the individual-level efficacy of these interventions can be ascertained with controlled trials, their effectiveness as public health interventions cannot. This is due to the fact that populations are different from one another both in terms of the fundamental transmission dynamics of HIV, but also in the availability of prevention resources. Identifying the optimal strategy of resource allocation must be based on a model of the underlying medical, biological, and social processes that captures the relevant features of the population.

Theoretical studies of intervention effectiveness are generally based on either state-space models represented by a system of ordinary differential equations (ODE) (often referred to as 'compartmental' models in the epidemiology literature, [START_REF] John | Qualitative theory of compartmental systems[END_REF][START_REF] Wassim | Nonnegative and compartmental dynamical systems[END_REF]) and agent-based models where large-populations of individuals with complex behavioral patters are directly simulated. The trade-off between these approaches exchanges verisimilitude in the agent-based formation for efficient computation in the compartmental formulation. Due to the computational challenges in calculating optimal allocation strategies, compartmental models are generally preferred [START_REF] Kerr | Optima: A model for HIV epidemic analysis, program prioritization, and resource optimization[END_REF][START_REF] Shattock | In the interests of time: improving HIV allocative efficiency modelling via optimal time-varying allocations[END_REF][START_REF] Punyacharoensin | Effect of pre-exposure prophylaxis and combination HIV prevention for men who have sex with men in the UK: a mathematical modelling study[END_REF][START_REF] Griffiths | Analytic solutions to compartmental models of the HIV/AIDS epidemic[END_REF]. In this paper, we present a method for calculating the optimal allocation strategy for a hypothetical health agency in a major US city, which receives funding to set up a long-term HIV prevention program for MSM in that city. Our analysis is based on a compartmental model of HIV transmission that includes natural history of infection, dynamic risk behavior, and partner preference components. For the sake of simplicity we restrict the interventions the health agency can allocate resources to TaP and PrEP. This in turn allows us to address the important question to which extent the investment into PrEP for MSM should be scaled up in the future, which is an open question in HIV prevention policy. While this type of model is in widespread use, our approach to modeling the enrollment of infecteds and high-risk susceptible into ART and PrEP, respectively, and our ability to find dynamic optimal resource allocation patterns over long time horizons is unique.

The described problem of determining the optimal allocation strategy can be formulated as an optimal control problem with a specific set of constraints. We developed an efficient numerical scheme to deal with this problem. The described approach is very general and can be applied to a wide class of optimization problems. To facilitate the reuse of the proposed scheme we provide a detailed description of all steps and indicate possible extensions and ramifications of the method. When solving optimal control problems, there are two classes of methods: indirect and direct ones. Indirect methods employ the firstorder optimality condition to formulate a two-point boundary-value problem for the system of 2N Hamiltonian equation, [START_REF] Boltyanskiy | Mathematical theory of optimal processes[END_REF], which is solved, either analytically or numerically (see, e.g., [START_REF] Pytlak | Numerical methods for optimal control problems with state constraints[END_REF][START_REF] Subchan | Computational optimal control: tools and practice[END_REF]). Analytical solutions to these problems are available only for systems of low order and will not be treated here. The use of numerical schemes is however also restricted. Numerical implementations of indirect methods are often numerically unstable and require a good initial guess which is in most cases not available, [START_REF] Uri M Ascher | Numerical solution of boundary value problems for ordinary differential equations[END_REF]. Furthermore, many non-standard constraints cannot be addressed within the framework of the classical optimal control theory. This is in particular true for the problem considered in the paper. Thus one has to resort to direct methods.

Direct methods attempt to directly minimize the cost function using constrained nonlinear programming. Within the class of direct methods, there are different possible approaches, most notable are shooting methods and simultaneous methods (see [START_REF] Anil | A survey of numerical methods for optimal control[END_REF][START_REF] John | Practical methods for optimal control and estimation using nonlinear programming[END_REF]). Shooting methods parametrize the controls and use numerical simulation to obtain the solution of the system's equations on the whole interval (single shooting) or on a number of subintervals (multiple shooting), [START_REF] Stoer | Introduction to numerical analysis[END_REF]. Nonlinear programming is then used to optimize the control parameters while obeying the constraints. The multiple shooting method has proven to be very efficient for solving many practical problems. However, the necessity of numerical integration of system's equations makes this method very time consuming. Simultaneous methods parametrize both the controls and the system's trajectory and determine the missing values of parameters by solving a (typically) large system of nonlinear algebraic equations. All unknown parameters are determined simultaneously thus giving the name to this class of methods. Along with system's equations, one describes the constraints using the introduced parametrization. This is referred to as the transcription procedure. The optimal control problem is thus formulated as a large scale nonlinear optimization problem, [START_REF] Garg | A unified framework for the numerical solution of optimal control problems using pseudospectral methods[END_REF]. A modification of the latter method was employed in the paper. It was used it to compute a set of optimal allocation strategies for a particular scenario of HIV propagation for different parameters of the treatment.

The paper is organized as follows: in Section 2, a model of HIV propagation is derived and analyzed, Section 3 describes the optimal control problem of resource allocation for HIV treatment and prevention while Section 4 describes the numerical scheme used to solve the optimal control problem. Finally, Section 5 presents and analyzes the numerical results.

Epidemiological model

Derivation of the model

We base our approach on a population balance model. This means that we divide the whole population of MSM in the respective city into a number of groups. All individuals within a given group are assumed to be identical in their evolution. The state variables of the model correspond to the number of people within each group. These are described in Table 1, together with other important notation used throughout the paper. The dynamics of each state variable can be described by the following differential equation:

Ż = V in -V out ,
where V in and V out describe the in-and out-flows. That is the number of people that enter or leave the respective group within the unit time interval. The structure of the flows within the system is shown in Fig. 1, the model parameters are summarized in Table 2. Transitions between states happen due to individuals becoming infected, progressing from acute to chronic stage, and

S H I AH I CH P T H α H φ H u P δ A u T + ūT x y S L I AL I CL T L α L φ L δ A u T + ūT y D δ C δ C ρ H ρ H ρ L µ Figure 1:
The dynamics of the system. In the upper darkgrey box all high-risk states are located, in the lower darkgrey box all low-risk states. All annotations on the arrows are relative transition rates except for u P and u T which are non-rate quantities that govern the rates of the respective transition. There are transitions between high-and low-risk in both directions for susceptibles, infecteds, and people treated with TaP. But people treated with PrEP only can adapt a low risk behavior and become low-risk susceptibles, while there are no transitions from low risk individuals into the group on PrEP treatment. The outflow µ applies to all groups equally.

Parameters used in the model were selected from studies of general MSM populations from the United States. Behavioral parameters including the proportion of the population that is high-risk, the contact rate ratio of high and low-risk individuals, and the rate at which high-risk individuals become low risk and vice versa were taken from an analysis of longitudinal sexual contact rate data of unaffected gay men in 3 large cities in the United States [START_REF] Romero-Severson | Dynamic variation in sexual contact rates in a cohort of HIVnegative gay men[END_REF].

Finally, we selected values of λ L and ūT such that the endemic equilibrium prevalence is about 20% and about 25% of infected individuals are on effective treatment, which is consistent with an average large MSM population in the United States [START_REF] Smith | Prevalence and awareness of HIV infection among men who have sex with men -21 Cities, US 2008[END_REF][START_REF] Rosenberg | Modeling disparities in HIV infection between black and white men who have sex with men in the United States using the HIV care continuum[END_REF]. However, there are several important aspects of the model which deserve particular attention. These are described below.

Mixing. To model infection events (which take place at rate φ H resp. φ L in our model), we assume that transmissions occur at three distinct sites:

• locations frequented exclusively by high-risk individuals 120

• locations frequented exclusively by low-risk individuals

• locations jointly used by both low-and high-risk individuals

The total rate of contact of high-resp. low-risk individuals is denoted by λ H resp. λ L . Hereby, both risk groups make a proportion π of their contacts at the site which is exclusively frequented by their own risk group and the remainder 125 of contacts at the mixing site.

To derive φ H and φ L , we omit the time dependency in the state variables for now. Then, the total number of contacts made at a given time is given by

θ = λ H N H + λ L N L with N H = S H + I AH + I CH + P + T H , N L = S L + I AL + I CL + T L .
The probability that, given a random individual has a contact, this contact is with a high resp. low risk individual is

η H = λ H N H θ , resp. η L = λ L N L θ .
Then, denoting the probability of transmission per contact by β A for acute-stage infecteds and β C for chronic-stage infecteds, the probability of transmission per contact at the mixing site is

130 σ = β A η H I AH N H + η L I AL N L + β C η H I CH N H + η L I CL N L = θ -1 [β A (λ H I AH + λ L I AL ) + β C (λ H I CH + λ L I CL )] .
Consequently, the per-capita instant rate of a high-risk susceptible becoming infected at the common site is

τ H = (1 -π)λ H σ.
Moreover, the per-capita instant rate of a high-risk susceptible becoming infected at the high-risk preferred site is

ψ H = πλ H β A I AH N H + β C I CH N H .
Therefore, the total per-capita transmission rate for high-risk susceptibles is given by

φ H = ψ H + τ H .
Analogously, the per-capita instant rate of a low-risk susceptible becoming infected at the common site is

τ L = (1 -π)λ L σ
and the per-capita instant rate of a low-risk susceptible becoming infected at the low-risk preferred site is

ψ L = πλ L β A I AL N L + β C I CL N L .
The total per-capita transmission rate for low-risk susceptibles is given by

φ L = ψ L + τ L .
Enrollment. For the enrollment on PrEP and TaP, we assume that both is done by randomly sampling individuals at locations where high-risk individuals resp. chronically infecteds are prevalent. In case a sampled individual turns out to be a high-risk susceptible resp. chronically infected, he is urged to enroll in PrEP resp. TaP.

That is, for PrEP we assume that there is a high-risk environment (HRE) where high risk individuals are overrepresented, like bars or sex clubs. The recruitment of the patience then would be carried out during the period in which the location is strongly frequented, i.e., typically in the evening on weekends. Then, denoting the probability that a random high resp. low risk individual is 140 in a HRE at a random moment during the recruitment period by

p H = P (HRE|R = H) p L = P (HRE|R = L),
The probability of a random individual encountered in a HRE at a random moment to be a high-risk susceptible is

P (S H |HRE) = P (HRE|R = H)P (S H ) P (HRE|R = H)P (R = H) + P (HRE|R = L)P (R = L) = p H S H N p H N H N + p L N L N = r b S H r b N H + N L =: ζ P ,
where r b = p H /p L . That is, r b is the odds of a high-risk person to go to a HRE compared to a low-risk person, which we assume to be r b = 0.8/0.2 in our work. Consequently, if u is the relative rate at which individuals are sampled at HREs and then put on PrEP in case they are high risk susceptibles, the absolute rate for transition from S H to P is uN P (S H |HRE) = uζ P N.

Using the same modeling approach for TaP, we obtain the absolute rate for transition from I CH to T H beyond the baseline to be vζ T,H N with

ζ T,H = r b I CH r b N H + N L .
Analogously, the absolute rate for transition from

I CL to T L beyond the baseline is vζ T,L N with ζ T,L = I CL r b N H + N L .
Dynamical system. Based on the considerations presented in the previous paragraphs, we obtain the following system of ODEs describing the dynamics of our

145 system: ṠH = α H -(φ H (X) + ρ H + µ)S H + ρ L S L + xP -u P ζ P (X)N ṠL = α L -(φ L (X) + ρ L + µ)S L + ρ H (S H + P ) İCH = δ A I AH -(ρ H + µ + δ C + v b )I CH + ρ L I CL + yT H -u T ζ T,H (X)N İCL = δ A I AL -(ρ L + µ + δ C + v b )I CL + ρ H I CH + yT L -u T ζ T,L (X)N İAH = φ H S H -(ρ H + µ + δ A )I AH + ρ L I AL İAL = φ L S L -(ρ L + µ + δ A )I AL + ρ H I AH ṪH = -(y + ρ H + µ)T H + v b I CH + ρ L T L + u T ζ T,H (X)N ṪL = -(y + ρ L + µ)T L + v b I CL + ρ H T H + u T ζ T,L (X)N Ṗ = -(x + ρ H + µ)P + u P ζ P (X)N, (1)
where

X = T H T L I CH I CL I AH I AL S H S L P is the vector of state variables; N is the sum of all states, N = 1, X , 1 is the column of ones; φ H (X), φ L (X), ζ T,H (X), ζ T,L (
X), and ζ P (X) are the non-linear (rational) functions of X which take on non-negative values for any X ∈ R n ≥0 ; u P and u T are the control inputs which correspond to the fraction of total population being involved either in PrEP (u P ) or in TaP (u T ), and all the remaining terms are non-negative constants. Below, we analyze an important property of the system (1) that will be used later on.

Nonnegativity

Consider the control system

Ẋ = f (X) + g(X, u), (2)
where f (X) : R n → R n , and g(X) : R n × R m → R n . The system (2) is said to be nonnegative if any solution starting at t 0 from [START_REF] Wassim | Stability and dissipativity theory for nonnegative dynamical systems: a unified analysis framework for biological and physiological systems[END_REF][START_REF] Wassim | Nonnegative and compartmental dynamical systems[END_REF] for more details).

X 0 ∈ R n ≥0 belongs to R n ≥0 for all t ∈ [t 0 , ∞), i.e., ∀X 0 ∈ R n ≥0 , ∀t ∈ [t 0 , ∞), X(t 0 , X 0 , t) ∈ R n ≥0 (see
Definition 2.1. Let f = [f 1 , . . . , f n ] : R n → R n . Then f is essentially nonneg- ative if f i (X) ≥ 0, for all i = 1, . . . , n, and X ∈ R n ≥0 such that X i = 0
, where X i denotes the i-th element of X.

We have the following result:

Theorem 2.1. The system (2) is nonnegative for any nonnegative control u(t) :

[t 0 , ∞) → R m ≥0 if f (X)
and g(X, u) are essentially nonnegative.

Proof. This result can be proved geometrically by examining the direction of the vector field on the bounding hyperplanes

H i = {X ∈ R n ≥0 |X i = 0}.
The essential nonnegativeness condition guarantees that on each H i the system's vector field points towards the positive orthant R n ≥0 .

In our case, we can readily observe that φ H and φ L as well as ζ T,L , ζ T,H , and ζ P are positive for all X ∈ R n ≥0 . Also, we have that the latter three functions turn to zero when, respectively, I CH , I CL or S H is equal to zero. Thus the right hand side of ( 1) is essentially nonnegative and the system's dynamics is non-negative either.

Optimal control problem

The optimal control problem is to minimize the total cost while respecting certain structural and budgetary restrictions. The instantaneous cost is defined as C(t, X(t), U (t)) and the total cost to be minimized is

J C (X) = t f 0 C(t, X(t), U (t))dt, ( 3
)
where t f is the time horizon which is chosen to correspond to the duration of the intervention. In the following we set the initial time to 0, but the problem can be easily modified to account for an arbitrary initial time. Note that a possible problem statement could include a discounting factor e -ρt , ρ > 0 which is typically used to describe our priorities: one may attach greater importance to decreasing the incidence in the near future while paying less attention to what will happen in the farther future. This factor could be used to discount some of the future benefits for shorter term benefits to account to for uncertainty in the validity of the model assumptions in the far future.

Remark 3.1. Note that the optimization problem with the cost functional (3) is well-posed if C(t, X(t), U (t)) is nonnegative for all admissible values of X(t) and U (t). Since the system (1) is nonnegative, an instantaneous cost function defined as a linear combination of the state variables (as it is typical in epidemiological applications) would satisfy this requirement.

The restrictions imposed on the system are twofold:

Restriction on the admissible control policies. The first class of restrictions is due to the structural limitations of the decision unit. Since the intervention is performed by a medical organization the control profile must be sufficiently regular. We assume that the set of admissible controls consists of piecewise constant functions with a fixed interval between two consecutive switches. This restriction can be easily relaxed in different ways: we can assume that the switches occur at some non-regular times, that the control is not piece-wise constant, but rather piece-wise linear (or even piece-wise continuous) and so forth.

Let T = {t i } nint i=0 , 0 = t 0 < t 1 < . . . < t nint = t f be the time instants at which control switches occur along with the initial and final time. We assume that for any 1 ≤ i ≤ n int , the duration of the respective interval is constant: t i -t i-1 = δt. In practice, δt is chosen to be a multiple of 1 year. The set of admissible controls is thus:

U : [0, t f ) → R m ≥0 s.t. U(t) = U i ∈ R m ≥0 , t ∈ [t i-1 , t i ), 1 ≤ i ≤ n int . (4)
The controls are assumed to be right-continuous at t i . In this way, a continuous control is described by a set of its discrete values {U i } i=1,...,nint . The goal of the optimization is to determine these values in order to minimize the cost function (3) while respecting the constraints.

Note that the controls are bounded by zero from below, but there are no upper bounds. The upper bounds are imposed implicitly as will be described below.

Dynamic budget allocation. The second class of constraints is due to the budgetary limitations. In practice, an intervention incurs large expenses which are compensated by the government only to some extent. We assume that at the beginning of each control interval [t i-1 , t i ) the government sets a baseline budget by estimating the expected expenses and allocates the money to be spent for the intervention starting from this baseline. This works as follows:

The total expenditures related to treating people with TaP or PrEP and the enrollment costs for TaP or PrEP are captured by the following cost function:

J B i (X, U i ) = ti ti-1 B i (X(s), U (s))ds,
where B i (X(s), U (s)) are certain positive defined functions. In our case, we formulate the optimal control problem to minimize the incidence of HIV infection. The instantaneous cost is thus defined as the incidence rate of HIV, i.e.,

C(X, U ) = S H φ H (X) + S L φ L (X).
Furthermore, when computing the budgetary restrictions we assumed

B i (X(t), U i ) = K (t) T [T H (t)+T L (t)]+K (t) P P (t)+K (e) T N (t)u i T (t)+K (e)
P N (t)u i P (t)

for i = 1, . . . , n int and t ∈ [t i-1 , t i ). Here, K (t)
T and K We compute this cost for the uncontrolled case to determine the baseline expenses, i.e., the expenses that the government would defray if there is no intervention. The assigned budget is allocated atop the baseline budget. Let Xi (t), t ∈ [t i-1 , t i ) be the uncontrolled (U i = 0) solution of (1) with initial condition Xi (t i ) = X(t i ). The dynamic budget constraint (DBC) is thus formulated as follows:

J B i (X, U i ) -J B i ( Xi , 0) ≤ B, i = 1, . . . , n int , (5)
where U i ∈ R m ≥0 and

J B i ( Xi , 0) = ti ti-1 B i (X(s), 0)ds.
The constraints (5) set an implicit limit to the set of admissible controls U.

In terms of optimal control theory such constraints can be classified as mixed integral inequality path constraints. There is in general no way to handle such constraints analytically, but they can be treated numerically as will be shown below.

Finally we formulate the resulting optimization problem as follows. Determine the discrete values of control

U i ∈ R 2 ≥0 , i = 1, . . . , n int s.t.                J C (X) = t f 0 C(t, X(t), U (t))dt → min J B i (X, U i ) -J B i ( Xi , 0) ≤ B, i = 1, . . . , n int X(t), t ∈ [0, t f ], satisfies (1) with X(0) = X 0 and u(t) = U i , t ∈ [t i-1 , t i ), Xi (t), t ∈ [t i-1 , t i ] satisfies (1) with X i (0) = X(t i ) and u(t) = 0, t ∈ [t i-1 , t i ). (6)

Numerical solution of the optimal control problem

Consider a dynamic system [START_REF] Alistar | Effectiveness and cost effectiveness of oral pre-exposure prophylaxis in a portfolio of prevention programs for injection drug users in mixed HIV epidemics[END_REF]. From now on we will follow the established convention and will assume that the state X(t) is a row vector. Both f (X) and g(X, u) are thus row-valued vector functions.

A collocation method interpolates the state and the control functions at a number of time points (called grid points), and requires the solution to satisfy the respective differential equation at the collocation points which may not necessarily coincide with the grid points. Some grid points can be used to ensure additional conditions on the solution, e.g., continuity.

Lagrange interpolation

The optimal state trajectory is a piece-wise smooth function whose first derivative is discontinuous at the points t i . Therefore it is natural to break it into n int intervals coinciding with [t i-1 , t i ) and interpolate on each interval separately using the basis of Lagrange polynomials, [6]

L i k (t) = ncp l=0, l =k t -τ i l τ i k -τ i l ,
where n cp is the number of collocation points τ i = {τ i k } k=0,...,ncp within the ith interval1 , i = 1, . . . , n int . Since the intervals are of equal length we assume that the number of grid points is the same for each interval.

Figure 2 shows a family of Lagrangian polynomials defined on a non-uniform grid. Notice that for any grid point there is only one polynomial that takes on a non-zero value (which is equal to 1) at this point. Let X(t) be the state trajectory for t ∈ [t i-1 , t i ). We define X i to be the [(n cp + 1) × n] matrix of the values of the state at times τ i k , i.e., X i kj = [X j (τ i k )], j = 1, . . . , n, k = 0, . . . , n cp . The interpolating polynomial for the jth component of the state over the ith interval is thus Xi

j (t) = ncp k=0 L i k (t)X j (τ i k ).
It is well known that for a regular (i.e., equispaced) grid the interpolating polynomial may fluctuate heavily between the interpolation points, especially close to the endpoints of the interval (this is referred to as the Runge phenomenon, see, e.g., [START_REF] Gautschi | Numerical analysis[END_REF]). To overcome this drawback one uses unevenly spaced grid points whose distribution density increases as we approach the endpoints of the interval. There are two standard choices for the grid points: Legendre and Chebyshev points which are zeros of Legendre or Chebyshev polynomials.

These polynomials belong to the class of orthogonal polynomials thus giving the name to the method (orthogonal collocations), [START_REF] Gautschi | Numerical analysis[END_REF][START_REF] Canuto | Spectral methods. Fundamentals in Single Domains[END_REF]. per interval. It is seen that the interpolation becomes nearly exact already with 5 grid points. In contrast to this, when using regular grid points the interpolating function may deviate from the approximated trajectory by several orders of magnitude.

There has been an extensive discussion regarding the merits and drawbacks of either of two choices (see, e.g., [START_REF] Lloyd N Trefethen | Is Gauss quadrature better than Clenshaw-Curtis[END_REF] and references therein). However, it seems that neither of the two is clearly superior to another one. In this study, we have chosen to use the Legendre points.

Collocation at Legendre-Gauss-Radau points

In the previous step the state was described using Lagrange interpolation polynomials which go through (n cp + 1) values of the state X j (t) at t = τ i k .

Determining the values X j (t) corresponds to determining an approximation of the solution of (2). To do so the derivatives of the interpolating polynomials are computed at points τ i k , k = 1, . . . , n cp , i.e., Ẋi j (τ i k ) and the computed derivatives are required to collocate with the right-hand sides of (2) computed for X j (τ i k ). Additionally, the continuity of the trajectory can be ensured by requiring the interpolating polynomials to be attached to each other at the knot points

τ i 0 = t i , i = 1, . . . , (n int -1): Xi-1 (t i ) = Xi (t i ).
To get the required distribution of points we use a particular class of Legendre points, referred to as the Legendre-Gauss-Radau (LGR) points, which are the roots of P ncp+1 (θ) + P ncp (θ), where P k (θ) is the k-th degree normalized Legendre polynomial, P k (1) = 1. The LGR points have the property that one of these points coincides with the left endpoint of the interval, i.e., the respective polynomial has a root at θ = -1 2 . That is, the LGR points are defined on the interval [-1, 1). Remark 4.1. Note that alternatively one can use the Legendre-Gauss (LG) or the Legendre-Gauss-Lobato (LGL) points which are determined as the roots of P ncp (θ) or as the roots of Ṗncp-1 (θ) together with {-1} ∪ {1}. The LG points lie completely within the interval while the LGL points include both endpoints. That is, the LG are defined on (-1, 1) and the LGL points on [-1, 1].

Any t ∈ [t i-1 , t i ) can be associated with θ ∈ [-1, 1) by the affine transfor- mation t = t i-1 + δt 2 (θ + 1). ( 7
)
Consider the (n cp + 1) LGR points θ k , k = 0, . . . , n cp with θ 0 = -1. One can map the LGR points θ k to the respective grid points τ i k using [START_REF] John | Practical methods for optimal control and estimation using nonlinear programming[END_REF], i.e., τ i k = t i-1 + δt 2 (θ k + 1). Following the procedure described above we use Lagrange polynomials to approximate the state trajectory at points τ i k . For the jth component of the state this results in a polynomial of degree at most equal to n cp :

X j (t) ≈ Xi j (t) = ncp k=0 X j (τ i k )L i k (t), t ∈ [t i-1 , t i ). (8)
Differentiating ( 8) and evaluating at the collocation point

τ i k we get Ẋi j (τ i k ) = ncp l=0 X j (τ i l ) Li l (τ i k ) = ncp l=0 X j (τ i l )D i kl , (9)
where D i is an [n cp × (n cp + 1)] differentiation matrix whose (k, l)-th element is the derivative of the Lagrange polynomial L i l at the collocation point τ i k , k = 1, . . . , n cp . Note that we do not collocate at the knot points τ i 0 = t i . The differentiation matrix D extends the idea of finite difference approximations of derivatives which are computed based on the trajectory evaluation at a finite number of points. In our case, the approximation of the derivative is a function of the state at all the n cp + 1 grid points. Remark 4.2. In [START_REF] Garg | A unified framework for the numerical solution of optimal control problems using pseudospectral methods[END_REF], it was shown that the differential matrix D turns out to be singular when the number of collocation points is equal to the number of grid points (here denoted by n gp ). The reason for this is that while the interpolating polynomial is degree n gp -1 its derivative is degree n cp -2 and so, requires only n gp -1 conditions to be uniquely determined. Thus the number of collocation points has to be one less than the number of grid points.

To overcome this difficulty, it is proposed in [START_REF] Garg | A unified framework for the numerical solution of optimal control problems using pseudospectral methods[END_REF] to use the set of LGR or LG points along with a boundary point {-1} or {1}. The interpolating polynomial is computed for the extended set of points while the collocation is carried out only at the LGR, resp. LG points.

The flaw of this approach is that it results in a set of grid points which is not longer produced by an orthogonal polynomial. This leads to a decrease in the accuracy of the resulting polynomial interpolation. The obvious remedy is to compute the interpolating polynomial for the whole set of orthogonal points while collocating at all but one point. The remaining point can be used to enforce the continuity condition as it is done in this paper.

The Lagrange basis polynomials defined for the LGR points can be written in barycentric form [START_REF] Berrut | Barycentric Lagrange interpolation[END_REF] as

L j (t) = P ncp+1 (t) + P ncp (t) (t -τ j ) Ṗncp+1 (τ j ) + Ṗncp (τ j ) ,
whence (see [START_REF] Shen | Spectral methods: algorithms, analysis and applications[END_REF] for the derivation)

Lj (τ k ) =              Ṗncp+1 (τ k ) + Ṗncp (τ k ) (τ k -τ j ) Ṗncp+1 (τ j ) + Ṗncp (τ j ) , j = k, τ k 1 -τ 2 k (n cp + 1)P ncp(τ k ) (1 -τ 2 k ) P ncp+1 (τ k ) + P ncp (τ k ) , j = k = 0, - n cp (n cp + 2) 4 , j = k = 0.
Remark 4.3. Note that the Lagrange basis polynomials are invariant with respect to the shift or the dilatation of the abscissa axis. Thus the approximate differentiation matrix does not change for different intervals provided the number and the type of grid points do not change. However, the use of the affine transformation (7) implies that the system's differential equations should be modified accordingly to take into account the transformed time variable. In practice, this means that the right-hand sides of (2) are to be multiplied with the correction factor δt/2.

Denoting by Ẋi kj the approximate values of the derivatives of X j (t) at collocation points τ i k we write compactly Ẋi = DX i . We also compute the derivatives of the state trajectory by evaluating the right-hand side of the ODE (2). We write

F ( Xi , U i ) kj = f j (X(τ i k )) + g j (X(τ i k ), U i ) and the resulting set of n cp • n constraints is hence DX i - δt 2 F (X i , U i ) = 0. (10)

Defect constraints

While (10) determine values of n • n cp • n int samples of the state trajectory, there are still n • n int free X's which can be used to ensure continuity of the trajectory. These are the values of X at the knot points t = t i , i = 0, . . . , n int -1. The initial value X 0 determines the first n points: X(t 0 ) = X 0 . The remaining values should be determined from the defect constraints which can be written as follows:

X i (t i ) -X i+1 (t i ) = 0 ∀i = 1, . . . , n int -1.
The second term in the left-hand side is X i+1 (t i ) = X(τ i 0 ) and the first term is

X i (t i ) = X i (t i-1 ) + ti ti-1 f (X(s)) + g(X(s), U i )ds. (11)
The integral in [START_REF] Canuto | Spectral methods. Fundamentals in Single Domains[END_REF] can be computed numerically using Gaussian quadrature:

ti ti-1 f j (X(s)) + g j (X(s), U i )ds = δt 2 ncp k=0 w k F j (X(τ i k ), U i ), (12)
where w j are the weights computed for the given distribution of grid points, [START_REF] Canuto | Spectral methods. Fundamentals in Single Domains[END_REF][START_REF] Gautschi | Numerical analysis[END_REF]. Now we have n • (n cp + 1) • n int constraints to determine the same number of discrete values of the state function. Solving these equations is equivalent to 315 getting a numerical solution of the respective differential equations. In other words, we reduced the procedure of integrating the system of ODEs to solving a system of nonlinear algebraic equations. Now we have the solution for given values of the control {U i }, i = 0, . . . , n int -1. Before proceeding to the optimization we have to formulate the constraints imposed on the controls.

Numerical approximation of dynamic budget constraints

To compute the dynamic budget constraints one has to solve the system's ODEs with zero controls for n int intervals with initial conditions determined by the solution of the controlled system. This problem can be formulated within the considered framework in the following way. Let X i 0 be state values at the grid points within the ith interval. We set X i 0 (t i-1 ) = X(t i i-1 ) and the remaining discrete values X i 0 (τ i k ), k = 1, . . . , n cp are determined from

DX i 0 - δt 2 F (X i 0 , 0) = 0. ( 13
)
The resulting set of inequality constraints is thus

δt 2 ncp k=0 w k B i (X i (τ i k ), U i ) -B i (X i 0 (τ i k ), 0) -B lim ≤ 0, i = 1, . . . , n int .

Constrained nonlinear programming problem

The cost function ( 3) is computed numerically using Gaussian quadrature and hence the optimal control problem (6) turns into the following constrained optimization problem:

                                         δt 2 nint i=1 ncp k=0 w k C(X(τ i k ), U i ) → min s.t. DX i - δt 2 F (X i , U i ) = 0 X(τ i 0 ) - δt 2 ncp k=0 w k F j (X(τ i k ), U i ) = 0, i = 1, . . . , n int , DX i 0 - δt 2 F (X i 0 , 0) = 0, i = 1, . . . , n int , X i 0 (t i-1 ) = X(t i i-1 ), i = 1, . . . , n int , δt 2 ncp k=0 w k B i (X i (τ i k ), U i ) -B i (X i 0 (τ i k ), 0) -B lim ≤ 0, i = 1, . . . , n int .
The above nonlinear optimization problem was implemented in Matlab with the use of fmincon function. The optimization was performed using Sequential quadratic programming (SQP) algorithm, see Sec. 5.1 for details.

Results

We calculated the optimal allocation strategy for an outbreak scenario in a large US city with 100,000 at-risk MSM individuals. That is, the introduction of HIV into a subpopulation with a low prevalence (e.g. people aged [START_REF] Garg | A unified framework for the numerical solution of optimal control problems using pseudospectral methods[END_REF][START_REF] Gautschi | Numerical analysis[END_REF][START_REF] Griffiths | Analytic solutions to compartmental models of the HIV/AIDS epidemic[END_REF][START_REF] Wassim | Stability and dissipativity theory for nonnegative dynamical systems: a unified analysis framework for biological and physiological systems[END_REF][START_REF] Wassim | Nonnegative and compartmental dynamical systems[END_REF][START_REF] Hall | Etimation of HIV incidence in the United States[END_REF][START_REF] John | Qualitative theory of compartmental systems[END_REF][START_REF] Kerr | Optima: A model for HIV epidemic analysis, program prioritization, and resource optimization[END_REF][START_REF] Kenneth K Mugwanya | Sexual behaviour of heterosexual men and women receiving antiretroviral pre-exposure prophylaxis for HIV prevention: a longitudinal analysis[END_REF][START_REF] Powers | Impact of early-stage hiv transmission on treatment as prevention[END_REF][START_REF] Punyacharoensin | Effect of pre-exposure prophylaxis and combination HIV prevention for men who have sex with men in the UK: a mathematical modelling study[END_REF] is considered. We assume the infection probability to be β A = 0.015 and β C = 0.001, i.e., acutely infecteds are 15 times more contagious than chronically infecteds [START_REF] Bellan | Reassessment of HIV-1 acute phase infectivity: Accounting for heterogeneity and study design with simulated cohorts[END_REF]. Furthermore, we consider a situation in which risk is static, that is, individuals do not change their risk behavior over time, ρ H = ρ L = 0. Hereby, 90 % of the MSM population exhibits a low risk behavior, the remaining 10 % display high-risk behavior. Hence, we set α L = 250 to be nine times α H = 28 so that the probability that a newly entering individual is high-risk is 10%.

We set µ = 1 360 , leading to the individuals staying 30 years in the system on average if there was no HIV-related removal. Here removal represents either natural death, becoming sexually inactive for medical or social reasons, or settling in a monogamous relationship of two uninfected individuals. There is little data available on the mixing patterns among high and low-risk individuals that we could use to fix π. However, we assumed that π = 0 which is a common assumption when mixing dynamics are unknown. Furthermore, we assume that individuals from the high risk group have ten times more sexual contacts than the ones from the low risk group based on analysis of longitudinal sexual behavioral data [START_REF] Romero-Severson | Dynamic variation in sexual contact rates in a cohort of HIVnegative gay men[END_REF]. Finally, we optimized over values of λ L and ūT such that at equilibrium the prevalence was 20% and the proportion of infected individuals on treatment was 25%, which is consistent with measured values [START_REF] Smith | Prevalence and awareness of HIV infection among men who have sex with men -21 Cities, US 2008[END_REF][START_REF] Rosenberg | Modeling disparities in HIV infection between black and white men who have sex with men in the United States using the HIV care continuum[END_REF], leading to λ H = 40.9, λ L = 4.09, and ūT = 0.00148. We also assume that treatment never fails and is never stopped, i.e., y = 0.

We consider four scenarios varying in the value of x: 0, 1 60 , 1 24 , 1 12 . That is, in one scenario PrEP never fails and is never canceled and in three scenarios this happens on average after 5, 2, and 1 year, respectively. These scenarios are meant to correspond to different distribution policies. More precisely, we assume that PrEP never fails and that a high risk individual, once identified, is prescribed PrEP either indefinitely or is only provided with it for 5, 2, and 1 year, respectively. The parameters of the cost function are given in table 3. The cost of enrollment was estimated by considering the total cost per enrollee of a similar intervention program implemented by the New York City Department of Health [START_REF] Blank | Reaching out to boys at bars: utilizing community partnerships to employ a wellness strategy for syphilis control among men who have sex with men in New York City[END_REF] plus the cost of the labs involved in determining infectious status. The cost of TaP and PrEP is taken from [START_REF] Alistar | Effectiveness and cost effectiveness of oral pre-exposure prophylaxis in a portfolio of prevention programs for injection drug users in mixed HIV epidemics[END_REF]. The trajectories of the controls are shown in Fig. 4, the ones of the number of individuals in the various states in Fig. 567. We can observe that the lower the value of x is, the lower is the total number of newly infecteds and the higher is both the overall values of u P relative to u T and the number of individuals on PrEP. These two dependencies are to be expected as we try to optimize the allocation strategy of the agency in charge of fighting HIV in the considered city and this agency does not cover the long-term cost of prescribing PrEP, but only the enrollment costs. Therefore, a lower x should render the overall intervention more effective and the agency should become more prone to employ PrEP more intensively.

K (t) T K (t) P K (e) T K (
For x = 0, elimination is nearly achieved, whereas for the other values of x this is not the case. That is, only the value x = 0 allows the agency to push the epidemic over the tipping point where the intervention leads to the infection cycle to break down. Moreover, the difference in the number of infecteds over time is much more pronounced between the scenario with x = 0 and the three scenarios with x > 0 than between the three scenarios with x > 0. Most extreme, the difference in the total number of infecteds is only 0.2 % between x = 1 24 and x = 1 12 . This is probably due to buffering effects also observed for other exogenous variables [START_REF] Powers | Impact of early-stage hiv transmission on treatment as prevention[END_REF] leading to a low influenceability of the system for these values of x.

Since only the consequences unfolding over the 50 year period considered are taken into account by the cost function, the choice of how to allocate the resources to enrollment into TaP and PrEP, respectively, becomes myopic to the end of the considered period. That is, for x = 1 24 and x = 1 12 resources are jammed into PrEP towards the end although letting the number of infecteds on treatment drop would ultimately backfire, i.e., leading to number of infecteds higher than necessary after a while. The same effect with the roles of TaP and PrEP switched can be observed for x = 1 60 . Moreover, since the number of infecteds or high-risk susceptibles declines over time when TaP or PrEP is favored, successful enrollment into TaP or PrEP, respectively, becomes more and more expensive, eventually favoring the other treatment. This can be observed for x = 0, x = 1 60 , and x = 1 24 . For x = 1 12 , TaP is favored over PrEP the whole period except for the switch at the very end.

It is therefore advantageous to modify the optimization problem to avoid such unrealistic results. However, as this reformulation is not always feasible, a possible remedy would be to to discard a period long enough at the end of the time interval for which the controls were optimized.

There are two heuristic approaches to accomplish this: (a) Calculating the costs for an only-PrEP and an only-TaP strategy and determining when the last switch between these strategies in terms of which strategy incurs less costs (normally there is only one) takes place. The length of the period discarded is then chosen to be equal to the duration until this switch. (b) Running the optimization procedure on intervals of different length and checking whether the same type of strategy switch taking place at the same distance to the end of the respective calculation period occurs. If yes, this switch is an numerical artifact and needs to be discarded.

Numerical implementation

The described optimization problem was solved both using the multiple shooting (not described in the paper) and the orthogonal collocation methods. In general, the orthogonal collocations approach is about 3-5 times faster than the multiple shooting one3 . However, the latter typically provides a better result albeit the improvement never exceeds a fraction of percent. Furthermore, with orthogonal collocations one may sometimes experience the situation when the program runs out of memory. This is overcome by a slight change of the initial guess. Practice shows that it is in general advantageous to choose an initial guess that provides a low value of the cost function while violating constraints. SQP algorithm recovers from the constraints violation in a couple of steps while keeping the cost function relatively small. If on the contrary, one chooses the initial guess in order to satisfy the constraints, the deviation from the optimal value of the cost function may turn out to be rather large. As the convergence of the algorithm is rather slow, it may take unnecessary many steps to achieve the optimal value. The slow convergence of the algorithm is explained by the particular nonlocal structure of the budget constraints. Since these constraints are computed along the to-be-optimized trajectory, any change of the control variables leads to the variation of the trajectory and thus to the re-computation of the constraints. As an example, a little change of the control during the first interval influences the budget constraints all over the whole optimization interval.

To speed up the computation the gradient of the cost function was supplied to the optimization algorithm. The Hessians were computed numerically using centered finite differences due to the complexity of the respective analytical expressions.

Conclusions

This paper presents a novel model describing an HIV propagation dynamics for a geographically concentrated MSM population, along with two control actions. The two controls represent allocation of funding to the two major drugbased interventions available for HIV, HAART / TaP and PrEP, by an health agency. Addressing this setting is of particular interest as it is heavily debated to which extent PrEP should enter the intervention portfolio of health care systems [START_REF] Punyacharoensin | Effect of pre-exposure prophylaxis and combination HIV prevention for men who have sex with men in the UK: a mathematical modelling study[END_REF][START_REF] Baeten | Antiretroviral Prophylaxis for HIV Prevention in Heterosexual Men and Women[END_REF][START_REF] Kenneth K Mugwanya | Sexual behaviour of heterosexual men and women receiving antiretroviral pre-exposure prophylaxis for HIV prevention: a longitudinal analysis[END_REF][START_REF] David | Preexposure prophylaxis will have a limited impact on HIV-1 drug resistance in sub-Saharan Africa: a comparison of mathematical models[END_REF]. A suitably modified orthogonal collocations method is applied to compute optimal control profiles for a realistic outbreak scenario in a large US city. Hereby, the influence of the effective cost of PrEP on the optimal allocation strategy and the dynamics of the epidemic is studied, by varying a parameter governing said cost. The obtained results show that the allocation pattern heavily depends on the cost of PrEP, rendering PrEP the dominant intervention or not employed at all depending on said cost. Moreover, whether elimination of HIV in the considered population is achievable also is dependent on the effective cost of PrEP.

Currently, we work on improving the developed algorithm in order to introduce it into the epidemiological community. A manuscript aimed at a public health health audience illustrating use of these methods is currently prepared.

100

  dying of AIDS (S • → I A• → I C• → D). Moreover, individuals change their risk behavior ((•) L ↔ (•) H , P → S L ), are put on PrEP or TaP treatment or cancel it (S H ↔ P , I C• ↔ T • ). Finally there is flow into the system due to individuals reaching an age of sexual activity and outflow from the system because of non-HIV related death or individuals becoming sexual inactive or settling in a monogamous, lifelong relationship.

PP

  are the monthly cost for treatment with TaP and PrEP, respectively, per patient. The coefficients K represent the costs for approaching and, if necessary, enrolling one patient into TaP and PrEP, respectively.

Figure 2 :

 2 Figure 2: A family of Lagrange polynomials on the interval [-1, 1]. The grid points are indicated by red circles on the x axis.

Figure 3

 3 illustrates this thesis. It shows a sample trajectory along with its Lagrange interpolation obtained using 3 and 5 non-uniformly distributed points a) b)

Figure 3 :

 3 Figure 3: Interpolation of a sample trajectory using different number of non-uniformly spaced grid points: a) 3 points, b) 5 points. The sample trajectory is shown in blue, its interpolation is in red. Grid points are indicated by circles on the x-axis. Dashed lines separate different intervals.

Figure 5 :Figure 6 :

 56 Figure 5: The trajectories of S H , S L , T H , and T L over 50 years.

Table 1 :

 1 The disaggregation is based on whether an individual is infected and in which stage, 95 what is his risk behavior, and whether he receives treatment and which one. In our model, the time unit is set equal to [month]. Notation used in the paper. The state variables have up to two indices. For variables with one index it indicates the risk status, for variables with two indices infectious stage and risk status are indicated.

	Size variables:
	n	Dimension of the state vector
	n int	Number of solution segments
	n cp	Number of collocation points
	m	Dimension of the control vector
	Points:
	t i , i = 0, . . . , n int Knot points
	τ i k , k = 0, . . . , n cp Grid points
	State vectors:
	X(t)	State evolution of the system,
	X(t) Xi (t)	Polynomial interpolation of X(t) through the grid points τ i k , The solution of the uncontrolled system on the ith interval,
	X i	Matrix of the state values at grid points τ i k
	State variables:
	S • (t)	untreated susceptible individuals
	I • • (t)	untreated infected individuals
	T • (t)	infected individuals on TaP
	P (t)	susceptible individuals on PrEP
	D(t)	individuals deceased due to AIDS
	N (t)	all individuals
	N • (t)	all individuals displaying a given risk behavior
	Indices:
	A/C	Infectious stage (acute/chronic)
	H/L	Risk status (high/low)

Table 2 :

 2

Model parameters. The variables are grouped into the ones directly governing the infection process and the treatment process, respectively, and other variables. The parameters β A , β C , and π are all probabilities, all other parameters are rates, measured in individuals per month.

Table 3 :

 3 Numerical values of the coefficients in the budget functional.

In the following, the upper index i will refer to the number of the respective interval [t i-1 , t i ), i = 1, . . . , n int . We will drop this superscript when the reference to a specific interval is not relevant. Furthermore, the lower index j will refer to the element of the state vector, i.e., j = 1, . . . , n.

One can also define LGR points which include the right endpoint instead the left one.

For the optimization time interval equal to 50 years the computational time for orthogonal collocations was typically within the range of 5e+3-10e+3 seconds while for multiple shooting this time was about 2e+4 -5e+4 seconds. Note that the computational time is strongly influenced by the length of the optimization interval and by the choice of the initial guess.


Bing Xia

Xin Zhao

Raymond De Callafon

Hugues Garnier

Truong Nguyen

Chris Mi

Accurate Lithium-ion battery parameter estimation with continuous-time system identification methods

Keywords: Lithium-ion battery, Equivalent circuit model, Battery management system, Continuous-time system identification, Instrumental variable

The modeling of Lithium-ion batteries usually utilizes discrete-time system identification methods to estimate parameters of discrete models. However, in real applications, there is a fundamental limitation of the discrete-time methods in dealing with sensitivity when the system is stiff and the storage resolutions are limited. To overcome this problem, this paper adopts direct continuous-time system identification methods to estimate the parameters of equivalent circuit models for Lithium-ion batteries. Compared with discrete-time system identification methods, the continuous-time system identification methods provide more accurate estimates to both fast and slow dynamics in battery systems and are less sensitive to disturbances. A case of a 2 nd -order equivalent circuit model is studied which shows that the continuous-time estimates are more robust to high sampling rates, measurement noises and rounding errors. In addition, the estimation by the conventional continuous-time least squares method is further improved in the case of noisy output measurement by introducing the instrumental variable method. Simulation and experiment results validate the analysis and demonstrate the advantages of the continuous-time system identification methods in battery applications.

Introduction

Lithium-ion batteries have been regarded as a promising candidate for energy storage in electric vehicles due to their acceptable energy and power density compared to other types of energy storage devices [START_REF] Lu | A review on the key issues for lithium-ion battery management in electric vehicles[END_REF]. To ensure safe and reliable operation of electric vehicles, the basic conditions of battery packs need to be monitored by the battery management systems (BMS), including voltage, current and temperature. Moreover, some essential states of the batteries also need to be estimated by BMSs, such as state of charge (SoC), state of heath (SoH), and state of power (SoP). Since these states are usually associated with complex electrochemical reactions and cannot be directly measured by physical sensors, researchers need to refer to model-based estimation algorithms.

Among the most commonly applied modeling approaches, the electrochemical model provides high accuracy in estimation since it originates from the first principles of underlying electrochemistry. However, the accurate estimation is achieved at the cost of high complexity in data processing, which requires a significant amount of memory space and computational power to cope with partial differential equations and unknown parameters. This makes this method impractical in real-time applications.

In contrast, the equivalent circuit model (ECM) is widely accepted in the application level because of its simplicity and ease in online implementation. Equivalent circuit modeling utilizes the system identification techniques to relate the input and output behavior of the battery with circuit elements [START_REF] Hu | A comparative study of equivalent circuit models for li-ion batteries[END_REF][START_REF] Sun | A novel dual-scale cell state-of-charge estimation approach for series-connected battery pack used in electric vehicles[END_REF]. An n th -order ECM of Lithium-ion batteries is shown in Fig. 1. The n-RC networks can be interpreted as various time domain characteristics of the physical processes and chemical reactions within the battery cell, the series resistance R 0 is the internal resistance of the battery, and the voltage source represents the electrochemical equilibrium potential at different states, also known as open circuit voltage (OCV) [START_REF] Monem | Lithium-ion batteries: Evaluation study of different charging methodologies based on aging process[END_REF].

All the aforementioned ECM parameters change as SoC changes, but the parameters can be regarded as constant when SoC variation is small, especially in the carefully designed character- ization profiles. The battery characterization determines the best estimates of all the parameters at different SoCs. In general, the OCV is obtained by experiments, and the current and voltage responses in the hybrid pulse power characterization (HPPC) test are used to determine the values of n-RC networks and R 0 at different SoCs [START_REF] Xiong | A data-driven multi-scale extended kalman filtering based parameter and state estimation approach of lithium-ion olymer battery in electric vehicles[END_REF][START_REF] Han | A comparative study of commercial lithium ion battery cycle life in electric vehicle: Capacity loss estimation[END_REF][START_REF] Chen | Loss-minimization-based charging strategy for lithium-ion battery, Industry Applications[END_REF][START_REF] Sun | A systematic state-of-charge estimation framework for multi-cell battery pack in electric vehicles using bias correction technique[END_REF].

v OCV C 1 R 1 R 0 i v bat t + - C n R n C 2 R 2
Different parameter identification methods can be employed to determine the unknowns. The most common approach uses the so-called indirect approach where the parameters of an equivalent discrete-time model are first estimated from which the continuous-time parameters are then derived in a second step. Numerical searching methods utilize global/local searching algorithms to find the best fit parameters, among which the genetic algorithm (GA) and particle swarm optimization are widely applied [START_REF] Seaman | A survey of mathematics-based equivalent-circuit and electrochemical battery models for hybrid and electric vehicle simulation[END_REF][START_REF] Fleischer | On-line adaptive battery impedance parameter and state estimation considering physical principles in reduced order equivalent circuit battery models: Part 1. requirements, critical review of methods and modeling[END_REF]. The disadvantages of such methods are: 1) the convergence to global optimization is not guaranteed given perturbations, and the result is sensitive to initial condition selection; 2) a constant searching time is not assured, which limits their application in real-time. A more robust method is to use curve fitting tools to fit constant current charging or rest periods of experiment data with exponential functions [START_REF] Chen | Accurate electrical battery model capable of predicting runtime and I-V performance[END_REF]. The limitation of this method lies in its ignorance of dynamic inputs in tests, thereby resulting into large voltage errors when dynamic inputs are applied. The standard discrete-time least squares (DT LS) method results in the best estimation of battery parameters with both static and dynamic inputs in terms of minimizing the residual. The DT LS method is faster and requires less memory than numerical searching methods. It has been adapted to recursive least squares (RLS) methods and RLS methods with forgetting factor for online estimation [START_REF] Fridholm | Robust recursive impedance estimation for automotive lithium-ion batteries[END_REF][START_REF] Tang | Li-ion battery parameter estimation for state of charge[END_REF][START_REF] Kim | Model-based condition monitoring for lithium-ion batteries[END_REF][START_REF] Xia | State of charge estimation of lithium-ion batteries based on an improved parameter identification method[END_REF].

As the wide application of the discrete-time model identification in Lithium-ion batteries, researchers have encoutered a limitation that the sampling time needs to be chosen with special care.

If the sampling time is too large, information on the fast dynamics will be lost, whereas too small a sampling time will result numerical problems, that is, the discrete poles accumulate within a small area close to the boundary of the unit circle [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF]. Since the battery system is a typical stiff system with both fast and slow dynamics [START_REF] Remmlinger | State-of-health monitoring of lithiumion batteries in electric vehicles by on-board internal resistance estimation[END_REF], the selection of sampling time is extremely difficult: rapid sampling is required to capture the fast dynamics, however, it will lead to large discrepancies in parameter estimation due to the existence of slow dynamics [START_REF] Garnier | The advantages of directly identifying continuous-time transfer function models in practical applications[END_REF][START_REF] Rao | Continuous-time models and approaches[END_REF]. Given limited storage resolution, when the discrete poles are closer to the stability boundary, the discrete system is more prone to numerical precision of the parameters and the location of the discrete poles [START_REF] Wu | Model parameter estimation approach based on incremental analysis for lithium-ion batteries without using open circuit voltage[END_REF]. It is also concluded by simulation in [START_REF] Yuan | Stability analysis for li-ion battery model parameters and state of charge estimation by measurement uncertainty consideration[END_REF] that the sensitivities to disturbance of the resistance and capacitance values in the RC networks increase quasi-linearly as sampling time decreases, which brings difficulties in accurate system identification.

The direct continuous-time system identification has been studied as opposed to the indirect discrete-time system identification [START_REF]Identification of continuous-time models from sampled data[END_REF][START_REF] Garnier | Direct continuous-time approaches to system identification. Overview and benefits for practical applications[END_REF]. Compared with discrete-time identification methods, the advantages of the continuous-time identification methods are 1) the physical world is continuous and it provides good insight of system properties; 2) the continuous-time identification method avoids discretization that gives rise to undesired high sensitivity issues, and thus can better deal with stiff systems [START_REF] Garnier | Direct continuous-time model identification of high-powered lightemitting diodes from rapidly sampled thermal step response data[END_REF]. More detailed differences between continuous-time and discrete-time identification are provided in [START_REF] Garnier | The advantages of directly identifying continuous-time transfer function models in practical applications[END_REF] [START_REF] Rao | Identification of continuous-time systems[END_REF].

To the authors' best knowledge, the continuous-time approaches have not been analyzed and compared with discrete-time identification methods in battery applications. The novelty of this paper is 1) the continuous-time system identification methods are utilized to identify system parameters of Lithium-ion batteries; 2) the advantages of continuous-time system identification are demonstrated in detail by analytical analysis, simulation and experiment; 3) The instrumental vari-able (IV) method is applied to improve the basic continuous-time least squares estimates. In this paper, the continuous-time and discrete-time n th -order ECM are first introduced. Then the system identifiability and sensitivity of the two methods are compared and discussed. Next, a case is studied to identify the parameters of a 2 nd -order ECM by continuous-time system identification methods. The continuous-time least squares solution is obtained after a state variable filter (SVF) is applied to pre-process the sampled data, and the estimation is further improved with the IV based estimator for coping with the measurement noise. Simulation and experiment results agree with the analysis and demonstrate the advantages of continuous-time parameter estimation methods in real battery applications.

Parameter identification of continuous-time and discrete-time battery models

Continuous-time battery model

The continuous-time transfer function of the n th -order ECM depicted in Fig. 1 can be written as

H(s) = V OCV (s) -V batt (s) I(s) = V(s) I(s) = R 1 R 1 C 1 s + 1 + R 2 R 2 C 2 s + 1 + • • • + R n R n C n s + 1 + R 0 (1)
where V OCV (s) = L{v OCV (t)}, V batt (s) = L{v batt (t)}, I(s) = L{i(t)}, and L{•} is the notation for the Laplace transform.

This function can be rewritten as

(s n + a 1 s n-1 + • • • + a n-1 s + a n )V(s) = (b 0 s n + • • • + b n-1 s + b n )I(s) (2)
where all the coefficients are functions of the unknown circuit parameters. Define

θ = [a 1 • • • a n-1 a n b 0 • • • b n-1 b n ] T = f (R 1 , • • • , R n-1 , R n , C 1 , • • • , C n-1 , C n , R 0 ) (3)
where f : R 2n+1 → R 2n+1 maps the circuit parameters to the equation unknowns.

Equation (2) can be directly written into a regression form as

v (n) (t) = ϕ(t)θ (4)
where v (n) (t) denotes the n th time-derivative of v(t) = v OCV (t)v batt (t), and ϕ(t) is the regression vector given by

ϕ(t) = [-v (n-1) (t) • • • -v(t) i (n) (t) • • • i(t)] (5)
In practice, the measured terminal voltage of batteries is usually contaminated by noises.

Therefore, an equation error e(t) is added to the given regression form. Thus a complete model is given by

v (n) (t) = ϕ(t)θ + e(t) (6)

Discrete-time battery model

To obtain the discrete-time battery model, by using the correct piecewise constant assumption for the current input in HPPC tests, the zero-order hold (ZOH) equivalent of (1) can be derived by [START_REF] Chen | Loss-minimization-based charging strategy for lithium-ion battery, Industry Applications[END_REF].

H(z) = (1 -z -1 )Z H(s) s (7)
where Z{•} is the notation for the z-transform.

Thus the continuous-time transfer function (1) yields to

H(z) = V(z) I(z) = R 1 (1 -e -Ts R 1 C 1 ) z -e -Ts R 1 C 1 + R 2 (1 -e -Ts R 2 C 2 ) z -e -Ts R 2 C 2 + • • • + R n (1 -e -Ts RnCn )
ze -Ts RnCn + R 0 [START_REF] Sun | A systematic state-of-charge estimation framework for multi-cell battery pack in electric vehicles using bias correction technique[END_REF] and it can be rewritten as

(1 + c 1 z -1 + • • • + c n z -n )V(z) = (d 0 + d 1 z -1 + • • • + d n z -n )I(z) (9)
where T s is the sampling interval, and the coefficients in [START_REF] Seaman | A survey of mathematics-based equivalent-circuit and electrochemical battery models for hybrid and electric vehicle simulation[END_REF] are functions of the unknown parameters as well as T s . Define

θ d = [c 1 • • • c n d 0 d 1 • • • d n ] T = g(R 1 , • • • , R n-1 , R n , C 1 • • • , C n-1 , C n , R 0 , T s ) ( 10
)
where g : R 2n+2 → R 2n+1 maps the circuit parameters and T s to the equation unknowns.

Similar to the continuous-time model, a regression form of (9) with a noise term (k) can be written as

v(k) = ϕ d (k)θ d + (k) (11)
where

ϕ d (k) = [-v(k -1) • • • -v(k -n) i(k) i(k -1) • • • i(k -n)] (12)
and v(k) denotes the sampled value of v(t) at time-instant t k .

System identifiability

The structural identifiability describes the uniqueness of parameters given a dynamic model with noise-free and persistent excitations [START_REF] Godfrey | Identifiability of model parameters[END_REF][START_REF] Bellman | On structural identifiability[END_REF]. The battery system is globally identifiable if there is a unique solution θ, θ ∈ Θ ⊂ R 2n+1 , for ( 6) and (11); the system is locally identifiable if there are a finite number of solutions in Θ; the system is unidentifiable if there are an infinite number of solutions in Θ.

The identifiability of Randles ECMs are discussed in [START_REF] Alavi | Structural identifiability of battery equivalent circuit models[END_REF]. The n th -order ECM discussed in this paper has similar structure as that of the Randles models, except for an additional capacitor connecting in series with the RC networks. Similar derivations hold for the ECM of interest in this paper as well.

Given the continuous system transfer function as [START_REF] Lu | A review on the key issues for lithium-ion battery management in electric vehicles[END_REF], if two sets of circuit parameters lead to the same transfer function of an ECM, then

R 1 R 1 C 1 s + 1 + • • • + R n R n C n s + 1 + R 0 = R * 1 R * 1 C * 1 s + 1 + • • • + R * n R * n C * n s + 1 + R * 0 , for all s (13)
A necessary condition for [START_REF] Tang | Li-ion battery parameter estimation for state of charge[END_REF] to hold is

(s + 1 R 1 C 1 ) • • • (s + 1 R n C n ) = (s + 1 R * 1 C * 1 ) • • • (s + 1 R * n C * n ), for all s (14)
Since an n degree polynomial is uniquely characterized by its n distinct roots, the assignments of (-1 14) need to be the permutations of (-1

R * 1 C * 1 ,• • • ,-1 R * n C * n ) in (
R 1 C 1 ,• • • ,-1 R n C n ). Assume the order of the distinct roots is strictly defined, i.e. R 1 C 1 < • • • < R n C n and R * 1 C * 1 < • • • < R n C * n . Since the functions of 1/(s + 1 R i C i ) in (14) are linear independent, R i C i = R * i C * i can be obtained and consequently R i = R * i from (13) for i = 1, 2, • • • , n.
Therefore, all the parameters in the model structure of (1) are uniquely determined. The one-to-one mapping of f in (3) ensures the uniquely determined θ in the reparametrization. Similar derivation holds for discrete time ECM identifiability. Therefore, the conclusion can be drawn that both continuous-time and discrete-time n th -order ECMs are globally identifiable when the order of the RC networks is strictly defined.

Sensitivity of discrete-time model poles

The time-constants of an ECM describe essential system properties of the circuit. Physically, they give the time domain characteristics of voltage variations due to different chemical or physical processes upon charge and discharge. The time-constants τ i in the ECM are calculated by

τ i = R i C i , i = 1, 2, • • • , n (15)
Note that the time-constants are closely related to the eigenvalues, or pole locations of the system as

τ i = 1 |λ i | , i = 1, 2, • • • , n (16)
where λ i is the i th eigenvalue of the system.

In order to characterize the fast dynamics of a system, the sampling frequency should be high enough. A rule of thumb [START_REF] Rao | Continuous-time models and approaches[END_REF] is to select a sampling time such that

f s ≥ 2 × 1 τ min ⇔ |λ max | T s ≤ 1 2 (17)
where λ max is the largest eigenvalue, τ min is its corresponding smallest time-constant in the system, T s is the sampling time and f s is the sampling frequency. In practical applications, the sampling frequency f s is usually set to be higher than the boundary value of 2/τ min to ensure sufficient sampling.

However, a high sampling rate gives rise to undesirable sensitivity issues as illustrated in Fig.

2.

It is well known that the negative real poles in the s-domain map to the interval between the origin and (1, 0) on the real axis of z-domain by [START_REF] Garnier | The advantages of directly identifying continuous-time transfer function models in practical applications[END_REF]. As λ i T s → 0, the poles in the z-domain approaches (1, 0).

z i = e λ i T s (18)
Re(z) Given [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF], the sensitivity of the continuous-time model time-constant with respect to z-domain pole location can be derived as [START_REF] Rao | Continuous-time models and approaches[END_REF].

Im(z) Re(s) Im(s)      s s T z e  1  2  n  Big differ ence Sma ll diff erence 0 s T   1 s T e  2 s T e  n s T e 
S = z τ ∂τ ∂z = e -Ts τ τ ∂τ ∂e -Ts τ = e -Ts τ τ τ 2 T s e -Ts τ = 1 |λ| T s = τ f s (19)
The result shows that the sensitivity is the product of sampling frequency and time-constant.

Without considering the hardware cost, high sampling frequency is desired so that the discretetime system is a good approximation of the continuous-time system. However, equation [START_REF] Rao | Continuous-time models and approaches[END_REF] implies that faster sampling yields to higher sensitivity of the discrete-time model identification to external disturbances and computation errors.

One can also observe from ( 19) that the sensitivity is larger for larger time-constants. In a practical on-board battery monitoring system, the sampling frequency can be 100 Hz, and typical time-constant values for 2 nd -order ECMs are 30 s and 1,000 s [START_REF] Jiang | Evaluation of acceptable charging current of power li-ion batteries based on polarization characteristics, Industrial Electronics[END_REF], thus the sensitivities can be locally as high as 3,000 and 100,000 respectively. It indicates that a small variation in z-domain pole estimation leads to a significant change in time-constant estimation. In the extreme condition, as λT s → 0, the sensitivity goes to infinity.

It is important to state that such a sensitivity issue stems from the discretization process, and it does not exist in the continuous-time model identification cases.

Prescaling in fixed-point storage

Most on-board vehicular microprocessors use fixed-point data storage systems. Such representation limits the information that can be stored during calculation due to precision loss and overflow.

In the continuous-time model identification process, given the stiffness of battery systems, the time-constants vary in one or two orders of magnitude. Since the to-be-estimated coefficients are closely related to the multiplication of different pole locations, these coefficients can be different in several orders of magnitude. If the magnitude of the coefficients is too small, it will lead to big discrepancies in storage, because of the fixed minimum resolution. In order to solve this problem, a parameter prescaling technique is used to maintain the coefficients in the same order of magnitudes, as

θ = Uθ (20)
where the prescaling matrix

U =                          u 1 0 • • • 0 0 u 2 . . . . . . . . . . . . . . . 0 0 • • • 0 u 2n+1                          ( 21
)
and θ is the parameter vector defined in (3).

For instance, if b n in θ has a magnitude in the order of 10 -9 , whereas the processor is 16 bits, i.e. the resolution is 2 -16 = 1.5 × 10 -4 , it is impossible to store b n in such processor for any further calculation. In this case, u 2n+1 in U can be set to 10 9 such that u 2n+1 b n can be stored properly. Then θ can be retrieved by

θ = U -1 θ (22)
In the discrete-time model identification, the dominating poles are close to 1, thus the to-beestimated coefficients have similar order of magnitudes. Therefore the prescaling is not necessary.

Summary of the differences between continuous-time and discrete-time model identification

The major differences of continuous-time and discrete-time model identification are listed in Table 1. In summary, the advantages of continuous-time over discrete-time model identification are a) The identified parameters do not change as T s changes, as shown in ( 3) and [START_REF] Fleischer | On-line adaptive battery impedance parameter and state estimation considering physical principles in reduced order equivalent circuit battery models: Part 1. requirements, critical review of methods and modeling[END_REF], where θ d is related to T s .

b) High sampling rate to avoid information loss does not lead to undesired high sensitivity issues.

c) The continuous-time approach includes inherent filtering, so it is more robust to measurement noises. This is to be discussed in detail in Section 3.

d) Additional advantages have been recently discussed and illustrated in [START_REF] Garnier | The advantages of directly identifying continuous-time transfer function models in practical applications[END_REF] and [START_REF] Garnier | Direct continuous-time approaches to system identification. Overview and benefits for practical applications[END_REF].

The disadvantage of continuous-time model identification are a) The prescaling is required for application in fixed-point storage.

b) More computational power is needed to filter the input/output. This is to be discussed in Section 3. Requires more computation in filtering.

Continuous-time identification of battery models

Conventional least squares based SVF method

The continuous-time model parameter can be obtained from the least squares solution of θ in [START_REF] Han | A comparative study of commercial lithium ion battery cycle life in electric vehicle: Capacity loss estimation[END_REF]. However, the time-derivatives of the input and output are usually not measured. One traditional approach to handle the time-derivatives measurement problem is to use the state variable filter (SVF) method (see e.g. [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF]), which generates filtered (smoothed) versions of the required time-derivatives. A typical SVF is in the form of

L n (s) = α s + α n ( 23
)
where n is the highest system order, and α determines the cut-off frequency of the SVF. The choice of α is recommended to be slightly higher than the estimated bandwidth of the system [START_REF] Garnier | Continuous-time model identification from sampled data: implementation issues and performance evaluation[END_REF].

Applying the SVF to ( 6) yields the linear model as [START_REF] Garnier | Direct continuous-time model identification of high-powered lightemitting diodes from rapidly sampled thermal step response data[END_REF].

v (n) sv f (t) = ϕ sv f (t)θ + e (t) (24)
where

ϕ sv f (t) = [-v (n-1) sv f (t) • • • -v sv f (t) i (n) sv f (t) • • • i sv f (t)], v (n) sv f (t) = L -1 {s n L n (s)V(s)}, i (n) sv f (t) = L -
= N k=1 ϕ T sv f (t k )ϕ sv f (t k ) -1 N k=1 ϕ T sv f (t k )v (n) sv f (t k ) (25)

Improved instrumental variable based SVF method

The least squares solution is unbiased when the equation error is uncorrelated to the regressor, whereas this is usually not true in real applications. In other words, θlssvf converges to the true parameter θ 0 with the assumption of E[ϕ T sv f (t k )e (t k )] = 0, but the least squares regressor ϕ sv f (t) contains the measured output voltage, which correlates with e (t) by [START_REF] Garnier | Direct continuous-time model identification of high-powered lightemitting diodes from rapidly sampled thermal step response data[END_REF]. In order to make the continuous-time parameter identification consistent for correlated equation errors, the instrumental variable (IV) method is introduced. The main idea of the IV method is to find an instrument ζ(t k ) whose components are uncorrelated with e (t k ), i.e. 1 N N k=1 ζ T sv f (t k )e (t k ) = 0 [START_REF] Young | Recursive estimation and time-series analysis: An introduction for the student and practitioner[END_REF]. The most common IV method utilizes an auxiliary model to generate a noise-free estimate of output, as

ζ sv f (t k ) = [-w (n-1) sv f (t k ) • • • -w sv f (t k ) i (n) sv f (t k ) • • • i sv f (t k )] (26)
where w sv f (t k ) is the noise-free estimate of voltage calculated by

w sv f (t k ) = L -1 L n (s)H(s, θlssvf )I(s) (27)
The continuous-time instrumental variable based state variable filter (CT IVSVF) estimator of

θ is given as θivsvf = N k=1 ζ T sv f (t k )ϕ sv f (t k ) -1 N k=1 ζ T sv f (t k )v (n) sv f (t k ) (28)
The flow chart for the implementation of n th -order CT IVSVF battery parameter estimator is summarized in Fig. 3.

It can be noticed that there is no further assumption or limitation of the operating conditions of batteries in continuous-time identification, so long as the battery is in its defined normal operating range. In addition, when dealing with stiff systems, the continuous-time identification method does not require as high storage resolution as discrete-time identification does, because it does not involve numerically dedicate transformation. Therefore, the continuous-time identification methods are expected to have better performance over discrete-time methods, especially when the system is stiff and the storage resolution is limited.

Model parameterization of the 2 nd -order ECM

Based on the analysis in the previous sections, the continuous-time identification methods can be applied to an ECM with any order. In the following sections, the 2 nd order ECM is selected because it is the simplest form to present a stiff system. Other than that, [START_REF] Nejad | A systematic review of lumped-parameter equivalent circuit models for realtime estimation of lithium-ion battery states[END_REF] concludes that the 2 nd order ECM is an optimum choice for implementation of most battery energy and power management strategies.

The continuous-time parameters for the 2 nd -order ECM identification are given as

a 1 = 1 R 1 C 1 + 1 R 2 C 2 a 2 = 1 R 1 C 1 R 2 C 2 b 0 = R 0 (29) b 1 = R 0 1 R 1 C 1 + 1 R 2 C 2 + 1 C 1 + 1 C 2 b 2 = R 0 + R 1 + R 2 R 1 C 1 R 2 C 2
For comparison purposes, the discrete-time parameters are given as

c 1 = -(e -Ts R 1 C 1 + e -Ts R 2 C 2 ) c 2 = e -( Ts R 1 C 1 + Ts R 2 C 2 ) d 0 = R 0 ( 30
)
d 1 = -R 0 e -Ts R 1 C 1 + e -Ts R 2 C 2 + R 1 (1 -e -Ts R 1 C 1 ) + R 2 (1 -e -Ts R 2 C 2 )
d 2 = R 0 e -( Ts R 1 C 1 + Ts R 2 C 2 ) -R 1 e -Ts R 2 C 2 (1 -e -Ts R 1 C 1 ) -R 2 e -Ts R 1 C 1 (1 -e -Ts R 2 C 2 )
Note that the continuous-time parameters are more complex functions of circuit parameters only, while the discrete parameters are functions of both circuit parameters and T s , as indicated in (3) and (10).

v LS algorithm L n (s) s n L n (s) L n (s) s n L n (s) IV algorithm Battery i Auxilliary Model w L n (s) s n-1 L n (s) svf v svf i svf i svf w ˆlssvf  ˆivsvf  ( ) n svf v ( ) n svf i ( ) n svf i ... ... ... ... ... ... ... ( ) n n L s s           ( ) n svf v ( 1) n svf w 

Simulation results

Simulation setup

Simulations are set up in MATLAB to compare the performance of continuous-time and discrete-time identification methods. In the simulation, the circuit of interest is assumed to have R 0 = 1 mΩ, R 1 = 0.3 mΩ, R 2 = 0.6 mΩ, τ 1 = 30 s and τ 2 = 1000 s, and the capacity of the battery is 40 Ah. It needs to be noticed that the simulated system is a typical stiff system, because it includes a fast and a slow time-constants.

The input excitation is shown in Fig. 6c. It consists of five pairs of discharge and charge pulses, each lasting for 20 s. The amplitude of the pulse pairs are 0.5C, 1.0C, 1.5C, 2.0C and 2.5C, respectively, where C is a measure of current at which a battery is discharged to its nominal capacity. The rest periods after each discharge pulses are 40 s, and the rest periods after each charge pulses are 60 s. After the five pair of pulses, a 0.5C 2% SoC discharge is applied, followed by a 1-hour rest.

Various sampling intervals are used in the simulation, which are set to be 1 s, 0.1 s, and 0.01 s, respectively. White Gaussian noises with different signal to noise ratios (SNR) are added to the voltage output in the simulation to emulate the measurement noise of the system, such that

v(k) = v det (k) + n(k) (31)
where v det is the deterministic (or noise-free) voltage signal, n(k) is the measurement noise added to the noise-free voltage signal, and v(k) is the noisy voltage measurement.

The SNR is defined as

SNR dB = 10 log 10 σ 2 v det σ 2 n ( 32
)
where σ 2 v det is the variance of v det and σ 2 n is the variance of the added noises. The standard deviation of noises are chosen to be 0.01 mV, 1 mV and 10 mV respectively. Therefore the corresponding SNR values are 62 dB, 22 dB and 2 dB, respectively.

The SVF cut-off frequency is 0.01 rad/s. Unless otherwise specified, the intermediate data storage for all the simulation is 16 bits, and the prescaling techinique introduced in Sec. 2.5 is applied to all simulation.

Monte Carlo simulations are used to evaluate the performance of the model identification methods, such that 100 different noise realizations are generated for each of the sampling intervals and noise levels. The Bode plots of the identified systems are presented for different methods in various conditions. Except DT LS, CT LSSVF and CT IVSVF methods, it is also interesting to see the performance of DT LS method with the same pre-filtering as in the continuous-time methods, as suggested by [START_REF] Fridholm | Robust recursive impedance estimation for automotive lithium-ion batteries[END_REF]. The Bode plots of estimated models based on discrete-time methods are presented in Fig. 4, and those based on continuous-time methods are presented in Fig. 5.

Discussion of simulation results

Fig. 4a and Fig. 4b give the Bode plots for the CT models identified by the indirect DT LSbased method. The identified models are clearly biased at the low frequency ends. As the noise level or sampling rate increases, the identified systems are more biased.

The Bode plots for the systems identified by the DT LS methods with both input and output filtering are given in Fig. 4c and Fig. 4d. The filter applied is the same as L 2 (s) with same cut-off frequency as SVF. It is hard to conclude solely from the two plots whether the pre-filtering improves the system identification, but at least the pre-filtering does not provide the desired system identification. Later simulation with 64 bit storage resolution in Fig. 4e and Fig. 4f demonstrates the effectiveness of the pre-filtering in high storage resolutions. The difference in the identification performance indicates that the pre-filtering is largely influenced by the rounding errors in discretetime identification.

Fig. 5a and Fig. 5b present the identification results of CT LSSVF method in various conditions. The performance of system identification is remarkably improved, which is due to the application of data pre-filtering and avoidance of discretization. At last, Fig. 5c and Fig. 5d show that CT IVSVF identification has accurate identifications in all situations, where the Bode plots of the identified systems overlap that of the true system with a high fidelity except when the SNR is equal to 2dB when a bias can be observed in the low-frequency part of the response.

In summary, most of the indirect discrete-time identification methods lead to large discrepancies in system identification, especially in the presence of noises and small sampling intervals; the simple CT LSSVF estimates are clearly less biased due to SVF pre-filtering and low sensitivity to

Experimental results

Battery tests

In order to identify the parameters of a Lithium-ion polymer cell, battery characterization experiments are conducted at room temperature (22 • C∼25 • C). The specification of the battery under tests is given in Table 2. The procedure of the characterization tests is shown in Fig. 6a.

The sampling rate in all the experiments is 1 Hz.

The capacity test consists of three charge/discharge cycles. In each cycle, the battery is charged with a constant current of 0.5C until the terminal voltage reaches the maximum charge voltage.

Then the voltage is kept at the maximum charge voltage until the charge current is below 1/20C.

After that, the battery is discharged at 0.5C until the minimum discharge voltage is reached. A 1 hr rest is applied after each charge/discharge operation. The charge/discharge capacity are calculated as the average value of the three cycles.

The charge/discharge OCV-SoC curves are obtained at 10% SoC step with 0.5C charge/discharge rate and 5 hr rest, as shown in Fig. 6b. It can be observed that the charge OCV curve is higher than the discharge curve due to the hysteresis effect. The maximum difference between charge/discharge curves is 22.6 mV. The average value of charge/discharge OCV curves is used in the model calculation for simplicity. The HPPC test profile is the same as mentioned in Section 4. The HPPC test is repeated between 10% and 90% SoC with 2% SoC step. The current input and voltage output of the HPPC pulses at SoC=50% are shown in Fig. 6d.

Finally, Urban Dynamometer Driving Schedule (UDDS) cycles with 10 minutes rest periods are applied consecutively to the battery cell to validate the identified model. The initial SoC for UDDS test is 90%, and the test is terminated after the first cycle when SoC drops below 20%. The current profile of UDDS is scaled such that the maximum discharge current is 100 A.

Discussion of experimental results

In order to demonstrate the advantages of continuous-time system identification, the performances of CT LSSVF and CT IVSVF estimators are compared with DT LS method and DT LS with pre-filtering. The cut-off frequency for both of the filters are 0.025 rad/s and the storage resolution is 16 bits.

The voltage outputs of the simulated models and experiments in UDDS are plotted in Fig. 7a, and the zoomed-in range between 5,900 s and 7,400 s are shown in Fig. 7b. The results show that i) The indirect discrete-time methods have the larger overall error.

ii) The DT LS method with pre-filtering has similar noise levels to the CT LSSVF method in some of the cycles, however, the systems identified are clearly biased in the 5 th , 6 th and 9 th cycle. The defect in robustness is mainly caused by its high sensitivity to unmodeled dynamics and colored noises in the experiment.

iii) The two direct continuous-time identification methods have better performance over other methods. Furthermore, CT IVSVF estimator gives better estimation because the newly introduced instrumental variable is less correlated to equation error. Based on the findings in simulation results, the DT LS method tends to estimate the poles of the system to be larger, then the corresponding time-constants are smaller. This can be observed in It needs to be noticed that DT LS presents acceptable results in experiments because the UDDS cycle mainly consists of high frequency contents, where it gives reasonable identification in the simulation as well. However, it has large bias in larger time-constant estimates, which will lead to larger discrepancies when the input is constant. This can be visualized by its larger discrepancies at beginning of rest periods between UDDS cycles.

The mean absolute error (MAE) and root mean squares error (RMSE) of the aforementioned estimators within the whole experimental period are listed as in Table 3. It validates that the continuous-time identification methods have better overall performance and CT IVSVF method gives voltage estimations with the smallest error.

Conclusion

This paper has presented the advantages of continuous-time system identification methods in battery ECM parameter estimation. The combination of fast dynamics induced by charge transfer and slow dynamics from diffusion causes the dynamics of a battery to exhibit stiff system behavior.

Continuous-time identification reduces the sensitivity of the model parametrization in case of a stiff dynamic system, leading to improvements on battery parameter estimation in simulation and experiment results, especially when battery dynamic parameters are stored in fixed precision.

The general modeling of both continuous-time and discrete-time n th -order ECM is first discussed. Then it is shown that the n th -order ECM is identifiable if the order of the time-constants is strictly defined. The comparison between continuous-time and discrete-time methods indicates that discrete-time identification methods are less robust due to undesired sensitivity issues in transformation of discrete domain parameters. In continuous-time parameter identification, SVF is used This research bridges the gap between battery parameter estimation algorithm development and real applications, where the system is stiff and the storage resolutions are limited. It is important to note that the offline continuous-time parameter estimation can be easily adapted online.

With a more accurate battery model, the estimation of other battery states can be improved accordingly, such as SoC, SoH and SoP. The future work includes the implementation of the improved recursive continuous-time system identification methods to battery online state estimation.

Figure 1 :

 1 Figure 1: n th order ECM for Lithium-ion batteries.

Figure 2 :

 2 Figure 2: Discrete pole converges to (1, 0) in z-plane as λT s → 0.

  1 {s n L n (s)I(s)}, and L -1 {•} denotes the inverse Laplace transform. Based on N sampled measurements of both input and output signals at time instants t k , k = 1, ..., N, the continuous-time least squares based state variable filter (CT LSSVF) estimator of θ is expressed as θlssvf

Figure 3 :

 3 Figure 3: Implementation flow of CT IVSVF battery parameter estimator.

  DT LS method for different SNR (T s =0.1 s, 16 bits).

  DT LS method for different T s (SNR=22dB, 16 bits).

  DT LS method with pre-filtering for different SNR (T s =0.1 s, 16 bits).

  DT LS method with pre-filtering for different T s (SNR=22dB, 16 bits).

  DT LS method with pre-filtering for different SNR (T s =0.1 s, 64 bits).

  DT LS method with pre-filtering for different T s (SNR=22dB, 64 bits).

Figure 4 :

 4 Figure 4: Monte Carlo simulation with different SNR and T s . Bode plots of CT models from indirect DT methods.

  CT LSSVF method for different SNR (T s =0.1 s, 16 bits).

  CT LSSVF method for different T s (SNR=22dB, 16 bits).

  CT IVSVF method for different SNR (T s =0.1 s, 16 bits).

  CT IVSVF method for different T s (SNR=22dB, 16 bits).

Figure 5 :

 5 Figure 5: Monte Carlo simulation with different SNR and T s . Bode plots of CT models from direct CT methods.

  OCV of the battery under test. HPPC pulses at SoC=50%.

Figure 6 :

 6 Figure 6: Battery characterization and validation experiments.

  Detailed voltage values and errors for 5,900 s and 7,400 s.

Figure 7 :

 7 Figure 7: Discrepancies of different identification methods.

Fig. 8

 8 Fig. 8 as well, where the time-constant estimates of DT LS are the smallest among all the methods.

Figure 8 :

 8 Figure 8: Time-constant estimates from experimental data with different methods.

Table 1 :

 1 Major differences between discrete-time (DT) and continuous-time (CT) model identification.

	DT

c) θ d changes as T s changes. θ remains unchanged. d) DT models are derived from CT, and go back to CT. Remains as CT.

e) Small T s or large τ make the system identification No undesired sensitivity issues.

Table 2 :

 2 Specification of battery under test.

	Material	Lithium-ion polymer
	Charge/discharge capacity	40.83/40.61 Ah
	Nominal voltage	3.7 V
	Maximum charge voltage	4.2 V
	Minimum discharge voltage 2.7 V
	rounding errors, and the CT IVSVF estimates are both accurate and robust in all situation.

Table 3 :

 3 MAE and RMSE of model voltage estimation in whole UDDS validation. time-derivative terms. After that, the IV method is applied to further increase the estimation accuracy by reducing the correlation between equation error and regressor.Simulation results show that the continuous-time identification methods demonstrate higher accuracy in stiff system estimations given different sampling intervals and noise levels, because of the inherent SVF pre-filtering and the avoidance of pole discretization. Characterization experiments are conducted and different identification methods are used to identify the battery parameters, including DT LS, DT LS with pre-filtering, CT LSSVF and CT IVSVF. The results indicate that the continuous-time identification methods results in the smallest mean absolute error and root mean square error. Among all, the best output voltage prediction is obtained by CT IVSVF method.

		DT LS DT LS filter CT LSSVF CT IVSVF
	MAE [mV]	8.59	8.71	6.34	4.27
	RMSE [mV] 9.35	10.36	6.96	4.95

Acknowledgement

The authors would like to acknowledge the funding support from Nanjing Golden Dragon Bus Co., Ltd., and the US Department of Energy under the Graduate Automotive Technology Education Center program.