University of Southern California Dissertations and Theses

Distributed parameter model based system identification and filtering in the estimation of blood alcohol concentration from transdermal alcohol biosensor data

(USC Thesis Other)

Distributed parameter model based system identification and filtering in the estimation of blood alcohol concentration from transdermal alcohol biosensor data

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Distributed Parameter Model Based System Identication and Filtering in the Estimation of Blood Alcohol Concentration from Transdermal Alcohol Biosensor Data by Jian Li A thesis submitted to the Department of Mathematics in partial fulllment of the requirements for the degree of Master of Science in the subject of Statistics University of Southern California Los Angeles, California August 2017 © 2017 Jian Li. All rights reserved. Acknowledgements This thesis is a summary of the work I have done in the past two years, which could not be possibly completed without the guidance, help, inspiration and love from many people around me, directly or indirectly. I would like to express my deepest gratitude to all of you. My advisor, Professor Gary Rosen, opened a new horizon to me, a pure engineer who never dove into mathematics before. Not only his academic teaching established a solid foundation of my skills in statistics, but also his professional advise guided me throughout the way of exploring my research. It was my honor to work with and co-author a few of papers with Professor Susan E. Luczak. Her expertise made up my lack of knowledge and experience in alcohol research. I also would like to thank Professor Chunming Wang for his great help and support for being a member of my committee. I want to thank all my current and former labmates - especially Chitresh Bhushan, Minqi Chong, Hossein Shahabi, Divya Varadarajan - for creating a wonderful working environment and your sel ess discussion. I must also acknowledge the USC Annenberg Graduate Fellowship Program and the National Institute of Health for providing generous nancial support for my graduate studies. I would like to thank all sta in the Department of Mathematics and the Ming Hsieh Department of Electrical Engineering for their administrative support. Finally and most importantly, I am grateful for the unlimited love and endless sacrice from my wife and parents, which made me strong and helped me chase my dream. ii Table of Contents Acknowledgements ii List Of Tables iv List Of Figures v Abstract vi Chapter 1: Introduction 1 Chapter 2: Estimating Blood Alcohol Concentration from Transdermal Alcohol Biosensor Measurements 4 2.1 A Quasi-Blind Deconvolution Problem . . . . . . . . . . . . . . . . . . . . . . . . . . 5 2.2 A Distributed Parameter Model for the Transdermal Transport of Ethanol . . . . . . 6 Chapter 3: Determining the Convolution Kernel or Impulse Response Function 9 3.1 Distributed Parameter Modeling and the Impulse Response Function . . . . . . . . . 9 3.1.1 Abstract Parabolic Input/Output Systems, their Associated Impulse Re- sponse Functions, and their Identication . . . . . . . . . . . . . . . . . . . . 10 3.1.2 Finite Dimensional Approximation, the Adjoint Method, and Dierentiating the Matrix Exponential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12 3.1.3 Application to the System Discussed in Section 2.2 . . . . . . . . . . . . . . . 14 3.2 Frequency Domain Techniques and the Impulse Response Function . . . . . . . . . . 15 3.3 ARMA Modeling and the Impulse Response Function . . . . . . . . . . . . . . . . . 16 Chapter 4: Result 19 4.1 Experimental Data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19 4.2 A Deconvolution Technique . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 4.3 Numerical Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 Chapter 5: Discussion and Concluding Remarks 24 Reference 26 iii List Of Tables 4.1 Five error statistics of the prediction using (a) ^ h PDF (t), (b)h FFT (t) and (c) ^ h ARMA (t) 21 iv List Of Figures 4.1 The 11 drinking episodes including the sensor measured transdermal alcohol concen- tration (TAC) and the contemporaneously measured breath alcohol concentration (BrAC) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.2 (a) The estimated impulse response ^ h PDF (t) (purple), ^ h FFT (t) (green) and ^ h ARMA (t) (light blue). (b) The estimated TAC from convolution with training BrAC using ^ h PDF (t) (purple), ^ h FFT (t) (green) and ^ h ARMA (t) (light blue). The true TAC is shown in blue. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22 4.3 Prediction of BrAC (red) via deconvolution from TAC (blue) using ^ h PDF (t) (purple), ^ h FFT (t) (green) and ^ h ARMA (t) (light blue) on the ten drinking episodes not used in training. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 4.4 Box plots of the ve statistics: the relative errors of peak, delay, slope 1, slope 2 and AUC of the prediction result using (a) ^ h PDF (t), (b) h FFT (t) and (c) ^ h ARMA (t) . . . 23 v Abstract Methods for the estimation of blood or breath alcohol concentration (BAC/BrAC) from biosen- sor measured transdermal alcohol concentration (TAC) are developed, evaluated and compared. Specically, a scheme based on a distributed parameter model with unbounded input and out- put for ethanol transport in the skin is compared to more conventional ltering/deconvolution techniques, one based on frequency domain methods, and the other on a time series approach using an autoregressive moving average (ARMA) input/output model. Our basis for comparison are ve statistics of interest to alcohol researchers and clinicians: peak BAC/BrAC, time of peak BAC/BrAC, the ascending and descending slopes of the BAC/BrAC curve and area underneath the BAC/BrAC curve. It can be shown that the ARMA-based method yields the best estimation of the peak while the distributed parameter model produces the best estimation of the time of the peak. The Fourier-based method has the least variance out of the three and is computationally very ecient. vi Chapter 1 Introduction Distributed parameter systems (i.e. partial and functional dierential equations) have been used for modeling, simulation, control, estimation, identication, and optimization in a variety of applications arising in engineering, biology, chemistry, physics, space science, and economics. The appeal of using these innite dimensional systems to describe the evolution of the state is based on a number of factors. These include: (1) their ability to yield high delity models with relatively low dimensional parameterization, (2) the fact that they typically result when modeling is based on rst principles and physics-based empirical observations, and (3) that there is a rich, functional analytic theory upon which the analyses of well-posedness (i.e. existence, uniqueness, continuous dependence), stability, long-term behavior, and nite dimensional and numerical approximation and convergence, can be based. In addition, although distributed parameter models result in an innite dimensional state equation, by exercising some level of mathematical precision they can readily be interfaced with nite dimensional control inputs and output observations, which are what are typically encountered in practice. There is a large body of research on the application of distributed parameter systems in control science reported on in the literature, but the overwhelming majority of it is highly theoretical. There are relatively few papers describing how these models perform when they actually are used to design a controller, identify a system, or simulate a process. More to the point, there is even less in the literature that compares this approach to using more traditional nite dimensional techniques. 1 2 We present and discuss our experience using a distributed parameter system as the basis for a data analysis system for a biosensor that measures transdermal alcohol concentration (TAC), or more specically, the number of ethanol molecules evaporating from the surface of the skin. The role of the data analysis system is to use the TAC signal from the biosensor to estimate the concentration of ethanol in the blood. Blood alcohol concentration (BAC) or breath alcohol concentration (BrAC) is currently the basis for all clinical research on the eects of alcohol on human physiology and behavior. BrAC serves as a proxy for BAC which is more dicult and invasive to measure. In addition, both BAC and BrAC are used as quantitative determinants of whether an individual is legally intoxicated or driving under the in uence of alcohol. A passive, non-invasive biosensor that provides researchers and clinicians with an accurate and continuous estimate of BAC in real time from the eld would be considered a major breakthrough and is currently the focus of several academic and commercial research eorts. Our approach here is to formulate the problem as a single input single output (SISO) system in which the alcohol in the blood is the source of the input and the biosensor measured transder- mal alcohol is the output signal. The objective of estimating the BAC then takes the form of a system identication/input determination, or blind deconvolution problem. We model the forward transport of ethanol through the skin using a one-dimensional diusion equation with input and output on the boundary. Linear semigroup theory and the innite dimensional variation of con- stants formula are used to obtain a parametric representation for the unknown forward convolution lter. We then show how the approximation theory for linear semigroups can be used to develop nite dimensional convergent numerical approximation schemes for estimating the convolution ker- nel and deconvolving BAC from TAC. Finally, we compare our approach to two more conventional approaches to the problem, a Fourier/frequency domain based approach and a scheme based on an autoregressive moving average (ARMA) input/output response surface. Our comparison is based on actual clinical breath and TAC data collected in the eld by one of the co-authors (S.E.L.) using a transdermal biosensor and a breath analyzer. An outline of the remainder of the paper is as follows. In Chapter 2 we discuss the problem of estimating BAC from TAC, its formulation as a blind deconvolution problem, and modeling it using Chapter 1. Introduction 3 a distributed parameter system. In Chapter 3 we develop our scheme for estimating the convolution kernel and the deconvolution of BAC from TAC. We also discuss the frequency domain and ARMA based approaches. In Chapter 4 we present our numerical results and a nal fth chapter contains discussion and concluding remarks. Chapter 2 Estimating Blood Alcohol Concentration from Transdermal Alcohol Biosensor Measurements Continuous monitoring of quantitative alcohol levels in naturalistic settings would be extremely valuable for advancing a wide range of alcohol research domains. These include how alcohol con- centrations relate to drinking motives, physical responses to alcohol, behaviors, decision-making, negative consequences, and environmental/situational factors over the course of a drinking episode and across drinking episodes. Such knowledge would further the understanding of alcohol-related health outcomes, disease progression, treatment ecacy, and recovery. In addition to the benet this would provide to researchers and clinicians, being able to quantitatively monitor real-time alcohol levels is now recognized as broadly applicable to lay individuals monitoring their personal health and for public safety. There is currently no known procedure for the accurate and unobtrusive eld collection of quantitatively interpretable alcohol consumption data. Collection of blood and urine samples in the eld is not practical and breath analyzers and drink diaries have either heavy subject burden or large time gaps between assessment points, hinder naturalistic drinking behavior, and may not be accurate even when individuals are trying to be compliant. Breath analyzer readings are frequently too high due to mouth alcohol or too low due to not taking deep lung breaths, while drink diaries are inaccurate due to not knowing the alcohol content of a drink or the amount consumed. 4 Chapter 2. Estimating BAC from TAC 5 2.1 A Quasi-Blind Deconvolution Problem Biosensor devices that measure TAC, the amount of alcohol diusing through the skin, are prov- ing to be a promising technology. Unfortunately, however, to date, they are primarily being used as abstinence monitors because TAC data does not consistently correlate with BrAC or BAC across individuals, devices, and environmental conditions. To wit, there is currently no well-established method for producing reliable estimates of BrAC/BAC (eBrAC/eBAC) from TAC data. A breath analyzer relies on a relatively simple model from basic chemistry (i.e., Henry's Law) for the ex- change of gases between circulating pulmonary blood and alveolar air. This model has been found to be reasonably robust across subjects and the environment. The transport and ltering of alcohol by the skin on the other hand, is physiologically more complex and based on a number of factors that dier across individuals (e.g., skin layer thickness, tortuosity) and drinking episodes within individuals (e.g., temperature, skin hydration, vasodilation). This means that, regardless of how reliable and accurate transdermal alcohol devices become at measuring TAC, the raw TAC data will never trivially map directly onto BrAC/BAC across individuals and drinking episodes. In order to address this problem, we have been developing a data analysis system to accompany the sensor hardware that converts TAC to BrAC or BAC. Our approach is based on a rst principles physics-based model (in the form of a distributed parameter system or partial dierential equation (PDE)) for the transport of ethanol molecules from the blood through the epidermal layer of the skin and its measurement by the sensor. The result is a forward input/output model that can be formulated as either a continuous or discrete time convolution in which the convolution kernel or lter, while unknown, is dened in terms of the parameters that appear in the physics based model. Obtaining an estimate for the BAC or BrAC can then be formulated as a blind or quasi blind (since the lter is dened parametrically) deconvolution problem. The current protocol for using the system involves simultaneously collecting calibration BrAC (by a breath analyzer) and TAC (by using the sensor) in either the lab or clinic for a single drinking episode before sending the patient or research subject out in the eld with just the TAC sensor. Estimated BrAC from the TAC for all subsequent drinking episodes recorded in the eld is then obtained by using the 6 2.2. A Distributed Parameter Model for the Transdermal Transport of Ethanol calibration data to t the parameters that determine the lter, and then using the t lter to deconvolve estimated BrAC from the TAC. 2.2 A Distributed Parameter Model for the Transdermal Transport of Ethanol Let '(t;) denote the concentration in moles/cm 2 of ethanol in the interstitial uid in the epidermal layer of the skin at depth cm and time t in seconds. Let L denote the skin thickness in cm. We model the transport of ethanol through this medium as a diusion process @' @t (t;) =D @ 2 ' @ 2 (t;); 0<<L; t> 0 (2.1) where D > 0 denotes the diusivity in units of cm 2 /sec. For boundary conditions, at the skin surface ( = 0), we model evaporation of the alcohol vapor using a Robin boundary condition; the ux (from right to left) is proportional of ethanol at the boundary of the epidermal layer, D @' @ (t; 0) ='(t; 0); t> 0 (2.2) where > 0 denotes the constant of proportionality in units of cm/sec. At the interface of the dermal layer (which has a blood supply) and the epidermal layer (which does not contain blood vessels) ( = L) we impose a Neumann boundary condition. That is that the ux (from right to left) is proportional to the BAC, D @' @ (t;L) =u(t); t> 0 (2.3) where the parameter > 0 characterizes the exchange of ethanol molecules between the blood and the interstitial uid in units of moles / (cm sec BrAC/BAC units), and u denotes the concentration of ethanol in the blood as given in BAC (or BrAC) units. The boundary condition Chapter 2. Estimating BAC from TAC 7 (2.3) serves as our input condition. We assume that there is no alcohol in the epidermal layer at time t = 0 which yields the initial conditions '(0;) = 0; 0<<L (2.4) We model the processing by the TAC sensor of the ethanol evaporating from the surface of the skin via a linear relation which serves as an output condition y(t) = '(t; 0); t> 0 (2.5) where denotes the constant of proportionality in units of TAC units cm 2 /mole. As it stands, the model given by the equations (2.1) - (2.5) is determined by ve parameters: D;L;; and . However, not all ve of the parameters are independent nor are they uniquely identiable from the input/output data. Without loss of generality and by converting to what are essentially dimensionless quantities, the number of unknown parameters to be t can be reduced to two, which we denote by the vector q = [q 1 q 2 ] T [1]. For simplicity, leaving the names of the variables unchanged (although their units certainly are), in this way, our input/output model becomes @' @t (t;) =q 1 @ 2 ' @ 2 (t;); 0<< 1; t> 0 (2.6) q 1 @' @ (t; 0)'(t; 0) = 0; t> 0 (2.7) q 1 @' @ (t; 1) =q 2 u(t); t> 0 (2.8) '(0;) = 0; 0<< 1 (2.9) y(t) ='(t; 0); t> 0 (2.10) The input/output system (2.6) - (2.10) is an example of a distributed parameter system. We note that there is a feature of the model that makes it somewhat non-standard. That is that the input and the output are on the boundary. In the next chapter, when the model is reformulated 8 2.2. Distributed Parameter Modeling abstractly as an abstract evolution equation in an innite dimensional Hilbert space, this results in the input and output operators being unbounded in the standard space where an initial boundary value problem such as (2.6) - (2.10) is typically set. This requires some care when we re-cast the problem using techniques from functional analysis, specically, linear semigroup theory. Also, reformulating the problem in discrete time is some help in this regard as well. Chapter 3 Determining the Convolution Kernel or Impulse Response Function In this chapter, we describe three techniques for obtaining and estimating the impulse response function for the input/output model for the dynamical system described in the previous chapter. We rst consider an approach based on the rst principles physics-based distributed parameter model discussed in Section 2.2. We then look at two more conventional approaches: a method based on frequency domain techniques and one based on an ARMA model. 3.1 Distributed Parameter Modeling and the Impulse Response Function In this section, we provide a brief summary of how tools from functional analysis and in particu- lar, the theory of linear semigroups of operators can be used to transform the distributed parameter model presented in the previous chapter into a linear input/output model in the form of a convolu- tion. We also describe how the resulting impulse response function, or convolution kernel or lter, which is dened in terms of operators on innite dimensional spaces, can be approximated via ma- trix representations for linear operators dened on a sequence of approximating nite dimensional subspaces. A more detailed presentation of the results discussed here including relevant theorems and proofs can be found in [2] and [1]. 9 10 3.1. Distributed Parameter Modeling 3.1.1 Abstract Parabolic Input/Output Systems, their Associated Impulse Response Functions, and their Identication LetV andH be Hilbert spaces that satisfy the dense and continuous embeddingsV ,!H ,!V , whereV denotes the space of continuous linear functionals that is dual to V . Leth;i denote the H inner product. For q2fQ;dg, a compact metric space, let a(q;;) :VV !R be a bilinear form satisfying the following three conditions: 1. (Boundedness) ja(q; 1 ; 2 )jjj 1 jjjj 2 jj; 1 ; 2 2V 2. (Coercivity) a(q; ; ) +j j 2 jj jj 2 ; 2V 3. (Continuous Dependence) ja(q 1 ; 1 ; 2 )a(q 2 ; 1 ; 2 )jd(q 1 ;q 2 )jj 1 jjjj 2 jj; 1 ; 2 2V; q 1 ;q 2 2Q Forq2fQ;dg letb(q;) :V !R andc(q;) :V !R be linear and continuous and consequently, therefore, it follows thatb(q;) =b(q) for someb(q)2V andc(q;) =c(q) for somec(q)2V . We then consider the input/output system in weak form given by h _ '; i +a(q;'; ) =b(q; )u; 2V; y =c(q;') (3.1) or equivalently h _ '; i +a(q;'; ) =hb(q); i V ;V u; 2V; y =hc(q);'i V ;V (3.2) whereh;i V ;V denotes the natural extension of theH inner product to the duality pairing between V and V . If we set W (0;T ) =f : 2 L 2 (0;T;V ); 2 L 2 (0;T;V )g and u2 L 2 (0;T ) it can Chapter 3. Convolution Kernel 11 be shown [3] that the system (3.1) or (3.2) admits a unique solution '2 W (0;T ) that depends continuously on u2L 2 (0;T ). It follows that W (0;T )C(0;T;H) and that y2L 2 (0;T ). For q2Q, the q-dependent bilinear form on VV , a(q;;) :VV !R, denes a bounded linear operator A(q)2L(V;V ) byhA(q) 1 ; 2 i =a(q; 1 ; 2 ), for 1 ; 2 2 V . Then, if we letH denote either of the spaces H or V , we can consider the linear operator A(q) to be the unbounded linear operator, A(q) : D q H ! H where D q = V in the caseH = V , and D q =f 2 V : A(q) 2Hg in the caseH = H. It can then be shown [4], [5], [6] that A(q) is a closed, densely dened unbounded linear operator onH and it is the innitesimal generator of an analytic semigroup of bounded linear operators,fe A(q)t :t 0g onH. For q 2 Q, dene the bounded linear operators B(q) : R ! V and C(q) : V ! R by hB(q)v; i V ;V = b(q; )v =hb(q); i V ;V v, and C(q) = c(q; ) =hc(q); i V ;V , respectively, for 2V and v2R. The input/output system can now be written formally in the standard form in H as _ x(t) =A(q)x(t) +B(q)u(t) y(t) =C(q)x(t); t> 0 (3.3) where the state x(t) = '(t;) 2 H. Assuming that the system is initially at rest (i.e. that '(0;) = 0), and using the fact thatfe A(q)t : t 0g is an analytic semigroup onH and therefore that e A(q)t 2 D q V , for 2 V , we obtain from the abstract variation of constants formula that y(t) = Z t 0 C(q)e A(q)(ts) B(q)u(s)ds = Z t 0 h(q;ts)u(s)ds (3.4) where h(q;t) =C(q)e A(q)t B(q); t> 0. Let a sampling time > 0 be given and consider zero order hold inputs of the form u(t) =u i , t2 [i; (i + 1)); i = 0; 1; 2;::: (typically u i =u(i); i = 0; 1; 2;::: , where u is a given continuous time input). Setx i =x(i) and lety i =y(i); i = 0; 1; 2;::: . It then follows in the usual way that x i+1 = ^ A(q)x i + ^ B(q)u i ; y i =Cx i ; i = 0; 1; 2;::: (3.5) 12 3.1. Distributed Parameter Modeling with ' 0 = 0 2 V , where ^ A(q) = e A(q) 2 L(V;V ), and ^ B(q) = R 0 e A(q)s B(q)ds 2 L(R;V ). Boundedness of the operators ^ A(q) and ^ B(q) follows from the fact thatfe A(q)t :t 0g is an analytic semigroup on V;H and V [4], [5], [3], [6]. If A(q) : D q V ! V is invertible with bounded inverse, it follows that ^ B(q) = R 0 e A(q)s B(q)ds =A(q) 1 e A(q)s B(q) 0 = ( ^ A(q)I)A(q) 1 B(q) and therefore that y i = i1 X j=0 C ^ A(q) ij1 ^ B(q)u j = i1 X j=0 C ^ A(q) ij1 ( ^ A(q)I)A(q) 1 B(q)u j ; i = 0; 1; 2;::: (3.6) or that y i (q) = i1 X j=0 ^ h ij (q)u j ; i = 0; 1; 2;::: (3.7) where ^ h i (q) =C ^ A(q) i1 ( ^ A(q)I)A(q) 1 B(q), i = 0; 1; 2;::: Given training data,f~ u j ; ~ y j g; j = 0; 1; 2;:::;N:, we formulate the identication problem as a nonlinear least squares t to data. That is we seek q2fQ;dg which minimizes the cost functional J(q) = N X i=0 jy i (q) ~ y i j 2 (3.8) where u i = ~ u i ; i = 0; 1; 2;:::;N. 3.1.2 Finite Dimensional Approximation, the Adjoint Method, and Dierentiating the Matrix Exponential Computing the value of the cost functionalJ(q) and its gradient for a given value ofq2fQ;dg as would be required to solve the optimization problem posed above requires the nite dimensional approximation of the innite dimensional operators that appear in the denition of ^ h i (q). For n = 1; 2;::: , let H n = spanf n j g n j=0 V , and let P n :H!H n denote the orthogonal projection of H onto H n with respect to the H inner product. We assume that thef n j g n j=0 are such that lim n!1 P n = in H for 2 H and in V for 2 V . For n = 1; 2;::: , and q2 Q, dene A n (q)2L(H n ;H n ) to be the nite dimensional linear operator whose matrix representation is given bybA n (q)c i;j =[h n i ; n j i] 1 [a(q; n i ; n j )], fori;j = 0; 1; 2;:::;n. Under certain conditions, Chapter 3. Convolution Kernel 13 in particular when A(q) 1 exists, and in condition 2 above is non-positive, it is not dicult to show that A n (q) = (P n a A(q) 1 ) 1 , where P n a is the orthogonal projection of V onto H n with respect to the inner producth;i a = a(q;;) on V . We set ^ A n (q) = e A n (q) and obtain the nite dimensional approximating discrete time input/output system given by x n i+1 = ^ A n (q)x n i + ( ^ A n (q)I)A n (q) 1 ^ P n B(q)u i y n (q) =C' n i ; i = 1; 2;::: (3.9) where x n 0 = 02H n , and ^ P n denotes the standard bounded extension of P n to V . It follows that y n i (q) = i1 X j=0 C( ^ A n (q)) ij1 ( ^ A n (q)I) ^ P n B(q)u j = i1 X j=0 ^ h n ij (q)u j ; i = 0; 1; 2;::: (3.10) where ^ h n i (q) = C( ^ A n (q)) i1 ( ^ A n (q)I) ^ P n B(q); i = 1; 2;::: . With the help of the Trotter-Kato semigroup approximation theorem, it can be shown that ^ h n i (q)! ^ h i (q) as n!1, uniformly in q for q2fQ;dg and uniformly in i for i in bounded subsets ofZ. We then seek q2fQ;dg which minimizes the cost functional J n (q) = N X 0 jy n i (q) ~ y i j 2 (3.11) where u i = ~ u i ; i = 0; 1; 2;:::;N. The cost functional J n is minimized iteratively. It is clear that for a given value of q, the value ofJ n (q) can now easily be computed. The gradient ofJ n (q) is computed using the adjoint method [7]. For i = 0; 1; 2;:::;N, set v n i = [2(Cx n i ~ y i ); 0;:::; 0] T 2R n+1 and dene the adjoint system z n i1 = [ ^ A n (q)] T z n i +v n i1 ; i =N;N 1;:::; 1; z n N =v n N (3.12) 14 3.1. Distributed Parameter Modeling The gradient of J n (q) can then be computed as ! 5J n (q) = N X i=1 [z n i ] T ( @ ^ A n (q) @q x n i1 @ ^ A n (q) @q ( ^ A n (q)I) ^ P n B(q)~ u i1 ^ P n B(q)~ u i1 ( ^ A n (q)I) @ @q ( ^ A n (q)I) (3.13) The tensor @ ^ A n (q) @q can be computed at the same time. ^ A n (q) is computed using the sensitivity equations. For t 0 and q 2 Q, set n (q;t) = e A n (q)t . Then n (q;) is the unique principal fundamental matrix solution to the initial value problem _ n (q;) =A n (q) n (q;); n (q; 0) =I (3.14) Setting n (q;t) = @ n (q;t)=@q, dierentiating with respect to q, interchanging the order of dierentiation, and using the product rule, we obtain _ n (q;) =A n (q) n (q;) + @A n (q) @q n (q;); n (q; 0) = 0 (3.15) Combining these two initial value problems and solving we obtain 2 6 4 @ ^ A n (q) @q ^ A n (q) 3 7 5 = 2 6 4 n (q;) n (q;) 3 7 5 =e 2 6 6 6 4 A n (q) (@A n (q)=@q) 0 A n (q) 3 7 7 7 5 2 6 4 0 I 3 7 5 (3.16) 3.1.3 Application to the System Discussed in Section 2.2 Let Q be a closed and bounded subset of R 2 endowed with the Euclidean metric, let H = L 2 (0; 1) together with the standard inner producth 1 ; 2 i = R 1 0 1 (x) 2 (x)dx, and norm denoted byjj, and let V be the Sobolev space V = H 1 (0; 1) together with its standard inner product hh 1 ; 2 ii = R 1 0 1 (x) 2 (x)dx + R 1 0 0 1 (x) 0 2 (x)dx and normjjjj. Then we have the usual dense Chapter 3. Convolution Kernel 15 and continuous embeddingsV ,!H ,!V , whereV denotes the space of distributions dual toV . The forms and functions a(q;;) :VV !R;b(q;) :V !R and c() :V !R are given by a(q; 1 ; 2 ) = 1 (0) 2 (0) +q 1 Z 1 0 0 1 (x) 0 2 (x)dx; 1 ; 2 2V (3.17) b(q; ) = q 2 (1), and c( ) = (0), for 2 V . It follows that b(q) = q 2 ( 1) 2 V and c(q) =2V , where denotes the Dirac delta distribution, or unit impulse at zero. With regard to nite dimensional approximation, for n = 1; 2;::: letf n j g n j=0 denote the set of standard linear B-spline on the interval [0; 1] dened with respect to the usual uniform mesh, fj=ng n j=0 , and set H n = spanf n j g n j=0 V (note the n j are the usual \pup tent" or \chapeau" functions of height one and support of width 2=n,b(j 1)=n; (j + 1)=nc T b0; 1c). IfP n :H!H n denotes the orthogonal projection of H =L 2 (0; 1) onto H n , it is well known (see for example [8]) that lim n!1 P n = in H for 2H and in V for 2V . 3.2 Frequency Domain Techniques and the Impulse Response Function In this section, we solve the kernel estimation problem using a signal processing approach. Non- parametrically we can rewrite equation (3.7) in a traditional convolution equation in the continuous domain as follows: y(t) =h(t)u(t) = Z t 0 h(t)u()d (3.18) The convolution Theorem [9] states that if h(t) has a Fourier transform H(f) and u(t) has a Fourier transform U(f), then h(t)u(t) has a Fourier transform H(f)U(f). Therefore, the convolution equation in Section 3.1.1 can be written in a frequency domain representation as Y (f) =Ffy(t)g =Ffh(t)u(t)g =Ffh(t)gFfu(t)g =H(f)U(f) (3.19) whereFfg is the forward Fourier transform. Then we can estimate the spectrum ofh(t), i.e. H(f), directly as follows ^ H(f) = Y (f) U(f) = Ffy(t)g Ffu(t)g (3.20) 16 3.3. ARMA Modeling and the Impulse Response Function The estimated spectrum ^ H(f) must be low-pass ltered with a certain cut-o frequency before taking the inverse Fourier transform back to the time domain for the following two reasons. First, the estimated spectrum ^ H(f) is very noisy due to the limited number of samples in the training dataset. Second, according to the Nyquist-Shannon sampling theorem [10], if the sampling fre- quency F s is twice higher than or equal to the maximum frequency of the actual signal, then that band limited signal can be exactly recovered by low-pass ltering with cut-o frequency at F s =2. In other words, the signal of interest in the frequency domain contains no more information beyond F s =2 provided that the Nyquist-Shannon condition is satised. That means the non-zero entries beyond F s =2 in ^ H(f) comes purely from the noise component of the signal and therefore we can safely reduce the noise via a low-pass ltering with cut-o frequency at F s =2 without losing any information about the signal. Mathematically, this procedure can be expressed as ^ H 0 (f) = ^ H(f)W LPF (f) (3.21) where W LPF (f) is the low-pass ltering window W LPF (f) = 8 > > < > > : 1; f <F s =2 0; fF s =2 (3.22) Finally, the estimated convolution kernel, denoted as h FFT (t), can be obtained from ^ h FFT (t) =F 1 f ^ H 0 (t)(f)g (3.23) whereF 1 fg is the standard inverse Fourier transform. 3.3 ARMA Modeling and the Impulse Response Function We can also model the underlying system as an ARMA process as follows y(t) = (t) (t) u(t) +e(t) (3.24) Chapter 3. Convolution Kernel 17 where(t) and(t) are the auto-regressive and moving average coecients respectively, e(t) is the estimation error or as it is sometimes called, the residue or noise, u(t) is the driving process or the input to the system, and y(t) is the observed process or the output of the system. In our case, BrAC is the u(t) and TAC is the y(t). The system can be identied and the impulse response can be obtained provided (t) and (t) are accurately estimated. There are various ways to estimate the parameters (t) and (t). For example, the AR coecients(t) can be estimated by solving the well-known Yule-Walker equation using Levinsons recursion. The MA coecients (t) can be estimated via Durbins algorithm. If the data are jointly Gaussian, then maximum likelihood estimation technique can also be applied to estimate the coecients. In this work, to better incorporate delays between BrAC and TAC due to the nature of the data, we adopt the method from [11] which is brie y summarized here. • We set maximal possible AR and MA order to be 60 and 120 heuristically based on the data. In fact, it turns out that our estimation results are not very sensitive to the choices of these two maximal orders. • We decrease the AR order one at a time forming a set of lower order models and use Minimum Description Length (MDL) criterion to evaluate and select the best one, the one with the lowest MDL score. This implicitly estimates the AR order while the MA order is xed at the maximal order of choice at this step. • We set the AR order to be the best one obtained in the previous step and remove MA coecients one at a time from the one with lowest S/N ratio dened in equation (18) in [11], forming a set of lower order models and then use MDL to select the best one. This step estimates the MA order, thereby establishing the overall ARMA order. • After nding the optimal orders, we t the ARMA model based on a least-squares approach to obtain (t) and (t). Once(t) and(t) are obtained from the steps above, the impulse responseh(t) can be estimated directly from y(t) by setting u(t) = (t), where (t) is the standard Dirac delta function. The 18 3.3. ARMA Modeling estimated ^ h(t), termed ^ h ARMA (t) for this method, is also low-pass ltered with the same cuto frequency F s =2 as described in Section 3.2. Chapter 4 Result In this chapter, we applied the three methods to the actual data described below in Section 4.1 and then evaluated their performance. We rst performed a test on the training session to check how well the three models could be t to the data. This was done by convolving the estimated ^ h(t) with the training BrAC and comparing the result against the training TAC. Second, to further test the performance of the estimation, we used ^ h(t) to deconvolve the testing TACs in another ten drinking episodes of data that were not used in training and then compared the deconvolution results with the contemporaneously collected BrACs. The deconvolution was done via solving an optimization problem described in Section 4.2. Then statistics were drawn based on the comparison results from the ten testing episodes of data and the numerical results are shown in Section 4.3. 4.1 Experimental Data One of the co-authors (S.E.L.) wore a WrisTAS™ 7 alcohol biosensor for 18 days while also collecting breath measurements. The WrisTAS™ 7 measures the local ethanol vapor concentration over the skin surface at 5-minute intervals. It looks like a digital watch. The participant consumed her rst drink in the laboratory with BrAC being measured and recorded every 15 minutes from the start of the drinking session until BrAC returned to 0.000. She then wore the TAC device in the eld and consumed alcohol ad libitum for the following 17 days. For each drinking episode, the subject would take BrAC readings every 30 minutes until the BrAC returned to 0.000. Fig. 4.1 shows the entire 18 day TAC signal along with the contemporaneous BrAC measurements. The 19 20 4.2. A Deconvolution Technique TAC measurements provided by the sensor are in units of milligrams per deciliter (mg/dl), while the BrAC measurements are in units of percent alcohol. 4.2 A Deconvolution Technique A deconvolution is a process that estimates the input u(t) given the output y(t) and the esti- mated kernel ^ h(t). The deconvolution result ^ u(t) is the answer to the question: what is the prole of BrAC given the measurements in the form of TAC? In order to fairly compare the three meth- ods, we need a deconvolution technique that is independent of all three kernel estimation methods, i.e. the deconvolution process is not biased toward any of the three kernel estimation methods. Therefore, we performed the deconvolution via solving the following optimization problem: ^ u(t) = argmin u(t) jjA(t)u(t)y(t)jj 2 l 2 + 1 jju(t)jj l 2 + 2 jj5u(t)jj l 2 ; s:t: u(t) 0 (4.1) where A(t) is the convolution operator (Toeplitz convolution matrix in discrete domain) formed from ^ h(t),5 is the gradient operator and 1 and 2 are two regularization parameters to impose smoothness of the data and its rst order derivative. We also restrict the result to be non-negative due to the nature of the data and the estimation problem. This optimization problem is a con- strained convex problem thus can be solved eciently via a non-negative least squares technique (Matlab function LSQNONNEG). 4.3 Numerical Results The three estimated impulse responses ^ h PDF (t), ^ h FFT (t) and ^ h ARMA (t) are shown in Fig. 4.2a. The convolution results on the single training episode of data are shown in Fig. 4.2b. Fig. 4.3 shows results of the prediction of BrAC in the other ten episodes of data which are totally unavailable during the training session. Fig. 4.4 shows the statistics of the relative errors of the ve measures as box-plots across ten episodes for (a) ^ h PDF (t), (b) ^ h FFT (t) and (c) ^ h ARMA (t), respectively. The associated table of statistics is shown in Table 4.1. Chapter 4. Result 21 Table 4.1: Five error statistics of the prediction using (a) ^ h PDF (t), (b) h FFT (t) and (c) ^ h ARMA (t) Peak Delay Slope 1 Slope 2 AUC Mean 0.325 0.387 0.579 0.567 0.333 ^ h PDF (t) Median 0.312 0.326 0.454 0.496 0.302 S.D. 0.268 0.256 0.417 0.468 0.169 Mean 0.291 0.418 0.456 0.519 0.302 ^ h FFT (t) Median 0.310 0.359 0.502 0.431 0.304 S.D. 0.182 0.315 0.294 0.519 0.129 Mean 0.319 0.421 0.528 0.552 0.286 ^ h ARMA (t) Median 0.264 0.294 0.501 0.433 0.320 S.D. 0.251 0.338 0.354 0.444 0.123 22 4.3. Numerical Results Figure 4.1: The 11 drinking episodes including the sensor measured transdermal alcohol concen- tration (TAC) and the contemporaneously measured breath alcohol concentration (BrAC) (a) (b) Figure 4.2: (a) The estimated impulse response ^ h PDF (t) (purple), ^ h FFT (t) (green) and ^ h ARMA (t) (light blue). (b) The estimated TAC from convolution with training BrAC using ^ h PDF (t) (purple), ^ h FFT (t) (green) and ^ h ARMA (t) (light blue). The true TAC is shown in blue. Chapter 4. Result 23 Figure 4.3: Prediction of BrAC (red) via deconvolution from TAC (blue) using ^ h PDF (t) (purple), ^ h FFT (t) (green) and ^ h ARMA (t) (light blue) on the ten drinking episodes not used in training. (a) (b) (c) Figure 4.4: Box plots of the ve statistics: the relative errors of peak, delay, slope 1, slope 2 and AUC of the prediction result using (a) ^ h PDF (t), (b) h FFT (t) and (c) ^ h ARMA (t) Chapter 5 Discussion and Concluding Remarks Qualitatively, we can see that ^ h FFT (t) has more curvature and is less smooth than ^ h ARMA (t), and than ^ h PDF (t), successively. In addition, the corresponding recovered TAC is closer to the true one using ^ h FFT (t) than it is using either ^ h ARMA (t) or ^ h PDF (t). Quantitatively we evaluate the TAC recovery performance by measuring the error, in l 2 norm sense, between the three estimated TAC curves and the true TAC. It turns out that ^ h FFT (t) yields an error of 32.388, ^ h ARMA (t) yields and error of 44.898 and ^ h PDF (t) yields an error of 65.373. This is the same trend observed from visualizations. To quantitatively test the performance of the predictions and compare the three kernels, we take the following ve measures: (1) the height of the peak (peak), (2) the time peak (delay), (3) the ascending slope (slope 1), (4) the descending slope (slope 2) and (5) the area under the curve (AUC). Then we dene a relative error of any of measurement as jM \ BAC M BAC j M BAC 100% (5.1) where M \ BAC and M BAC are the measurements obtained from the estimated BrAC and the true BrAC. Overall, the three methods are very comparable. The ARMA-based method yields the most accurate estimation of the peak while the PDE-based method produces the best estimation of the delay. The Fourier-based method has the least variance out of the three, especially for the estimations of the two slopes, although it has slightly larger bias. It is also worth noting that the Fourier-based method is very computationally ecient but theoretically has innitely many 24 Chapter 5. Discussion 25 parameters to estimate. On the other hand, the PDE-based method requires the estimation of only two parameters but requires more time to t the data. The ARMA-based method lies in between those two in terms of computational eciency and the number of parameters. Reference [1] Z. Dai, I. Rosen, C. Wang, N. Barnett, and S. E. Luczak, \Using drinking data and phar- macokinetic modeling to calibrate transport model and blind deconvolution based data anal- ysis software for transdermal alcohol biosensors," Mathematical Biosciences and Engineering, vol. 13, pp. 911{934, jul 2016. [2] I. G. Rosen, S. E. Luczak, and J. Weiss, \Blind deconvolution for distributed parameter systems with unbounded input and output and determining blood alcohol concentration from transdermal biosensor data," Applied Mathematics and Computation, vol. 231, pp. 357{376, 2014. [3] J. L. Lions, \Optimal Control of Systems Governed by Partial Dierential Equations," 1971. [4] H. T. Banks and K. Ito, \A Unied Framework for Approximation in Inverse Problems for Distributed Parameter Systems.," tech. rep., DTIC Document, 1988. [5] H. T. Banks and K. Kunisch, Estimation Techniques for Distributed Parameter Systems. Springer Science & Business Media, 1989. [6] H. Tanabe, Equations of evolution, vol. 6. Pitman Publishing, 1979. [7] A. F. J. Levi and I. G. Rosen, \A Novel Formulation Of The Adjoint Method In The Optimal Design Of Quantum Electronic Devices," SIAM Journal on Control and Optimization, vol. 48, no. 5, pp. 3191{3223, 2010. [8] M. H. Schultz, \Spline Analysis," p. 192, 1973. [9] R. N. Bracewell, The Fourier transform and its applications, vol. 31999. McGraw-Hill New York, 1986. [10] C. E. Shannon, \Communication in the Presence of Noise," Proceedings of the IRE, vol. 37, no. 1, pp. 10{21, 1949. [11] M. H. Perrott and R. J. Cohen, \An ecient approach to ARMA modeling of biological systems with multiple inputs and delays,"IEEETransactionsonBiomedicalEngineering, vol. 43, no. 1, pp. 1{14, 1996. 26

Abstract (if available)

Abstract Methods for the estimation of blood or breath alcohol concentration (BAC/BrAC) from biosensor measured transdermal alcohol concentration (TAC) are developed, evaluated and compared. Specifically, a scheme based on a distributed parameter model with unbounded input and output for ethanol transport in the skin is compared to more conventional filtering/deconvolution techniques, one based on frequency domain methods, and the other on a time series approach using an autoregressive moving average (ARMA) input/output model. Our basis for comparison are five statistics of interest to alcohol researchers and clinicians: peak BAC/BrAC, time of peak BAC/BrAC, the ascending and descending slopes of the BAC/BrAC curve and area underneath the BAC/BrAC curve. It can be shown that the ARMA-based method yields the best estimation of the peak while the distributed parameter model produces the best estimation of the time of the peak. The Fourier-based method has the least variance out of the three and is computationally very efficient.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Finite dimensional approximation and convergence in the estimation of the distribution of, and input to, random abstract parabolic systems with application to the deconvolution of blood/breath al...

PDF

Obtaining breath alcohol concentration from transdermal alcohol concentration using Bayesian approaches

PDF

An abstract hyperbolic population model for the transdermal transport of ethanol in humans: estimating the distribution of random parameters and the deconvolution of breath alcohol concentration

PDF

Population modeling and Bayesian estimation for the deconvolution of blood alcohol concentration from transdermal alcohol biosensor data

PDF

A nonlinear pharmacokinetic model used in calibrating a transdermal alcohol transport concentration biosensor data analysis software

PDF

Determining blood alcohol concentration from transdermal alcohol data: calibrating a mathematical model using a drinking diary

PDF

Prohorov Metric-Based Nonparametric Estimation of the Distribution of Random Parameters in Abstract Parabolic Systems with Application to the Transdermal Transport of Alcohol

PDF

Linear quadratic control, estimation, and tracking for random abstract parabolic systems with application to transdermal alcohol biosensing

PDF

M-estimation and non-parametric estimation of a random diffusion equation-based population model for the transdermal transport of ethanol: deconvolution and uncertainty quantification

PDF

Physics-informed machine learning techniques for the estimation and uncertainty quantification of breath alcohol concentration from transdermal alcohol biosensor data

PDF

Estimation of random input to semi-linear abstract parabolic systems with application to quantitative description of drinking behavior based on trans-dermal alcohol concentration

PDF

A Bayesian approach for estimating breath from transdermal alcohol concentration

PDF

Delta Method confidence bands for parameter-dependent impulse response functions, convolutions, and deconvolutions arising from evolution systems described by…

PDF

Simultaneous parameter estimation and semi-blind deconvolution in infinite-dimensional linear systems with unbounded input and output

PDF

Covariance modeling and estimation for distributed parameter systems and their approximations

PDF

Functional connectivity analysis and network identification in the human brain

PDF

From least squares to Bayesian methods: refining parameter estimation in the Lotka-Volterra model

PDF

A comparison of GLM, GAM, and GWR modeling of fish distribution and abundance in Lake Ontario

PDF

Nonparametric estimation of an unknown probability distribution using maximum likelihood and Bayesian approaches

PDF

Reconstruction and estimation of shear flows using physics-based and data-driven models

Asset Metadata

Creator Li, Jian (author)

Core Title Distributed parameter model based system identification and filtering in the estimation of blood alcohol concentration from transdermal alcohol biosensor data

School College of Letters, Arts and Sciences

Degree Master of Science

Degree Program Statistics

Publication Date 07/17/2017

Defense Date 06/22/2017

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag blind deconvolution,distributed parameter systems,filtering,OAI-PMH Harvest,system identification,transdermal alcohol biosensor

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Rosen, Gary (committee chair), Luczak, Susan (committee member), Wang, Chunming (committee member)

Creator Email jli981@usc.edu,silencer1127@gmail.com

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c40-399407

Unique identifier UC11264510

Identifier etd-LiJian-5519.pdf (filename),usctheses-c40-399407 (legacy record id)

Legacy Identifier etd-LiJian-5519.pdf

Dmrecord 399407

Document Type Thesis

Rights Li, Jian

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA