CC BY
1
Math-Net.Ru
I. V. Беш^Ыи, Идентификация передаточной функции посредством минимизации рассогласования оценок состояния адаптивного и оптимального фильтров, Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки, 2023, номер 3, 544-572
001: 10.14498^^2037
Использование Общероссийского математического портала Math-Net.Ru подразумевает, что вы прочитали и согласны с пользовательским соглашением
http://www.mathnet.ru/rus/agreement
Параметры загрузки:
IP: 109.252.33.182
29 сентября 2024 г., 12:14:08
Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki
[J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2023, vol. 27, no. 3, pp. 544-572 ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu2037
MSC: 93E35, 93A30, 65C20
Transfer function identification by minimizing the adaptive vs. optimal filter state estimates mismatch
I. V. Semushin
Ulyanovsk State University,
42, L. Tolstoy st., Ulyanovsk, 432017, Russian Federation.
Abstract
The article is concerned with a further development of the Active Principle of parametric system identification in the class of linear, time-invariant, completely observable models. As the identification target model, the optimal Kalman filter (OKF) is designated that is present, no more than conceptually, in the system's discretely observed response to a training excitation of the white noise type. By modifying the physically given structure into the standard observable model in both the observed response and the Adaptive Kalman Filter (AKF), a so-called Generalized Residual (GR) is constructed equaling the mismatch between the adaptive and the optimal filter state estimates plus an AKF-independent noise component. By virtue of this modification, the GR mean square becomes a new model proximity criterion for these filters. Minimizing this criterion via conventional practical optimization methods produces exactly the same result (AKF = OKF) as would be obtained by minimizing the theoretical criterion being, unfortunately, inaccessible to any AKF numerical optimization methods. The article presents a detailed step-by-step procedure explaining the above solution in terms of a parameterized transfer function. For the sake of clarity and for stimulating real world applications of the approach, the article employs the transfer function model of a twisted-pair line in a typical xDSL system. The implementation challenges of theoretical provisions of the method are discussed. The issue of extending the proposed approach to the problems of identifying linear models for nonlinear systems is outlined in the directions for further research.
Keywords: LTI model, complete observability, Kalman filter, adaptive filter, indirect performance index, implementation challenges.
Received: 23rd June, 2023 / Revised: 5th August, 2023 / Accepted: 21st August, 2023 / First online: 28th September, 2023
Mathematical Modeling, Numerical Methods and Software Complexes Research Article
© Authors, 2023
© Samara State Technical University, 2023 (Compilation, Design, and Layout) 3 ©® The content is published under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/) Please cite this article in press as:
Semushin I. V. Transfer function identification by minimizing the adaptive vs. optimal filter state estimates mismatch, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2023, vol. 27, no. 3, pp. 544-572. EDN: XUQAGA. DOI: 10.14498/vsgtu2037. Author's Details:
Innokentiy V. Semushin https://orcid.org/0000-0002-3687-1110
Dr. Techn. Sci., Professor; IEEE Member; Professor; Dept. of Information Technology;
e-mail: [email protected]
1. Introduction. The theory and practice of system identification (SI) in their more than half a century of history have received a powerful development reflected in hundreds of thousands of scientific publications around the world. As Gianluigi Pillonetto and Lennart Ljung note in their recent paper [1], 'Despite its long history, such research area is still extremely active.' Indeed, even in a nonlinear setting, research is being done on how to deal with the presence of nonlinear distortions in systems by using linear SI techniques [2].
The abundance of publications in this field signaled the need for some serious cleanup work in order to single out the truly independent concepts. According to Ljung, in SI there are two independent and universal key concepts: the choice of a Parametric Model Structure, PMS, and the choice of a Model Proximity Criterion, MPC, the latter is the criterion of fit indicating erroneousness of a model with respect to a target [3]. Looking generally at what takes us in the identification process from observed data to a validated model, there are four main components: ' (1) The data itself, (2) The set of candidate models, (3) The suitability criterion, and (4) The validation procedure' [4].
Indeed, at the heart of SI—or, in the modern AI terminology, of system mathematical model machine learning—is the principle of fitting the response data of an adaptive predictive model to the data of the real system response, which actually exists as a 'black box,' under conditions of the same excitatory (learning or training) input for them, by some predefined cost function. Nevertheless, the question of interest remains: Given the PMS, HOW to use the available data to predefine the MPC?
In the SI community, the impressive Prediction Error Framework, PEF, [5] reflects a generally accepted understanding of this issue. Such a view has been expressed [6] on more than one occasion: 'All existing parameter identification methods can be seen as special cases of this prediction error framework.' At that, existing PEF methods fit the adaptive model in the system response space, not in the state space. This is due to the fact that the useful concept of 'state space' is intended for purely theoretical work to formulate and minimize the system model optimality criterion, which we call direct performance index, DPI. This limiting feature is generated by the certainty that it is impossible to overcome the obvious barrier, namely, the inaccessibility of state space elements in explicit form and, hence, of DPI. In fact, DPI cannot be accepted as MPC for identification algorithms.
Putting this barrier overcoming on the agenda, this article proposes an alternative solution to the HOW question posed above. That is, in formulating the research question, the intention here is to form an Indirect Performance Index, IPI, and then organize the minimization of IPI so that it is equivalent to minimizing the discrepancy between the internal states of the adaptive model that is available, and the internal states of the optimal model, which is only theoretically known as Optimal Kalman Filter, OKF, but is not accessible because of parameter uncertainty and, moreover, is sought as the target result of the parametric optimization of the Adaptive Kalman Filter, AKF.
Thus, the alternative approach considered in this paper should, as is conceivable, minimize the discrepancy between the AKF and OKF state estimates. This would be reasonable, since the notion of 'state' is intended to exhaustively characterize the behavior of an object. Moreover, such minimization, if implemented in practice, corresponds one-to-one to a theoretically optimal filter designing. This
feature prevents the SI algorithm deviating from theoretical results of the OKF design and, therefore, from the bias errors inherent in some other SI methods.
To make such a solution feasible, the real system observed output is represented as if it were generated by the desired but hidden from us optimal filter rather than the given physically structured system. In the interest of realizing such a conceptual vision, the system is supposed to be completely observable to gain access to the ability to change the basis of the system's internal states formally without changing its input-output description, which is called the transfer function, TF, in the class of linear time-invariant, LTI, systems the article addresses to.
Building an LTI model for a dynamic system is usually being made either in the frequency domain (by a TF) or in time domain (by differential or difference equations) [7] to answer the challenge of reducing model uncertainty. The LTI SI theory and practice have reached a high degree of maturity and are frequently used in many disciplines where an object of interest exists, for example, in mechanical [8], electrical [9], electronic [10], chemical [11], civil [12], and even in biomedical [13,14] applications. Besides, the point is that the object of interest for which it is necessary to parameterize the model in the form of TF can be not real, but fictitious. The most striking example of this is the construction of a dummy filter forming a model stationary random process from a white-noise process, which should approximate by its correlation function the experimental correlation function of a real process in a real system. An example may be identification of a parameterized instrumental errors model of a multi-component inertial navigation system [15].
The novelty of this article is that it encourages the application of the approach in the real world where it has not previously been considered. As an original example and for clarity, it uses the twisted-pair model in a typical Digital Subscriber Line, xDSL system [16]. As known, there exists a computational cost reduction challenge to solve the crosstalk precoding problem and this problem cannot be solved without knowing the direct and cross channel TFs, DCTFs and CCTFs [17]. There are several solutions to this engineering task in the literature of recent years, to exemplify [18-23]. Most of them are similar in that they propose to estimate channel TFs in the frequency domain, which is quite understandable since TF itself is a function of signal frequency and is to be known for each tone to eliminate the crosstalk phenomenon. Such methods of TF estimation narrow the field of their possible application, reducing it to the xDSL technology, where they are recognized to be effective. In contrast, this work proceeds from the fact that the problem of TF estimation can be solved in a more general formulation, considering it as a problem of parametric LTI model identification for a dynamic system in the state space time domain.
This article addresses the following research issues.
® First, we intend to overcome the obvious problem that the state vector of a dynamical system is unattainable explicitely, or, put this differently, is beyond of reach in a perfectly measured form, just as a signal disturbed by noise in filtering problems is, by definition, immeasurable in a pure form. The same applies to optimal state estimators, since the optimal filter remains machine-unrealizable or, put it tentatively, covert in the mesurement data until the necessary parameters are identified.
® The solution to overcome this unattainability barrier in [24] considered indi-
vidual cases of a'priori uncertainty level constraints under which the solution works. We aim to show that this solution is in fact feasible with no limitations on the size of a'priori parametric uncertainty of the system model. We must make sure that this method of solution makes it universally applicable under the only unencumbered condition: the LTI model under study is completely observable and can be considered adequate to reality. We want to show that this quality of solution is achievable by converting the model into a standard observable form, SOF, to gain a solution in a computer-implementable tool.
® As for common xDSL applications, we have to check whether it is possible to use time-domain formulas instead of traditional frequency-domain formulations to estimate the DCTFs or CCTFs, and show how to do so for any frequency (or tone) of interest in the channel operating frequency range.
® Further, to organize the computational process with its numerical robustness and also to translate all decisions into a software design, reasonable suggestions are needed.
© Finally, a determination has to be made about the novelty of this work in terms of its results, advantages, and limitations, and concerning objectives of further research in the proposed direction.
Consideration of these issues constitutes the main content of this article. Section 2 is devoted to an illustrative example for which the IPI-based LTI system identification method may be of practical interest. Section 3 presents a formal statement of the problem with two generalizations. Section 4 explains a detailed procedure of how to identify a parametric OKF estimator in terms of a parametrized TF. Section 5 discusses three practical challenges associated with implementing the solution: (1) organizing the computation time; (2) ordering the computation in terms of its numerical robustness; and (3) scheduling the work for a software project. The final Section 6 summarizes the work, describes the limitations, and outlines possible research on the approach.
2. An illustrative example. Only within this example, symbol f is to designate the signal frequency in the electronic RsLsGsCs-circuit of Fig. 1 that mimics a very short—of length Al—section of a twisted-pair line in the typical xDSL system [17, Chapter II]. The circuit can help the DCTF evaluation for transmission line of I full length. Primary transmission line parameters are R = R(f), L = L(f), G = G(f) and C = C(f) being functions of frequency f can be seen as expressed through the secondary cable parameters for standard twisted pairs depending on cable diameter, material and design. Taking these values from [17, p. 19] yields the parameters of a short section: the section resistance Rs = R-Al, the section inductance Ls = L-Al, the section conductivity Gs = G-Al, and the section capacitance Cs = C-Al.
Remark l.The subscript s written in roman typestyle for lowercase index should not to be confused with the below Laplace variable s. It serves to remind that the quantity refers to a twisted pair section as shown in Fig. 1. The square brackets below denote the dimensionality of physical quantities in units of SI.
The formula called Generalized Ohm's Law (GOL) in the complex domain defines the impedance of an electronic two-terminal element as across variable (voltage) divided by through variable (current), both in terms of the Laplace transform. When it is coupled with Kirchhoff's Current and Voltage Laws (com-
monly shortened to KCL and KVL), one has a sufficient set of tools for analyzing circuits. By writing KCL and KVL for the circuit in Fig. 1 (a), the equivalent linear Two-Port Network (TPN) shown in Fig. 1 (b) is obtained, yielding the following equations (1) in terms of the Laplace transform variables:
■ V1 ' \M*) B(s)] 'V2'
h C (s) D(s)_ .J2.
(1)
A(s) â 1 + Fs(Sy, [A(s)j = 1,
Fs(s) â (Rs + sLs)(Gs + sCB); [FB(s)] = 1,
B(s) â (Rs + SLS); [B(s)] = Q,
C(s) â (Gs + sCs); [C(s)] = S = Q-1,
D(s) â 1.
Remark 2. In this notation, referring to variables and their transforms interchangeably, Laplace transforms are distinguishable by the use of an uppercase letter or (in more detail) the complex-valued argument (s = a +jw) with u, [w] =Rad/s, equaling angular velocity 2nf and f, [f] =s-1, meaning the frequency variable,
j = V=1.
Define the dimensionless TF of a AMength line section as Ts(s) = V2/Vq. Applying KVL yields Vq and KCL h in (2)
Vo /1
+ [Zo + (Rs + sLB)]h V2/Z2 + (Gs + sCs)V2
(2)
12 13
It follows that Ts(s) is the quantity inverse to
l + [Zc + (Rs + sLs)] [Z-1 + (Gs + sOO] . Together with (1), it leads to (3)
Ts(s) = {A(s) + Zo [Z2-1 + C(s)] + Z-1B(s)}
1
(3)
From now on, consider a mathematically idealized experiment with, condition (i), a voltage source V0 connected to the section input thus assuming Z0 ^ 0;
and with, condition (ii), a voltmeter having a very large inner impedance Z2 ^
connected to the section output to measure V2. In this scenario, Ts(s) ^ Ts(s) = 1 °
A- ( s). Ts(s) is the section intrinsic transfer function, SITF, we are interested in to move closer to the reality of xDSL multi-user transmission, xDSL-MUT.
As known from [17, Chapters II and III], the DCTF denoted by H(f, I) is frequency-dependent and changes with the cable length . When the transmission line is connected to a source Vs with source impedance Zs and terminated with load impedance Z^, this H(f, I) is expressed by (4)
h (f, i) = Zs + Zl
( Zs + ZL )cosh(jl) + Z*smh(jl) Z* = Z +
A ^ , ZS • ZL
(4)
Z
through the characteristic line impedance Z* defined as
Z a [R + №fL
= y GTWfC (5)
and the propagation constant 7 = 7(f) calculated by
7 (f) = V(R + №fL)(G +]2vfC). (6)
If the line terminates ideally at Z* (5), so that Zl = Z = Zs, the channel transfer function simplifies [17, Chapters II and III] to
H(f, l)= e-^K (7)
Now, noticing a similarity between (6) and Fs(s) in (1), we obtain Fs(s) = (Al)272(s) after substitution s =]2kf. Hence
7 (f ) = m-1^[Ts(s)]-i\s=i2,f -1. (
s=]2nf
If one manages to evaluate the expression under the square root sign in (8) as a complex-valued magnitude in dependence on frequency , then the problem is solved for any value desired. Thus, the solution is to parametrically identify the SITF, that is, Ts(s) to use it in (8) and then substitute in (7).
Remark 3. Everywhere, the fact that a value {■} is unknown and so is to be estimated is reminded by the notation {■} with the overscript °. When moving later to the solution, we change marking the estimated parameters to the commonly used {•}, instead of true {■}.
Directly from equations (1) and/or Fig. 1 (a), the following expression
Ts (s) = 2 (9)
s2 + a\s + a0
is obtained with the following three parameters
cq â 1/(LSCS),
ào â (RsGs + 1)/(LBCs) = ul,
a Rs G s ai = -j- + 77- = 2(un,
M = s-2 [ ao] = s-2
[ ài] = s
-i
(10)
Here are some intermediate values
Wn =
c =
a1
[Wn] = S , [<] = 1,
D = ao -ai/4 = ^ (1 - (2), [D] = s-2,
X
a\
c
2Vd yr—C2 '
1x1 = 1
(11)
introduced through the basic parameters (10) for further convenience, checking (2 < 1 for the literature obtainable secondary cable parameters [17, Table 2.1].
3. Problem statements. Given the specific case A with TF (9), good-quality estimates are required for parameters (10) of the numerator and denominator of this TF. In the most general case B, given an LTI 'black-box' as an nth order ordinary differential equation (nth order ODE), good-quality estimates are required for the numerator and denominator parameters of the corresponding TF. The solution is sought for the below A and B cases.
A. Illustrative example (9). The DSL environment is a multi-user transmission environment enabling a Central Office and the Customer Premises Equipment (CPE) to communicate in the downstream (from the CO to the different users) or in upstream (opposite) directions. The CO and the locally distributed CPEs are connected via twisted pair lines, each line belonging to one user. The twisted pairs are physically close to each other because they are bundled in a cable binder. Electromagnetic coupling between lines results in mutual interferences at all modems operating within the same cable [25]. These interferences known as crosstalk channels must be mitigated or, better, canceled.
Of two different kinds of crosstalk, namely near-end crosstalk (NEXT) and far-end crosstalk (FEXT), the latter represents the largest performance limiter in the xDSL system. A variety of suggestions have been made to reduce the impact of FEXT.
Most DSL and discrete multi-tone transmission (DMT) scenarios use the decomposition-based zero-forcing precoding (DBZF) to deal with FEXT. In DBZF, the transmit vector signal is pre-perturbed by the [N x N] precoder matrix P defined for each tone, where N is the number of users. For each tone, this number may be in thousands, matrix P is the inverse of the normalized (i.e. unit-diagonal) channel matrix H-¿m. Formally, Hnorm is the channel matrix H pre-multiplied by matrix HJ^, the latter being a diagonal matrix composed of inverse transfer coefficients of direct channels [16; 17, pp. 34-35]. Hence, for downstream transmission with efficient precoding, i.e. full crosstalk cancellation, it is necessary to know
the channel [N x N] matrix H, which consists of the DCTFs (on the diagonal) and crosstalk channel transfer functions, CCTFs (off the diagonal).
For a very short section (see Fig. 1), the generally accepted model DCTF is given by formula (9). A large number of such DCTFs cascade to form a DCTF of the entire line. It is shown in Fig. 2 with approximations dR = R-dl & Rs, dL = L-dl & Ls, dG = G-dl & Gs, and dC = C-dl & Cs, iout = kn + diin and Vout = Vin + dvvn, given Al & dl. Therefore, the resulting DCTF will be of a higher order, while remaining a proper fractional-rational function of s. As for the CCTF, solutions for its modeling include various approaches that have a solid physics basis but high computational complexity [17, p. 28]. Nevertheless, the adopted CCTF model does not escape the proper fraction form we move to now.
d R d L
^wvm-
■"in AC-
^out
d R d L
v
■d G
-^wvm-
Vmt
dC:
dG
Figure 2. Equivalent lumped RLCG-circuit of a 2-wire transmission line
B. General case. In the most general form, the transfer function of a channel to the jth output from the fth input is defined as follows:
T^s) =
cmsm + Cm-lSm 1 +-----h CIS + CQ
sn + an-isn-1 +-----+ äis + a®
(12)
where m + n + 1 < 2n + 1 parameters may be unknown. With (12), a DSL system is thought of as a MIMO—specifically, [N x N]—system, for which the crosstalk is modeled as an input rather than noise and the acronym MIMO stands for Multiple-Input Multiple-Output (Fig. 3).
Thus, the ¿th input Ui(s) causes a direct response Zi on the (j = f)th output and creates crosstalk contributions Zj on all other, (j = f)th outputs, plus an external noise Vj (s) in every j th channel:
N
Yj (s) = £ T,l(s)Ul(s) + V) (s), j = 1, 2,...,N.
(13)
i=1
The fact we are seeking to solve the inverse problem of recovering (12) from (13) dictates the only possible identification scenario (cf. Fig. 3): feed only one, namely ¿th input Ui(s) per single, namely ¿th identification session, into the MIMO system:
Y3 (s) = Tji(s)Ui (s) + V, (s), j = 1, 2,...,N. (14)
(As for the uppercase variables notations in (13) and (14), cf. Remark 2.) The time of each ¿th session (14) needs to be spent to determine the ¿th column T (s)^ = [fji(s)], j = 1,2,... ,N of matrix T (s) = [T^s)], and the whole scenario will require repeating N sessions: i = 1, 2,... ,N as in (14).
4. Problem solution framework. Focusing the research on the class of linear constant-coefficient ODEs to describe the wide range of LTI dynamical
Feedback Line i m Zl(t)
CE
SI
CI
Tu(s)
VI®
CPE1
Channel
i(t)
CPEi
from i to j
vN(t)
TNi(s)
zN{t)
Viv(i)
CPE,
N
Feedback Line i
Figure 3. The distributed MIMO channel estimation structure. Legend: CE - Channel Estimation; SCO - System Central Office; SI - System Information; CI - Channel Information; CPE -Customer Premises Equipment; N - the number of customers, j = 1, 2,... ,N
systems, we first state that the choice of the excitation signal Ui(t) (cf. Fig. 3) is extremely important for parameter system identification. Gaussian white-noise random excitations w(t) are very popular among practitioners because they seem to be simple to design. We also stick to this choice, assuming Ui(t) = Wi(t). However, using random-phase multisines for Ui(t) is also possible, given the design of the amplitude spectrum of the multisine is such that the equivalence between the random-phase multisine and the Gaussian random noise concerning the system behavior is guaranteed. Such signals are known as Riemann-Equivalent Excitation Signals, REESs [2, p. 44]. Using random input excitations makes the system output under study a stochastic process.
4.1. Cauchy form ODE system. Since there is no uniformity in the structure of matrices for the general Cauchy form and little else can be said about this form without additional knowledge of the particular dynamical system, we assume that the output, i. e. measurement data are generated by a completely observable physical system whose observability index is designated p. Hence, we focus the attention on the SOF among the known three standard system forms [26, pp. 2832]. The SOF provides a sort of unified approach to TFs of general form (12), not just (9). Besides, using SOF is beneficial to the below solution.
Given (9), using the notation mentioned in Remark 3 yields the following system of equations
d_ dt
X\ = 0 1 X\ + 0
X2 —o>o — ä\ X2 Co_
w(t)
(15)
F
with w(t) as a stationary input voltage v\(t) (cf. Fig. 1). Let w(t) be REES, that is Riemann-equivalent to the Gaussian white-noise excitation with the correlation function Rww(t) = Q5(t) in terms of Dirac's delta function 5(r) with some Q > 0, M = V2-s where Q is possibly given. Next, assume that the output voltage
v2(t) = x\(t) is measured with random error v(t) of Gaussian type with correlation function Rvv(t) = R5(r), R > 0, [R] = V2- s, to obtain the measurement data as
H
X\ X2
+ v(t).
(16)
SOF model (15)+(16) corresponds to conditions (i) and (ii) of the experiment above mentioned on page 549. Its characteristic polynomial q(s) = s2 + a\s + a0
O O O O O ° 1
has the discriminant —D < 0. Besides, Ts(s) = H®s(s)r with ® s(s) = (Is—F)-i, in matrix notation. The inverse Laplace transform of $s(s) yields the continuous-time state transition matrix
<hi(t) <h2(t) fai(t) fait)
with its entries
fai(t) = e-c^nt cos(tVD) + xsin(tVD)
$i2(t) = e
Vd
sin(tVD),
fai(t) = -ulfa2(t),
fait) = \cos(tVD) — xsin(tVD)
Given (12), it leads to the general SOF
X1 0 1 ■■ ■0 X1 bi
X 2 = 0 0 ■■ ■0 X2 + b2
% n— 1 0 0 ■■ ■ 1 bn—1
, - (1Q —a1 ■ ■ ■ — (ln—1 Xn . bn _
H
y(t) =[ 1 0 ■ ■ ■ 0 0 ] x(t) + v(t)
H
(17)
(18)
w(t),
instead of (15)+(16), where b1, b2,..., bn satisfy the following equation
0
0
Cm
Co
1 0 0
an-i 1 0
än-2 än-i 1
(12 äs ■ ■ ■
a1 a2 ■ ■ ■
^n-
0 bi
0 b2
0 k
0 bn—1
1 . bn _
fl
and Fji is the Frobenius companion matrix for the characteristic polynomial q(s) = sn + an-1sn-1 + ■ ■ ■ + a1s + a0 of (12). Acting as before (17) yields the
OOOO O ° 1
check relation Tji(s) = H$ ji(s)Tji, in which $ ji(s) = (Is — Fji)-1 serves to find ^ji(t) as the inverse Laplace transform of ^¿(s), quite similar to (17). Check: the Tji(s) thus found must coincide with (12).
Remark 4. What will be done in the next Subsec. 4.2. and Subsec. 4.3 based on the preceeding Subsec. 4.1 for the illustrative example given in (9) can be repeated similarly for the general case given by (12), furnishing the results with the subscript ji. We omit these details and ji subscripts due to the obviousness of the technique. We also omit subscript s as is done at the transition to (17).
4.2. Discrete-time model (DTM). Belonging of {■} to the discrete-time model is indicated below by the subscript d as in {-}d. Before making the change, it is necessary to reasonably choose the sampling interval T. Obviously, the sampling rate 1/T must be much higher than the natural frequency fn = 1/Tn = wu/(2k) of the system to be able to track the system behavior. This requirement means Tfn ^ 1. From the other side, (2, cf. (11), must remain less than one for the consideration to stand. As a result, requirement T ^ Tn in the transition to the discrete-time model means that parameter
d 4 e-^T
(19)
appearing in (18) at t = T must lie within the sufficiently wide boundaries of inequality e-2w ^ d < 1, that is be less than one, but possible insignificantly less. Given (9) and its continuous-time model (15)+(16), the DTM
x(tl+i) = $dx(ti) + wd(ti), = <P(T), 1 y(ti) = [1 0}x(ti)+ Vd(ti) 1 (20)
H )
yields by the standard method [26]. Here it is checked that (20) is observable with the observability index p = n = 2 and has physical dimentionalities [x\(ti)} = V, [x2(ti)] = V- s-1. The discrete white-noise wd(ti) in (20) is a zero-mean process
A fti+1 °
Wd(U) 4 - T)Tdß(r)
Jti
definded via the Brownian motion ß(t) related formally with w(t) in its differential dß(r) 4 w(t)dr. The covariance [2 x 2]-matrix of wd(ti) is [26]
fti+1
° A I O O O rp ° rr\
Qd 4 / - T)VQVT^T(tl+1 - T)dr. (21)
Jti
For the illustrative example (cf. 550), four entries of (21) are calculated directly
using (10) and (11); the result is in
qii =
4|| [1 - (d2 + 2X^DcPii(T)MT))]/
di2 = C-^^2(T )= Q21,
Q22 = [1 - (d2 - 2xVD^i2(T)fa2(T))]. 4 (,
(22)
Remark 5. Calculations like (22) are technically trivial, so they are omitted here. They may seem complicated if done manually. Manual work can be avoided by using Maple to obtain the result quickly, easily and accurately. Additionally, although dimensionality analysis does not guarantee the correctness of result, it can be an auxiliary tool as in this case: [qn] = V2, [q 12] = Vs-1, [q22] = Vs-2.
To finalize formulating DTM, it is worth going from Wd(U) to a dimensionless vector quantity (d(ti) for which wd(ti) = Ld(d(ti) with a matrix Ld such that
o o o ™
Qd = LdLT by the lower triangular Cholesky decomposition [27, p. 40]. From (22),
hi = /qü, [ I ii ]=V,
hi = qi2/lii, [12i] = V ■ s-i, (23)
I22 = \! q22 - hi2, [12i] =V ■ s-
i
being three non-zero real-valued entries of [2 x 2]-matrix ¿d. The model (20) now takes the final form
x(U+i) = $dx(U) + LdU(ti), = kT), ] y(ti) = [1 0}x(ti)+ Vd(U). 1 (24)
h )
As a result, the discrete white sequence £d(ti) in (24) has the unit covariance matrix, and the measurement discrete white sequence vd(ti) may have some unknown
covariance Rd > 0. Vector 6 = \61 = c® | d2 = a0 | O3 = a1 | d4 = Q | 05 = EdJ .
Vector § = \01 = do l 02 = a0 l §3 = = Q | 05 = Ed]T will be the estimator
for . The explicit dependence of the in (24) matrices on can be easily traced from the above formulas.
4.3. Standard Observable Discrete-time Model (SODM). Turning back to the general solution of the problem, case B, let us introduce
M = [H tKh $d)T|■■■|(н $r1)T]T,
the observability matrix for a linear n-dimensional one-way derivable DTM. It is invertible as the observability index p is supposed to equal n. Performing a nonsingular basis transform in the state space by relation x* = Mx, we obtain the Stan-
o A 00 o 1 o A°°1
dard Observable Discrete-time Model, SODM, with = M$dM-1, H* = HM-1,
and L* A MLd. Let us note that we use {*} as a superscript or subscript for any magnitude {■} belonging to the SODM, keeping in mind that the transfer function does not change when the basis changes nonsingularly. For this general case, notice Remark 4. For the specific case of (20), (24), we obtain the following SODM:
x* (ti+i) = (ti)+ LCd(ti), y(U) = H+x* (ti)+ vd(U),
0 d2
o A O O O -i
A M$dM-i =
2d cos(T VD)
H* A HM-i = [1 0],
hi
hi$ii + hi<Pi2
A
0
l22<fti2
u A MLd =
fiii A 0ii(T), $i2 A 0i2(T), 022 A 022(T). ,
(25)
Every SODM thus obtained has matrix in the form of the Frobenius companion matrix, and matrix H* with its first element equal to 1 and the rest to zeros.
4.4. OKF as the Target Model (OKF-TM). Using the above technique culminated in (25) and taking into account Remark 4, one obtains the following unique Optimal Kalman Filter-Target Model, OKF-TM, in the SODM basis for the general case arising from (12):
x*(U+ilU) = Q^falti), x*(tilti) = x*(tilti-i) + (tilti-i)
i),
(26)
together with
y(ti) = H+x* (tilt—) + u(tilti-i), V(tilti-i) a y(ti) - H*x*(tilU-i) is defined as Innovation Sequence, IS,
O 0 0 rp / o o o rp o . _ 1
K* = P-Hj(H+P-H? + Rd) i,
O OpO OOO^O .r-r-^ O O .r-r-^
P- = [P- - KAP--] + L^.
(27)
We aim for a parametric identification of the steady-state OKF-TM (26)+(27). In this filter, ^t|t-i is a white-noise Gaussian sequence (WGS), and the last two equations in (27) form a Discrete-time Algebraic Riccati Equation, DARE. Note the IS behaves like an WGS because d in (26)+(27) is assumed to be a true, albeit unknown, real parameter vector of some dimension q: 0 e R9.
Thus, algorithm (26)+(27) is a set of steady-state Kalman filter equations optimal for the true parameter 0. It is written under the unrealized assumption that d is known and that steady-state operation of this algorithm has been achieved by a theoretically assumed numerically stable DARE solution.
Remark 6. The preceding contains the correct characterization of v(ti\ti-i) provided that the mathematical model on which the filter (26)+(27) is based accurately represents the real behavior of the system.
Thus, we can imagine—and consider this representation fair reasonable and therefore bearable—that the observed output y(ti) provided in fact by the real system is as if were generated, very conventionally, by the target model, that is, by the optimal Kalman filter, as presended in the first line of (27).
4.5. Concurrent Candidate Models (CCM). Since d is unknown, one can only use its estimated value d, which is located in some 0 space. Where the object of interest exists with no specific constraints, the real-valued space 0 = R is formed by all possible values a[ ] of the estimator vector a. Here and below denotes the order number of value a in some scanning trajectory over the space 0: 0[j] G 0, j = 0,1, 2,..., Jtota\, i.e. over an imaginary set of Concurrent Candidate Models, CCM. The CCM set plays the role of Machine Learning Models if one prefers to use Machine Learning terminology. When implementing a numerical iterative filter optimization method capable of sequentially converging to the OKF-TM (26)+(27), j has the meaning of the method step number, since it is common to test suboptimal models sequentially, their total number Jtot3\, even if we admit, quite theoretically, the possibility of testing them in parallel (i.e. synchronously). It is important that in both variants, sequential or parallel, of the target model (26)+(27) identification, it is possible and even expedient to base the work on the same observational data y(ti) supplied (conditionally as said in Remark 6) by the target model (26)+(27), using processing and analyzing the responses to these data of the suboptimal models under test as candidates for the role of the target, that is, optimal, model.
Assuming that a has taken a particular a[ ] value in 0, imagine that instead of optimal Kalman filter (26)+(27), we have managed to implement a suboptimal steady-state Kalman filter we refer to as th Standard Observable Kalman Filter, the jth SOKF, or SOKF(0[j]), for short. The latter is the jth candidate model
g)(tm|ti) = $g*(tllti), 9j(UlU) = |ti-i) + K*.rjj(UIti-i), y(ti) = H*gj(tiI ti-i)+ rjj (til ti-i), r]j(tilti-i) 4 y(ti) - Hg)(tilti-i) is defined as the jth Residual Sequence, RSj,
K*3 = p-HT{H*p-HT + Rdj)-i P- = $[p- - ] $T + L*JLT3
(28)
with = $*(§[j]), Rdj = Rd(0[j]), and L*j = L*(9[j]). The model is intended to participate in testing to come as close as possible to the optimal filter (26)+(27), provided that the target model is also in the CCM set.
However, what does it mean: 'we have managed to implement (28) ?' In the real case scenario, this means that when trying to test candidate models sequentially, i.e. SOKF(0[j]) after SOKF(d[j — 1]), SOKF(d[j + 1]) after SOKF(%j), and so on, we must solve DARE, i.e. the last two equations in (28) at each such step. Doing this job iteratively for each SOKF(6l[j]), we introduce the local notation (i) for the iteration number, (i) = (0), (1),..., (Idare), where Idare denotes the
final iteration number, and compute as follows, labeling the computed quantities with j :
(0) P*j-l(lDARE + 1) (29)
and then
k*, (0 = P-^{H+P-uH + Rdj y1
J
P-(,+1) = $Kw - (-«] + (30)
, (i) = (0), (1),..., ( Idare).
The final value k^. (iDARE) should be used as K-j in the second equation of (28),
that is, := k^.(/DARE), and the final value P-^^E+1) as the starting point
p-+m = P-(lDARE+1) by the (29) type but now for the SOKF(%' + 1]) at
the (j + 1)th optimization step if any, over the CCM set. 'Real-case scenario' means that these Riccati iterations should be stopped at IdaRE when a reasonable convergence criterion is satisfied. It also means that by the time of the final iteration ( i) = (IDaRE) we assume that the so iterated filter (28) has reached the desired steady-state operation defined by equations (28).
Iterations (30) may and should be performed on an accelerated time scale in the form of known numerically robust algorithms, e.g. [28], for each value j, in other words, at each jth step of the numerical approximation to the optimum, that is to the target algorithm (26)+(27). This numerical optimization should be performed by a single AKF scanning sequentially the elements of the theoretically unbounded CCM set.
4.6. Predictors to form the AKF. We supplement the jth candidate model (28) with the predictors and make them operate as follows:
9*j(t i+h\ti) = $ g*(t i+h-i\ ti) Vj (ti+h\U) = H*g* (ti+h\ti)
h = 1, 2,... ,p (31)
where is the total observability index of the system.
Remark 7. In an n-dimensional system with m outputs, any % 'th output of m outputs can be assigned a partial observability index pi. The sum of partial observability indices is always equal to the dimensionality n of the system if only the system has the property of complete observability. The total observability index p of the system is defined as the greatest of the partial indices. The case p < n is possible if only m > 1. In the problem under consideration, m = 1, therefore p = n elsewhere in what follows. NevertheJess, we distinguish between the notations p and n, intending further work to extend the solution to the case where the number m of system outputs exceeds one. Only then p may occur less than n.
The fact that fjj(ti+h\ti) in (31) and beyond depends on d[j] can also be denoted by the subscript bearing in mind the equivalence of the two possible
notations: yj (ti+h\ti) = y§y](ti+h\ti), h = 1, 2,..., p. In (32) that follows for the
case of (20)+(25) when p = 2, the first line comes from (31) while the second equa-
tion in (32) comes from the first three lines of the target expressions (26)+(27):
9~*(t ml U), x*(ti+il ti) +
yj (t i+l \ t i) ' H
yj (ti+2 \ ti)_
y(t i+l) " H '
y(t i+2) HA.
1
v(t i+ll ti) v(t i+2l ti+l)
(32)
The composite (stackable) vectors opening expressions (32) in the specific case = 2 are to be redefined when turning to the general case. Their definitions follow using notation p for the total observability index:
(ti+i1 - [Vi(ti+l\ li) I I yi(ti+v\ li)]
y(t+) - [y(t + )I---Iy(t+)]T, ti+i — (tí+1, tí+2, . . . , ti+p) .
T
(33)
For the case of a single-output completely observable linear n-dimensional DTM, we have p = n. Thus, we obtain the advantages of changing to the SOF, viz.,
H () H HK
H* ^f = i,
00 , T
(h^- ^f = 1.
(34)
Equations (34) are true regardless of j and the non-trivial entries in the Frobenius matrices $ * as defined in (25) for (26) and as commented for (28).
4.7. The generalized residual (GR). What follows is the general case of using p = n in the key relations (34) as a result of computing these composite (stackable) vectors:
Vö\j]( lf 0 = 1■ ti (t + \t i),
y{ti$) = I-X*(t+lit i) + 1 0 ■ ■ ■ 0' HJ+K 1 ■■■ 0
0 0 r>— 1 0
HM 1K+
0 0 n—9 0
HM 2K
v(t i+ll ti) v(t Î+2I ti+l)
1 \y(ti+p\ti+p-
-
^(ti+p\ti+p—l) (ti+l\ti)
( )
(35)
Expressions (35), obtained at intervals equal to the system's observability index p, show that the discrepancy between the system outputs y and the corresponding predicted data yQ^t1^]U) both expressed in SODM terms, contain a valuable but explicitly unavailable mismatch between the states x*(ti+1lti) of the optimal filter (26)+(27), which is latently present in the discretely observed system output in response to the learning excitation w(t), and g*(ti+1lti) being
computed at the jth iteration of SOKF(%]) (28)+(31).
Remark 8. The mismatch of the optimal filter (26)+(27) compared to the suboptimal filter (28)+(31) is the difference between the object state estimates given by the optimal and suboptimal filters.
Naming the difference between the second and first lines in (33), or equally in (35), the Generalized Residual, GR, calculated to be the (n x 1) vector process
GR: efa U) â y(f+) - yê[j] (t^\ti) e R™
(36)
and introducing the notion of adaptive filter state estimation error, or better to say, the concept of Adaptive vs. Optimal Filter State Estimation Mismatch,
AOFSEM: e*è[j](tm\ti) â x*(ti+l\ti) - g*(t+\ti) e R
(37)
yields the key result:
Theorem 1. Let GR be calculated as (36) and AOFSEM, which does not have a computer-manipulable representation, defined as (37). Based on the fact that the Direct Performance Index
DPI â ^ (9) â || e*m(ti+i\ti)\\<\ e Rl,
(38)
or in other words, Expected Direct Cost Function, EDCF, is not explicitly available to optimize the suboptimal filter (28)+(31), we introduce the Indirect Performance Index
IPI4 Jlp+l0) 4 E{||4b1(ii+?|tOH2} G R1, (39)
or in formal words, the Expected Indirect Cost Function, EICF. Then minimizing the IPI (39) by any numerical optimization method in 9 = 9[j] G Q at each discrete time ti+p is equivalent to minimizing the DPI (38) in 9 = 9[j] G Q at time ti+1:
{minê JlZ(9)} ^ {wins J™0) = O}.
(40)
Proof. Given definitions (36) and (37), relations (35) show that
ti ) = ^ i+l^ i)+5
^(ti+plti+p— 1 ) 1 \ti)
( )
with ö
,(ti-\-p\t'
+ p\ii-{-p—1)
(ti+ l\ti)
(9) defined in (35) being, first, independent of the estimated
value 9 = 9[j] e Q and, second, uncorrelated with error ë*^(ti+l\ti) (37) since
the stackable vector
(ti + P \tr'
+p\ti+p— 1
/(t i|i.) formed by the white-noise IS in (26)+(27) is
separated by one sample time interval T from all preceding IS values that determine error ti+1\U) (37). It is this circumstance, together with the theoretical
fact that IS has the properties of a white-noise sequence, that entails statement
IPI = DPI + Const*
where Const* equals E
(i+p\ti+p—1) (ti+l\ti )
( )
a value
2
independent of 6 = d[j] £ 6, which definitively proves statement (40). It is also easy to verify that
^mN = U) (41)
by virtue of the predictors (31) and first line in (34). □
Remark 9. The data y(ti+p) defined in (33) and used in (36) does not depend on d[j], so only these measurement data can be collected and stored in computer memory before running the method's algorithm aimed at analyzing any SOKF(0[j]) in terms of its state (behavior) closeness to that of the target optimal filter and ability of further diminishing the mean square discrepancy between these states.
5. Method implementation challenges. An attempt to implement the given solution with its advantages poses several challenges concerning the organization of computation time, calculation sequence, and numerical stability. Let us briefly discuss these challenges.
5.1. Computation time organization. As noted above, it is possible and even expedient to calculate parameter estimates in the accelerated off-line mode, that is, after the accumulation of measurement data in a database. We come to the contents of the database having defined the function to be minimized as follows.
By shifting IPI (39) back by p time points, we determine the Expected Indirect Objective Function, EIOF,
f(§) 4 E {||el{j](tl_p+1IU-p)\\2} (42)
to be minimized in 9 £ R9. For practical work, we have to turn to the Averaged Indirect Objective Function, AIOF (43)
__m
ipi ^ m) 4 j—y, || 4w(ti-P+iitl-P)\\2 4 fM 0 = m (43)
k=0
considered as a real-valued function fM( 6 = d[j}) to be minimized in parameter § £ Rq instead of (42), the latter being the shifted back EICF (39).
Remark 10. The use of the upper index (k\ k = 0,M 4 0,1,...,M, from (43) on as the sample number is especially justified when averaging M + 1 sample paths of the process, if necessary. Otherwise, it is sufficient to let M = 0 and so release of using (k\ M may be as large as desired positive integer number for better-averaging. We relate 11, the time the computer starts up to process the data, to the real time ti by which the data is ready.
As seen from (43), each (k)th (ti-p+111■i-p)-path must result from a one-to-
one time-mapping—no more than in the algorithmic computations computer-time formal representation—of interval ti_p+1lto the real-time (fc)th segment
ti-P+i 4 (k+i)p+i, k-(k+i)p+2,..., k-(k+i)p+p} 4 i■_ ^.+1^+1 , 1 (44)
k = 0,1,...,M j
of a set of points. We imagine an entire record ( M + 1)p-length of all data y(t l-(M^p^ = { y(t i-(M+l)p+l), y(ti-(M+i)p+2^), ••• , y(t í-l), y(t 0} in the Measurement Data Base, MDB, as composed of ( M + 1) portions (45)
y(k)(fi-p+i) = y(f-*!+i)P+J , k = 0,1,...,M (45)
obtained in real time but referenced to the computer-time stackable p-vectors y(k)(tl-p+i) 4 [y(k) (t-+i) | y(k) (t,-p+2) | • • • | y(k)(tl)]T ,
k = 0,1,...,M.
Remark 11. The idea of notation (46) explained as regards the (k)th sample path l-p+i\ti-P) notation by relating the t\-p+i\ti-p interval—for the algo-
rithmic computations in computer-time—to the real-time segment (44) should be clear below when applied to other quantities as well. One can always see how realtime data, such as (45), relates to the same data, such as (46), when the latter is stored in the MDB for further computer processing. Additionally, the upper index (k indicates that the (k)th sample data is considered as being in the MDB.
(46)
Given (35), (36), and (37), let us write down all (k)-sampled time-shifted
values ^ • (fc) •
£e\i] №-p+i\= y(k(ti-p+i) — y(kj]{ti-p+i\t-) (47)
k =0,1,...,M
of the GR, (36), stored in the MDB to compute /m{&) and, respectively, all
e§\j\(i-'p+i\( \ti-p+i\h-p) — 9j( )(U-p+i\ti-p) (48)
k = 0, 1, . . . , M
AOFSEM, (37), to go from (43) to optimization algorithms.
According to Theorem 1, the (k)th sample value (47) of the random vector (36) varies, if we consider it in the Mean Square, MS, sense, by Const# remaining constant during the numerical scanning—does not matter sequentially or in parallel—of the parameter space 0, from the (k)th sample value (48) of the p-points delayed random discrepancy (37) between (A) state x*(k(ti-p+i\ti-p) of the target optimal filter (26)+(27) and (B) state g*(k\ti-p+i\U-p) of the suboptimal filter (28), that is from the error committed by (B) if used as the jth estimator of (A) based on the (k)th sample data (45)=(46). The theorem also proves equation (41), which we use hereafter as the p-point delay for any (k)th sampling path thus obtaining (49)
-p+i\U-p) = gkj(ti-p+i\ti-p). (49)
Remark 12. The delay by (k + 1)p points of discrete time from the current moment ti = 11 in (45), and in (47), (48) as well, if we relate (47), (48) to the
real time, is irrelevant because of the stationarity assumption of the system understudy and the MS sense used.
Note in passing that, by virtue of Remark 7, vectors in the (k)th sample defined by expressions (45), (47), and (48) have dimension n, due to equality p = n in the situation of scalar system output, m = 1.
Assuming that the 6 parameter has no explicit constraints, the unconstrained optimization is applicable. The most straightforward way of doing this gives rise to what is known as Newton's method (50) [29, pp. 44-79]:
(a) solve G(d[j]) A = -V§fM {, = 0[j]) for vector A;
(b) set §[j + 1] = §[j] + A(0 = §[j]) with A(0 = §[j]) 4 A
Hereafter, gradient operator is used as
.1
V (o) = & &
— (o)
(o)
ddq ( )
T
-t- (o)---
r=1,q
(50)
(51)
to be applied to each scalar element of vector or matrix (o); if (o) = fu{0 = 9[j}), we have (50) with the matrix of second partial variables G(9[j]) known as Hessian matrix and defined by V2-fT(d = d[j]) 4 fM(§ = §[j])T). When M = 0,
(50) is treated as a pure stochastic approximation of Newton's method.
The following gradient descent optimization
o[j + i]=o[j] -nv^m{o = e[j])
(52)
can be a reasonable alternative to (50) with a lower computational burden. It uses a small enough step size 7j G R+ to guarantee /m(0 = 0[j + < fu{0 = 9[j})
and thereby performs the transition to the next SOKF(0[j + 1]), as shown in Fig. 4 and mentioned in Subsec. 4.5 to test candidate models sequentially. In this figure, the arcs directed by arrows from points ti-p+l, ti-p+2, and ti to ti-p
on the horizontal line tell us, cf. (31), that all predictors g*(k\ti-p+hIti-p), their numbers h = 1,2,... ,p, involved in obtaining the estimates y^k)(ti-p+hIti-p) =
H*gf\t i-p+h It i -p) are conditioned, in the probabilistic sense, on the entire measurement history y 4 {y(ti-p-i),y(ti-p-{l-1)),..., y(U-p-i), y( ti-p)}, (theoretically, I ^ to) which is incorporated in the predictors and precedes their calculation immediately after time ti-p.
Thus, upon a closer look at what the theory requires, we realize that it requires the cost of time intervals (44) depicted in Fig. 4 by boxes as ti-k^+l(d[j]), then
tf-p+M + 1]), and so on, as related to SOKF(%]), then SOKF(§[j + 1]), and so on during the AKF step-by-step accelerated optimization.
5.2. Computation scheme and numerical robustness challenge. Let us
perform the necessary calculations for the above algorithm (52). Using (49) to obtain
V,
(J£e[j](t ti-p+1
I ti-P)f)
-V
(>;(fc)(i i-H-i
IU-P )]T)
(53)
T
ti—p p+1 tj—p+2' * * tj—1 tjJ
1]), feo.Miti-,,
Figure 4. A timing diagram for minimizing the average (43) of the indirect performance index (39) by a gradient sequential method (52) using the AIOF as defined in (43)
and further (43) as (o) in (51) yields (54), the Objective Function Gradient, OFG:
vêfM{ ê = ê[j])
M
m+Ï^ v* № l-'+>iu-> >Hx
-2
k=0
M
£ô\k](^ ^-p+i1 ^i-p)
M + 1
k=0
£ <Kk\ ('i-'P+i1 ^i-P)
(54)
Controlling whether the OFG norm reaches a neighborhood of zero, that is checking it for values 'greater than/equal to' or 'less than' (<) a small threshold 5 is a convenient criterion for continuing or terminating the procedure:
\vêfM (ô = ê[j])\
continue
< S.
stop
(55)
The identification procedure will continue repeating operations numbered (Fig. 5) by blocks ®, ®, ®, ®, ©, and © with the same data stored in block ® after updating the parameter estimate in block ®:
0[j + 1] := 0[j] or will stop with capturing the result in block ®:
$FINAL := 6[j] .
(56)
(57)
The most important thing about these repeatable operations (cf. Fig. 5) is that the calculations in block © - the AKF, and in block © - the AKF Sensitivity Model being an algorithm to compute values of partial derivatives of state vector estimates and covariance matrix elements relative to the ^-vector parameter estimates d[j], must be numerically stable and robust with respect to ill-conditioned models. In this regard, it should be noted that practical projects using Kalman filtering, KF, since the very first one [30], have opened a wide field for research and development on giving KF algorithms, including Riccati equations, required numerical stability and robustness properties. The fundamental ideas of Bierman [27] about matrix factorization served as a powerful impetus.
For the method of Active System Identification developed in this article and earlier works by the author, the greatest contribution from Russian scientists was made by Julia Tsyganova and Maria Kulikova in their dissertations [31,32] and numerous publications, partly co-authored [33, 34] with the author of this synthesis paper. More references on the pioneering titles can be found in a recent survey [35]. In there, one can find discussions and current developments for the efficient and robust computation of derivatives on the parameters of discrete filter equations, including a set of vector-valued filter sensitivity equations and a set of matrix-valued Riccati-type sensitivity equations arising in implementing the (steepest) gradient descent method (52); the necessary experimental proofs are there also available. In promising software projects, it is strongly recommended to create modern implementations of © and ® blocks (in Fig. 5) based on orthog-onalized array methods for parametric identification of discrete linear stochastic systems [31].
m
v{k)(tlP+1)
9j{ti-p+i\ti-p)
-E> *{k) e\j] ^ ® -p+i 'i—p)
Vol M{6 = = %1)
-v,
•§ ([9f\ti.p+1\ti.p)]T)
Ö[7 + l]=ö[?]-7iVÄ/M(fl = ö[j]) continue
m -=o\j+i]
Öfinal := 6[j]
Figure 5. Generalized Parametric Identification Scheme by the Active Principle. Legend: ® = (45) = (46); © = (31)^(29)^(30)^(31)^(49); ® = (47); ® = (53); © = (54); © = (55);
® = (52); ® = (56); ® = (57)
5.3. Sequence of work for a software project. If there is an object of interest in the application domain, an appropriate mathematical model must first be written for parametric identification of the TF, considered to be an adequate description based on the laws of Physics. An example of such preprocessing operations is given in Sect. 2. It is highly expedient to perform further actions on the application of this method not manually, but in the automated mode on the Maple software, according to the guidelines in Remark 5. The calculation procedure recommended for identifying the TF of an object according to the method outlined above is shown graphically in Fig. 6. This diagram can be useful when creating a specialized software tool to solve similar problems if such a project arises. Application-specific (AS) calculations for the problem taken in this article as an illustrative example are dictated by the following intermediate quantities: © parameters wn, (, D, and % in (11); © matrix ^(t) in (17) with its entries in (18); © parameter d in (19); © matrix Qd in (21) with its entries in (22); © matrix Ld with its three non-zero entries in (23); © matrix in (24); and © matrices and L* in (25).
These quantities are all continuous and differentiable functions of elements Oj, of vector 9 introduced after (24) for this application problem. Following the form of these functions, we apply not the true parameter value , but its estimated value 6 to have continuous and differentiable, over d, functions. These properties make it possible to compute the gradient V^J^t?= @[j]) in algorithm (52)
tuning the parameter ,.
Figure 6. Flow-diagram of works for a software project to implement the Active Principle of parametric system identification in the class of linear, time-invariant, completely observable models. AS = application-specific and SM = standard matrix calculations. The paper section numbers and equation numbers (within parentheses) are shown
6. Conclusions. Returning to the problem issues posed at the beginning, we believe the goals have been achieved.
The incompatibility of the theoretical (direct) criterion of optimality of the system model and practical methods of optimization has shown itself to be the most difficult or even impossible obstacle to overcome under conditions of uncertainty directly. In this article, incongruity between theoretic perception of the system optimality and practical on-computer optimization is overcome through the construction of an indirect proximity criterion of the adaptive and optimal system models to each other, which can become a practical tool for parametric system identification. This is done here for a class of linear time-invariant models characterized by a transfer function, relying on Kalman filter theory. As proven here, it is necessary and sufficient to modify a physically specified structure into the discrete-time standard observable model in both the observed data and the adaptive Kalman filter to implement this idea.
The advantage of this modification is twofold. First, the restrictions on the permissible size of a'priori parametric uncertainty are removed because the measurement channel parameters are transfered formally into the modified state equation. Second, and most importantly, we replace the direct, but unattainable for identification, objective function with the indirect objective function, which is equivalent to the original direct function, but—to our satisfaction—attainable and suitable for conventional optimization methods. The work of this preliminary modification is difficult to perform manually. Fortunately, it is greatly facilitated and accurately performed using well-known symbolic computation tools (Maple).
We prefer to denote the above-used concept of 'equivalence' of the two objective functions by the new term 'equimodality,' which simply means the coincidence of their minimizers in the argument space. This coincidence is important because ensures that minimization of the constructed practical indirect cost function does lead to the unbiased state estimates along with the unbiased parameter estimates, as it should be in the optimal filter.
The inclusion of a new illustrative example from modern digital communication technology increases the visibility of the approach and thereby encourages its extension to broad real-world applications.
Theoretical questions should not overshadow the practical side of the case. Therefore, this paper also includes practical critical issues: the organization of calculations in the computer time-scale, the structural algorithmic construction of the identification procedure, and the planning of the corresponding design work.
As a limitation of the work, it could be pointed out that underlying its results are the assumptions that the system to be modeled is linear and time-invariant. There is nothing to argue against this, except for the famous aphorism by George E.P. Box [36, p. 2]: "Models, of course, are never true, but fortunately it is only necessary that they be useful."1 Accordingly, if the system output measurement shows the presence of nonlinear distortions, it may mean that the suitability of the proposed method to identify a nonlinear model with a dominant linear kernel defined by the Best Linear Approximation concept, BLA, based on [2] propositions, should be considered and recommended for extended study.
1https://www.tandfonline.com/doi/epdf/10.1080/01621459.1979.10481600? needAccess=true&role=button
Competing interests. The author declares no competing interests.
Authorship contribution and responsibility. The author has approved the final
version of the manuscript.
Funding. The research has not had any funding.
Acknowledgments. In the preparation of this paper, the author has benefited from the
discussions with Dr. Alexandru Murgu (University of Cape Town) who provided extended
comments and gave valuable suggestions on the manuscript.
References
1. Pillonetto G., Ljung L. Full Bayesian identification of linear dynamic systems using stable kernels, Proc. Natl. Acad. Sci. USA, 2023, vol. 120, no. 18, e2218197120. DOI: https://doi. org/10.1073/pnas.2218197120.
2. Schoukens J., Vaes M., Pintelon R. Linear system identification in a nonlinear setting: Nonparametric analysis of the nonlinear distortions and their impact on the best linear approximation, IEEE Contr. Systems Magaz., 2016, vol.36, no. 3, pp. 38-69. DOI: https:// doi.org/10.1109/MCS.2016.2535918.
3. Ljung L. Convergence analysis of parametric identification methods, IEEE Trans. Auto. Contr., 1978, vol.23, no. 5, pp. 770-783. DOI: https://doi.org/10.1109/TAC.1978. 1101840.
4. Ljung L. System identification, In: Signal Analysis and Prediction, Applied and Numerical Harmonic Analysis; eds. A. Prochazka, J. Uhlir, P. W. J. Rayner, N. G. Kingsbury. Boston, MA, Birkhauser, 1998, pp. 163-173. DOI: https://doi.org/10.1007/978-1-4612-1768-8_ 11.
5. Ljung L. System Identification: Theory for the User. Upper Saddle River, N.J., Prentice Hall, 1999, xxii+609 pp.
6. Gevers M. System identification without Lennart Ljung: what would have been different?, In: Forever Ljung in System Identification; eds. T. Glad, G. Hendeby. Lund, Sweden, Stu-dentlitteratur AB, 2006, pp. 61-85.
7. Schoukens J., Pintelon R., Rolain Y. Time domain identification, frequency domain identification. Equivalencies! Differences?, In: Proc. 2004 American Control Conf., vol. 1 (30 June 2004 - 02 July 2004). Boston, MA, USA, 2004, pp. 661-666. DOI: https://doi.org/ 10.23919/ACC.2004.1383679.
8. Oomen T., van Herpen R., Quist S., et al. Connecting system identification and robust control for next-generation motion control and a wafer stage, IEEE Trans. Contr. Syst. Tech-nol., 2014, vol. 22, no. 1, pp. 102-118. DOI: https://doi.org/10.1109/TCST.2013.2245668.
9. Dedene N., Pintelon R., Lataire P. Estimation of a global synchronous machine model using a multiple-input multiple-ouput estimator, IEEE Trans. Energy Conver., 2003, vol. 18, no. 1, pp. 11-16. DOI: https://doi.org/10.1109/TEC.2002.805198.
10. Wei Y., Mantooth A. LLC resonant converter—frequency domain analysis or time domain analysis, In: 2020 IEEE 9th Int. Power Electronics and Motion Control Conf. (IPEMC2020-ECCE Asia; 29 November-02 December, 2020). Nanjing, China, 2020, pp. 552-557. DOI: https://doi.org/10.1109/IPEMC-ECCEAsia48364.2020.9367734.
11. Rivera D. E., Lee H., Mittelmann H. D., Braun M. W. Constrained mulisine output signals for plant-friendly identification of chemical process systems, J. Process Contr., 2009, vol. 19, no. 4, pp. 623-635. DOI: https://doi.org/10.1016/j.jprocont.2008.08.006.
12. Peeters B., Ventura C. E. Comparative study of modal analysis techniques for bridge dynamic characteristics, Mech. Syst. Signal Process., 2003, vol.17, no. 5, pp. 965-988. DOI: https://doi.org/10.1006/mssp.2002.1568.
13. Westwick D. T., Kearney R. E. Identification of Nonlinear Physiological Systems. Pis-cataway, N.J., Wiley-IEEE Press, 2003, xii+261 pp. DOI: https://doi.org/10.1002/ 0471722960.
14. Semushin I., Tsyganova J., Kulikova M., et al. Identification of human body daily temperature dynamics via minimum state prediction error method, In: Proc. 2016 Europ. Control Conf. (ECC2016; June 29-July 1, 2016). Aalborg, Denmark, 2016, pp. 2429-2434. EDN: YVHZMP. DOI: https://doi.org/10.1109/ECC.2016.7810654.
15. Semoushin I. V. Identifying parameters of linear stochastic differential equations from incomplete noisy measurements, In: Recent Development in Theories and Numerics (International Conference on Inverse Problems; January 9-12, 2002); eds. H. Yiu-Chung, Y. Masahiro, C. Jin, L. June-Yub. River Edge, NJ, World Scientific032.62073, 2003, pp. 281290. DOI: https://doi.org/10.1142/9789812704924_0026.
16. Semushin I. V., Tsyganova J. V. Reducing computational complexity for DBZF precoding in xDSL downlinks, J. Phys.: Conf. Ser., 2018, vol. 1096, 012159. EDN: OHQTID. DOI: https:// doi.org/10.1088/1742-6596/1096/1/012159.
17. Düngen M. Crosstalk Mitigation Techniques for Digital Subscriber Line Systems, Dissertation (Dr.-Ing.). Hamburg, Technische Universität Hamburg, 2016, 160 pp. DOI: https:// doi.org/10.15480/882.1293.
18. Begovic A., Skaljo N., Behlilovic N. A simple model of copper local loop for use in current DSL application, In: Proc. 2015 23rd Telecommunications Forum Telfor (TELFOR; November 24-26, 2015). Belgrade, Serbia, 2015, pp. 114-117. DOI: https://doi.org/ 10.1109/TELFOR.2015.7377427.
19. Begovic A., Skaljo N., Behlilovic N. An example of modeling local loop transfer function in DSL environment, In: Proc. ELMAR-2014 (September 10-12, 2014). Zadar, Croatia, 2014, pp. 1-6. DOI: https://doi.org/10.1109/ELMAR.2014.6923362.
20. Rodrigues R. M., Sales C., Klautau A., et al. Transfer function estimation of telephone lines from input impedance measurements, IEEE Trans. Instrum. Meas., 2012, vol.61, no. 1, pp. 43-54. DOI: https://doi.org/10.1109/TIM.2011.2157431.
21. Foubert W., Neus C., Boets P., van Biesen L. Modeling the series impedance of a quad cable for common mode DSL applications, In: Proc. 2008 IEEE Instrumentation and Measurement Technology Conf. (IMTC 2008; May 12-15, 2008). Victoria, BC, Canada, 2008, pp. 250-253. DOI: https://doi.org/10.1109/IMTC.2008.4547040.
22. Neus C., Boets P., van Biesen L. Transfer function estimation of digital subscriber lines with single ended line testing, In: Proc. 2007 IEEE Instrumentation and Measurement Technology Conf. (IMTC 2007; May 01-03, 2007). Warsaw, Poland, 2007, pp. 1-5. DOI: https://doi. org/10.1109/IMTC.2007.378980.
23. Bostoen T., Boets P., Zekri M., et al. Estimation of the transfer function of a subscriber loop by means of a one-port scattering parameter measurement at the central office, IEEE J. Sel. Areas Commun., 2002, vol.20, no. 5, pp. 936-948. DOI: https://doi.org/10.1109/ JSAC.2002.1007376.
24. Semushin I. V. Adaptation in stoshastic dynamic systems—Survey and new results II, Int. J. Commun. Netw. Syst. Sci., 2011, vol.4, no. 4, pp. 266-285. DOI: https://doi.org/10. 4236/ijcns.2011.44032.
25. Stolle R. Electromagnetic coupling of twisted pair cables, IEEE J. Sel. Areas Commun., 2002, vol.20, no.5, pp. 883-892. DOI: https://doi.org/10.1109/JSAC.2002.1007371.
26. Maybeck P. S. Stochastic Models, Estination, and Control. Vol. 1, Mathematics in Science and Engineering, vol.141. New York, Academic Press, 1979, xix+423 pp. DOI: https:// doi.org/10.1016/s0076-5392(08)62169-4.
27. Bierman G. J. Factorization Methods for Discrete Sequential Estimation, Mathematics in Science and Engineering, vol. 128. New York, Academic Press, 1977, xvi+241 pp. DOI: https://doi.org/10.1016/s0076-5392(08)x6052-7.
28. Sima V. Algorithms for Linear-Quadratic Optimization, Pure and Applied Mathematics, Marcel Dekker. New York, Marcel Dekker, 1996, vii+366 pp. DOI: https://doi.org/10. 1201/9781003067450.
29. Fletcher R. Practical Methods of Optimization. Chichester, John Wiley & Sons, 2000, xvii+436 pp. DOI: https://doi.org/10.1002/9781118723203.
30. Potter J. E., Stern R. G. Statistical filtering of space navigaton measurements, In: Proc. 1963 AIAA Guigance and Control Conf. (AIAA 1963; August 12-14, 1963). Cambridge, MA, 1963, pp. 1-5. DOI: https://doi.org/10.2514/6.1963-333.
31. Tsyganova J. V. Orthogonalized Array Methods for Parametric Identification of Discrete Linear Stochastic Systems, Diss. Dr. Sci. (Phys. & Math.). Ulyanovsk, Ulyanovsk State Univ., 2017, 400 pp. (In Russian). EDN: IMPPAX.
32. Kulikova M. V. Methods of Calculation of the Logarithmic Likelihood Function and its Gradient in Kalman Filtering Algorithms, Diss. Cand. Sci. (Phys. & Math.). Ulyanovsk, Ulyanovsk State Univ., 2005, 131 pp. (In Russian). EDN: NNNPKP.
33. Semushin I. V., Tsyganova J. V. Numerically implementing the API based solution to learn a plant dynamics for stochastic control concurrently, In: Proc. 2020 European Control Conference (ECC 2020; May 12-15, 2020). St. Petersburg, 2020, pp. 1105-1110. EDN: FWJUQS. DOI: https://doi.org/10.23919/ECC51009.2020.9143870.
34. Kulikova M. V., Semoushin I. V. On the evaluation of log likelihood gradient for Gaussian signals, Int. J. Appl. Math. Stat., 2005, vol.3, no. 5, pp. 1-14.
35. Tsyganova J. V., Kulikova M. V. On modern array algorithms for optimal discrete filtering, Vestnik YuUrGU. Ser. Mat. Model. Progr., 2018, vol.11, no. 4, pp. 5-30 (In Russian). EDN: YOTRJJ. DOI: https://doi.org/10.14529/mmp180401.
36. Box G. E. P. Some problems of statistics and everyday life, J. Am. Stat. Assoc., 1979, vol.74, no. 365, pp. 1-4. DOI: https://doi.org/10.1080/01621459.1979.10481600.
Вестн. Сам. гос. техн. ун-та. Сер. Физ.-мат. науки. 2023. Т. 27, № 3. С. 544-572 ISSN: 2310-7081 (online), 1991-8615 (print) d https://doi.org/10.14498/vsgtu2037
EDN: XUQAGA
УДК 681.5.015
Идентификация передаточной функции посредством минимизации рассогласования оценок состояния адаптивного и оптимального фильтров
И. В. Семушин
Ульяновский государственный университет,
Россия, 432017, Ульяновск, ул. Л. Толстого, 42.
Аннотация
Статья посвящена дальнейшему развитию активного принципа параметрической идентификации системы в классе линейных, инвариантных во времени, полностью наблюдаемых моделей. В качестве целевой модели идентификации выбран оптимальный фильтр Калмана (ОФК), который не более чем концептуально присутствует в дискретно наблюдаемом отклике системы на обучающее возбуждение типа белого шума. Путем модификации физически заданной структуры в стандартную наблюдаемую модель как в наблюдаемом отклике, так и в адаптивном фильтре Калмана (АФК), строится так называемый обобщенный остаток (ОО), равный рассогласованию между оценками состояния адаптивного и оптимального фильтров плюс независимая от АФК шумовая составляющая. В силу этой модификации средний квадрат ОО становится новым критерием близости модели для этих фильтров. Минимизация этого критерия с помощью обычных практических методов оптимизации дает точно такой же результат (АФК = ОФК), как и минимизация теоретического критерия, который, к сожалению, недостижим для любых методов численной оптимизации АФК. В статье представлена подробная пошаговая процедура, объясняющая вышеуказанное решение в терминах параметризованной передаточной функции. Для наглядности и стимулирования применения подхода в реальном мире в статье используется модель передаточной функции линии витой пары в типичной системе xDSL. Обсуждаются проблемы реализации теоретических положений метода. Вопрос о распространении предложенного подхода на проблемы идентификации линейных моделей для нелинейных систем обозначен в направлениях дальнейших исследований.
Математическое моделирование, численные методы и комплексы программ Научная статья
© Коллектив авторов, 2023 © СамГТУ, 2023 (составление, дизайн, макет)
3 ©® Контент публикуется на условиях лицензии Creative Commons Attribution 4.0 International (https://creativecommons.org/licenses/by/4.0/deed.ru) Образец для цитирования
Semushin I. V. Transfer function identification by minimizing the adaptive vs. optimal filter state estimates mismatch, Vestn. Samar. Gos. Tekhn. Univ., Ser. Fiz.-Mat. Nauki [J. Samara State Tech. Univ., Ser. Phys. Math. Sci.], 2023, vol. 27, no. 3, pp. 544-572. EDN: XUQAGA. DOI: 10.14498/vsgtu2037.
Сведения об авторе
Иннокентий Васильевич Семушин А https://orcid.org/0000-0002-3687-1110 доктор технических наук, профессор; член IEEE; профессор; каф. информационных систем; e-mail: [email protected]
Ключевые слова: ЦП-модель, полная наблюдаемость, фильтр Кал-мана, адаптивный фильтр, косвенный показатель качества, проблемы реализации.
Получение: 23 июня 2023 г. / Исправление: 5 августа 2023 г. / Принятие: 21 августа 2023 г. / Публикация онлайн: 28 сентября 2023 г.
Конкурирующие интересы. Автор заявляет об отсутствии конкурирующих интересов.
Авторский вклад и ответственность. Окончательная версия рукописи одобрена автором.
Финансирование. Исследование выполнялось без финансирования.
Благодарности. При подготовке этой статьи автору были полезны дискуссии с доктором А. Мургу (Университет Кейптауна), который дал расширенные комментарии и ценные предложения по рукописи.
