Physically structured sequential data modeling: integration of qualitative and quantitative research

Abstract

However broad the fuzzy methods acceptance may be in the industry, mastering of crisp methods for finding optimal solutions to complex real-life problems, as well as for reducing the computational demands of their implementation is an absolute necessity for the intelligent systems modeling. Given the existence of two different approaches to research (qualitative vs. quantitative), this paper updates on the problem of searching the sufficient and adequate model for integrating of qualitative (fuzzy) and quantitative (crisp) methods in science and practice of data engineering.

1. Introduction

Fuzzy sets and fuzzy technology have been applied extensively to many areas where classical modeling, reasoning, and computing cannot be made deterministic and unambiguous. As we say: with the lapse of time, the plot deepens [1]. Lotfi Zadeh: “As the complexity of a system increases, our ability to make precise and yet significant statements about its behavior diminishes until a threshold is reached beyond which precision and significance (or relevance) become almost mutually exclusive characteristics.”

However, no generally accepted exact definition of complexity exists yet. In the authors’ opinion, this idea means complex systems’ behaviors exhibit uncertainty as of the very first property. Uncertainty is what accounts first of all for the complexity of the system. This property may manifest itself in various forms: stochastic vs possibilistic, recoverable vs unrepairable, estimable vs unevaluable, stochastically vs pre-existing, isolated vs spatially organized. Consequently, there are many ways to act under uncertainty.

In research and practice, there are two critical approaches used for knowledge acquisition: qualitative vs quantitative. As the first one is typical for fuzzy sets theory and technology, the latter remains an exclusive prerogative of crisp science. Notwithstanding these methods appear to be very different in many attributes, the time has come when the integration of them is needed. Increasingly popular fuzzy technology in the industry (FTI) and other human activities requires that qualitative data and decisions can be coded and reported quantitatively [2]. Various models of such integration are possible. As Elena T. Carbone, who is Associate Professor and Chair, Department of Nutrition, UMass Amherst, illustrates, essential models are: (i) Qualitative methods are used to help develop quantitative measures and instruments. (ii) Qualitative methods are used to help explain quantitative findings. (iii) Quantitative methods are used to embellish a primarily qualitative study. (iv) Qualitative and quantitative methods used equally and in parallel.

In research, we always have to fit our questions and expected results with methods. For example, the human body daily temperature, without regard to the measurement site, is a non-linear process characterized by moment-to-moment complex variability and even stochasticity [3]. Human body temperature indications are used in ambulance and clinical practice on a full scale to make correct health survey and good treatment decisions. Daily local thermometry data must be evaluated not only in such fuzzy terms as “growth,” “precarity,” “equilibrium,” or “decline” a physician can see firsthand but be supplemented by quantitative information including results of more sophisticated analysis and predictions.

However, the teaching of “Complex systems modeling” as an advanced course of studies for master’s degree in Information Systems reveals one negative situation: Classic methods of mathematical modeling are full of shallow, quick fixing, subliminal words for those students who shortly before were taught an intensive course of fuzzy methods. It seems to be a part of a general incompatibility factor typified by knowledge representation mismatch in cognate disciplines. It is as if facilitators of fuzzy methods, on the one hand, and instructors of applied mathematics, on the other, are two competing fractions of the same scientific research community. We stand witness of today’s prevailing commonsense assumption that fuzzy methods and technology, in and of itself, will solve the serious challenges we face in artificial intelligence, automation, data analytics, robotics, advanced computing power and the Internet of Things. This thinking seems to be not only wrong but harmful, particularly to those as-yet-malleable hearts and minds.

In our research, we strive to fix the problem of searching the effective and adequate model of integrating qualitative (fuzzy) and quantitative (crisp) methods in teaching science and practice of data engineering. Keeping to the subject, we organize our paper in the following order. Section 2 provides two illustrative examples of data modeling task where crisp mathematics is a vital necessity. Section 3 phrases the problem in three precise classificatory generalizations as they may sound in physically motivated data fitting or time series analysis. Sections 4 and 5 present the critical research insights into solutions for the examples of Section 2. The last section concludes the paper with some recommendations for the instructors and researchers in the field.

2. Exemplifying the Problem Example 1

Suppose the ore agglomerate analyzer we operate with has indicated the presence of n fissionable materials in the in-situ rock, however at an unknown ratio. Material kinds in this ore agglomerate may be known (case C1) or hidden (case C2). Geiger-Mueller counter readings yk , y(kτ ) taken at scheduled times tk , kτ, k = 1, 2, . . . , m , m > n with the constant sampling interval τ are the only data available to determine the mass mixing ratio (in C1) and also ingredient kinds if needed (in C2).

To get things rolling, we write the radioactive decay law for a kind of material as M (t) = M0 exp (−λt) with M0 being the initial mass. With d , exp (−λτ ) , we have Mk , M (tk) = M0dk . Exclusively for definiteness in application to our example, let n be equal 3. Then we obtain yk ≈ M0,1d1k + M0,2d2k + M0,3d3

k vk , yk − x1d1k + x2d2k + x3d3

k k = 1, 2, . . . , m > 3 ( 1 ) with vk defined as a discrepancy between the G-M readings and the general decay law. This value is also termed as a residual in data fitting practice. ConTsidering vector x = x1 T· · · xn T ∈ Rn (with n = 3 for this example), vectors y = y1 · · · ym and v = v1 · · · vm , both in Rm , and also (m × n) -matrix H obvious in view of its transposed form H T =  ddd231 ddd223212 ··· ··· ··· ddd231mmm  , we avail of the opportunity to write

y = Hx + v in the standard linear algebra notation.

NB 1: In experimental design theory, matrix H is known as the experiment design matrix.

Three unknown model parameters x1 , M0,1 , x2 , M0,2 , and x3 , M0,3 gathered into a column vector x are to be estimated from the experimental data y for either of the two cases: C1 and C2: C1 Ingredient kinds and hence their half-life periods are known. Thus, we infer di , exp (−λiT ) to be given as some real numbers 0 < di < 1 , i = 1, 2, 3 . This case is a trivial one (linear regression analysis). We use the three-variable model ( 1 ) written as ( 2 ) to obtain the off-line Ordinary Least Squares (OLS-optimum) estimator xˆ = (xˆ1 xˆ2 xˆ3) from the normal equations:

HTH xˆ = HTy.

C2 This case is far more complex. The above OLS solution turns out to be impossible since the observation matrix H is now inherently ambiguous. The task becomes nonlinear, and the question of ‘How to find an efficient solution for C2?’ deserves to be asked.

Example 2

For another example, consider the human body daily temperature data (HBDTD) as having a periodic nature depending on the time of day (figure 1). ( 2 ) ( 3 ) ( 4 ) ( 5 )

The biophysical law dictating such a circadian behavior is fundamentally unknown, however looking at figure 1, a researcher may decide to consider the process approximately as generated by a deterministic harmonic oscillator (DHO) around θ? , the human average daily temperature. Starting from this point of view, he/she approximates the thermometry data yk , y(kτ ) taken at m scheduled times tk , kτ, k = 1, 2, . . . , m by the following expressions: yk − θ? ≈ AN sin(kωNτ + φN) = aN sin(kωNτ ) + bN cos(kωNτ ) vk , (yk − θ?) − (aN sin(kωNτ ) + bN cos(kωNτ )) k = 1, 2, . . . , m > 2, aN = AN cos φN, bN = AN sin φN Now vector y ∈ Rm appears in equation ( 2 ) with the (m × 2) -matrix

 sin(1ωNτ ) H =  sin(2ωNτ )  . . .

sin(mωNτ ) cos(1ωNτ )  cos(2ωNτ ) 

. . .  cos(mωNτ ) times vector xT = x1 x2 of n = 2 unknown quantities: x1 , aN and x2 , bN .

NB 2: In the above expressions, subscript N stands for marking the value as “natural” one.

NB 3: This example corresponds to case C1 if the natural circular frequency ωN , 2π/TN can be taken as a known steady-state value, for example by letting TN be equal 24 hours (i. e. 1440 min). However, in some of the most complicated cases, ωN , 2π/TN is by no means known quantity. This situation brings us back to the nonlinear case C2.

NB 4: Use of the DHO seems to be an oversimplistic, too easy approach. A more realistic model of the human body circadian temperature rhythm should include, as a minimum, a random component in observations—for example, in the form of Wiener stochastic process (Brownian motion [4, p. 151])—and also replace the DHO by a Gaussian harmonic oscillator (GHO [4, pp. 167, 169]). It is the case which was given in research and to which a full and thorough examination has been made [5].

3. Classificatory Problem Generalizations (PGs)

• PG-1 The first generalization is plain to see:

How to move from the offline OLS-optimum estimator ( 3 ) to the online one? This question is not a problem at all for the linear case C1. Recursions for the OLS-optimum estimator are well-known either as the Information Processor or Covariance (Kalman Filter, KF) Data Processor [6]. However, when having to deal with the nonlinear case C2, we will be forced to keep PG-1 while also moving to the next problem generalization. • PG-2 The other generalization lies in answering the question:

How to justifiably represent the data source (DS in figure 2) whose output y in ( 2 ) is the observed data to be processed to obtain the OLS-optimum estimator xˆ for the model parameter x— the notations used up to this moment including figure 2 implies that we take them in the sense of relations ( 1 ) to ( 5 )—by a dynamical (state-space) gray-box (physically structured) model (DGBM, eventually stochastic)? We mean to describe that DG-BM by the following two discrete-time equations: xk = Φ(θ)xk−1 + Bd(θ)uk−1 + Gd(θ)wd,k−1 yk = Hxk + vk k = 1, 2, . . . ( 6 ) taking into account that from this time on we are switching to the notations peculiar to the classic Stochastic Models, Estimation and Control (SMEC) theory [4], as well as to the modern Adaptive System Identification (ASI) literature: [7], [8], and [9].

This tradition means that we start considering PG-2 from the assumption that a modeling technique has produced an adequate data description in the form of a linear stochastic difference equation— the first one in ( 6 )—to describe the propagation of state x(·) , with discrete-time noise corrupted linear measurements—the second equation in ( 6 )—classified as the available experimental (or reallife) data. Here x(·) , x(t) is an n -vector state process, one sample of which would generate a system state time history: xk−1 would be the system state at time tk−1 propagating in the direction to xk through an n -by- n state transition matrix Φ(θ) ; u(·) , u(t) is an r -vector of piecewise continuous deterministic (premeditated) control input, one sample uk−1 of which would act upon the state through an n -by- r deterministic input matrix Bd(θ) at times between tk−1 and tk . In summary, ( 6 ) is a gray-box-model parameterized by an unknown p -vector θ , with the usual Gaussian description of the initial state and zero-mean white noise processes of s -vectors wd,k−1 entering the model through an n -by- s noise input matrix Gd(θ) , and of m -vectors vk acting as

Data y DS

Model

x

Hx

v = y − Hx + − mxin kvk2 random measurement errors:

E {x0} = x¯0(θ),

E n[x0 − x¯0(θ)][x0 − x¯0(θ)]To = P0(θ) E n

T o = Qd(θ)δij , wd,iwd,j

E nvivjTo = R(θ)δij δij , • PG-3 The third generalization implies answering the question:

How to identify the parameterized gray-box state-space model ( 6 )–( 7 ) with a minimum of computational effort and acceptable quality level as regards the unknown vector parameter θ dictated by the laws of physics governing the data source?

The following text is aimed at seeking answers to the research questions declared in PG-2 and PG-3 as supported by Examples 2 and 3.

4. The solution of PG-2 and PG-3 for Example 1

Having M (t) = M0 exp (−λt) we wonder: What is the differential equation whose solution is M (t) ? Because M˙ (t) = −λM (t) , an immediate answer for PG-2 follows: M˙ (t) + λM (t) = 0 with M (0) = M0 . For model ( 6 ), we have: both uk−1 and wd,k−1 vanish, Φ(θ) = diag [d1 d2 d3] , di , exp (−λiT ) , with θ , [d1 d2 d3]T and H = [1 1 1] . For solving PG-3 by our Active Principle of Adaptation (APA) method [9], we move from this physically data model (PhDM) to the standard observable data model (SODM) by the similarity transform x? = W?x with the observability matrix W? , HT (HΦ(θ))T (HΦ(θ)2)T T which is the third-order Vandermonde matrix V3 = V3(d1, d2, d3) , det V3 = (d2 − d1)(d3 − d1)(d3 − d2) in this example. We find Φ?(θ) = W?Φ(θ)W?−1 to be the SODM state transition matrix

 0 Φ?(θ) =  0

? −a0 1 0

? −a1

 −10a2?  ,  −−aa1?0? == d(d11dd−213d(dd23d)2(13+d3d−3)d1) + (d2d−23d(d1d)1(32+d3d−3)d2) + (d1d−33d(d3d)1(33+d3d2−)d2)

 −a2? = (d2−d1)(d3−d1) + (d1−d2)(d3−d2) + (d3−d1)(d3−d2) Once Φ?(θ) is known to be the companion matrix for characteristic equation q(λ) , det[Iλ − Φ?(θ)] , in this case q(λ) = λ3 + a?2λ2 + a1λ + a0? , and given the fact that q(λ) does not alter after a similarity ? ( 7 ) transform, i. e. q(λ) , det[Iλ − Φ(θ)] , in this case q(λ) = (λ − d1)(λ − d2)(λ − d3) , it is not hard to get the simplest form of parameters in ( 8 ) (that can be readily verified to be true by straightforward reducing ( 8 ) to ( 9 )): The APA-method has been proved [9] to be computationally efficient in providing estimates for a ? ? companion matrix, in this case estimates aˆi , i = 0, 1, 2 for parameters ai , i = 0, 1, 2 in ( 9 ). As a result, one might substitute aˆi? for ai? in ( 9 ) to extract estimates θˆi as a solution to ( 9 ) for the sought parameters di , θi , i = 0, 1, 2 . It is fortunate that this impractical way is not necessary to do (in this example). Indeed, the task of finding parameters di from system ( 9 ) should be replaced—when Φ(θ) = diag [d1 · · · dn] —by solving polynomial equation q(λ) = (λ − d1) · · · (λ − dn) = 0 ; this computational procedure is tried-and-true in many math problem solvers. After one has found (estimated) di , matrix H for system ( 2 ) may be considered known, normal equations ( 3 ) prove to be linear and can be solved either offline or preferably online (i. e. sequentially), see PG-1 in Section 3 for more details.

5. The solution of PG-2 and PG-3 for Example 2

NB 5: In HBDTD stochastic model construction, we follow one of our recent research [5].

Let HBDTD be modeled in continuous time by an ωN = 2π/TN rad/min, TN = 24 h deterministic harmonic oscillator [x1(t), x2(t)]T whose output x1(tk) = aN sin(kωNτ ) + bN cos(kωNτ ) occurs in the first line of ( 4 ) as disturbed by an Ornstein–Uhlenbeck process x3(t) , i. e.

dx1(t) = x2(t) dt dx2(t) = −ωN2x1(t) dt dx3(t) = −(1/T )(x3(t) − θ?) dt + σp2/T dβ˚(t)     physical state equations with unknown parameters T > 0 and σ > 0 and β˚(t) being the standard (zero-mean and unitdiffusion) Wiener process, t0 →lim−∞ β˚(t0) = 0 (a.s.) , where θ? is a human average daily temperature. Let further the current body temperature x1(t) + x3(t) be measured every τ minutes with a sensor whose accuracy is limited by an error v(tk) ; process v , {v(tk); t ∈ Z} is modelled by a discretetime white noise with a covariance R : y(tk) = x1(tk) + x3(tk) + v(tk) measurement equation

The various equivalent state representations can be related through similarity transformations to the original model ( 10 ) and ( 11 ). One result as it is in [11] follows: 3dDRCM, 3 -dimension ( n = 3 ) Discrete-time Real-valued Canonical (Jordan) Model  cos ωNτ xk+1 =  sin ωNτ 0 − sin ωNτ cos ωNτ

0 |  +  | } xk + vk, k = 1, 2, . . .

(14) with the sampling interval τ (we assume that τ = 5 min ), λ , 1/T , d , e−λτ . In (13) and (14), wd,k is a discrete-time zero-mean covariance Q = 1 white noise; ut = θ? ; ωN = 2π/TN rad/min, TN = 24 h ; vt is a zero-mean covariance R measurement error we assume R = (0.1)2 , and vector θ , [ λ σ ]T consists of two unknown parameters λ and σ . For any arbitrarily selected experiment, θ takes on the value ˚θ= [ ˚λ ˚σ ]T , which is to be identified.

To identify the unknown value of parameter θ , we base our solution on using the APA method reported theoretically and computationally in [5]. To validate the results, we plan the experimental conditions as shown in Table 1 [5].

Given any fixed (true) system parameter ˚θ 6= 0 , the system was simulated for 2N measurements. The data set was divided into two segments, each representing one day. We used the first segment for identification and the second for validation. The validation metric chosen was the variance-accountedfor (VAF) index defined as follows [12]:

VAF = max 1 − var(yk − yˆk) , 0 var(yk) × 100% (15) where yk was the validation data, yˆk was the estimated output signal; var(·) denotes the variance of a quasi-stationary signal.

Results demonstrated by figure 3 and summarized in Table 2 show that the computed estimate θˆ comes close to the actual parameter value ˚θ thus providing an adequate model approximation of the output signal yk by the estimated output yˆk . As figure 3 suggests, our HBDTD stochastic model (13), (14) is capable of capturing the human body temperature circadian variations after the unknown model parameter ˚θ has been estimated based on the first segment data—with our identification method—and then replaced by the finally computed estimate θˆ to compare estimates predicted for the next day, i. e. yˆk with the second segment data, i. e. yk for computing VAF as defined by (15).

NB 6: All the works related to computational issues, computer programming the framework and experimental validation of the method with the simulated data in the abovementioned facts have been contributed by professors Julia V. Tsyganova of Ulyanovsk State University, Maria V. Kulikova of CEMAT, Instituto Superior Te´cnico, Universidade de Lisboa, and Andrey V. Tsyganov of Ulyanovsk

Human body temperature level [ oC] ˚θ= [˚λ , ˚σ]T, the true value of θ θˆ = [λˆ , σˆ]T, the finally computed estimate Relative estimation error δ = ||θ − θˆ||/||θ||

˚ ˚ VAF State Pedagogical University, the co-authors of [5]. The research was supported by professor Andrey B. Peskov of Institute of Medicine, Ecology and Physical Education, Ulyanovsk State University on behalf of Ulyanovsk Regional Hospital.

Recurring to the former subject declared in Example 2, Section 2: How to identify the unknown quantities aN and bN for model ( 4 )? , we come to the following solution to the question.

Note that in this case expression ( 2 ) assumes a slightly modified shape

y = H x + xˆ3 + v with H from ( 5 ), x , abNN , xˆ3 = xˆ3(1τ ) xˆ3(2τ ) · · · xˆ3(mτ ) T where xˆ3(kτ ) is the estimator for the Ornstein–Uhlenbeck process x3(t) at t = kτ , k = 1, 2, . . . , m . It is fair to say if take into account the physically structured model given by ( 10 ) and ( 11 ). Those estimates are produced by the KSRCF (cf. figure 3) using the HBDTD model with θˆ = [λˆ , σˆ]T , the finally computed estimate substituted for the unknown ˚θ = [˚λ , ˚σ]T , the true value of θ (cf. Table 2). As a result, to obtain the OLS-optimum estimator xˆ = [xˆ1 xˆ2]T for x , aN , one has to solve the normal equations whose bN general form ( 3 ) has to be modified, in the light of (16), to equation (17): xˆ = H T(y − xˆ3). (16) (17) • Examples presented in the paper show that finding the optimal solutions to complex problems cannot be distilled only to the qualitative approaches peculiar to fuzzy contexts. • Increasingly, we come to the conclusion that the quantitative methods must be used to get a deeper insight into a primarily qualitative—fuzzy methods based—study. • A researcher, especially at the beginner level, must be educated to use qualitative and quantitative methods convergently. • The integration of qualitative and quantitative methods as exemplified in this paper seems to be instructive to research strategy for physically structured sequential data modeling.

Acknowledgments

This research is funded by Russian Foundation for Basic Research and the government of the Ulyanovsk region of the Russian Federation, grants 18-47-730001 p a and 18-41-732002 p mk.

References

[1] Zimmermann H-J ( Hans-Ju¨rgen) 2001 Fuzzy Set Theory- and Its Applications ISBN 978-94-010-3870-6 ISBN 978-94- 010-0646-0 (eBook) DOI 10 .1007/ 978 -94-010-0646-0 (Springer Science+Business Media, LLC) -4th ed

[2] Carbone E T 2014 Using Qualitative & Quantitative Research Methods to Answer Research Questions (Webinar for UMass Medical Center , Worcester, MA)

[3] Kelly G 2006 Body Temperature Variability (Part 1): A Review of the History of Body Temperature and its Variability Due to Site Selection, Biological Rhythms , Fitness, and Aging Alternative Medicine Review vol 11 no 4 pp 278 - 293

[4] Maybeck P S 1979 Stochastic Models , Estimation, and Control vol 1 (Academic Press New York-San Francisco-London)

[5] Semushin I V , Tsyganova J V, Kulikova M V , Tsyganov A V , Peskov A B 2016 Identification of Human Body Daily Temperature Dynamics via Minimum State Prediction Error Method Proc . 2016 European Control Conference (ECC) (Aalborg: Denmark ) pp. 2429 - 2434

[6] Bieman G J 1977 Factorization Methods for Discrete Sequential Estimation (Academic Press New York-San FranciscoLondon)

[7] Merce `re G , Prot O , Ramos J A 2014 Identification of Parameterized Gray-Box State-Space Systems: From a Black-Box Linear Time-Invariant Representation to a Structured One IEEE Trans . on Automat. Contr . vol 59 no 11 pp 2873 - 2885

[8] Semushin I V 1985 Identification of Linear Stochastic Objects from Incomplete Noisy Measurements of the State Vector Automation and Remote Control (The USSR Academy of Sciences ) vol 46 no 8 pp 975 - 985

[9] Semushin I V 2011 Adaptation in Stochastic Dynamic Systems-Survey and New Results II Int . J. Communications , Network, and System Sciences vol 4 no 4 pp 266 - 285

[10] Semushin I V 2011 Computational Methods of Linear Algebra and Estimation (Ulyanovsk : Ulyanovsk State Technical University Publishers) Access mode: http://venec.ulstu.ru/lib/disk/2013/119.pdf

[11] Semushin I V , Tsyganova J V, Skovikov A G 2013 Identification of a Simple Homeostasis Stochastic Model Based on Active Principle of Adaptation Proc. 15th Applied Stochastic Models and Data Analysis Int. Conf. (ASMDA- 2013 ) (Mataro (Barcelona): Spain ) pp. 775 - 783

[12] Verdult V , Verhaegen M 2002 Subspace Identification of Multivariable Linear Parameter-varying Systems Automatica vol 38 pp 805 - 814