Lecture
In the study of automatic control systems (ACS) for aircraft engines (AE), a class of ACS with measurement constraints is of particular interest. Analysis of known identification methods showed that most of them are oriented toward objects with full measurability. In practice, however, it is not always possible to fully measure the state variables [3, 7, 22]. Thus, ACS-AE objects are characterized by limited measurability, which requires a more complex identification algorithm and a separate study of the measured object. For hydromechanical units, which serve as the actuating mechanisms of ACS-AE, the main dynamic coordinates are the displacements of the moving elements. Such a selection is reasonable, since displacements are the output variables of the typical dynamic functional elements of hydraulic systems. However, not all displacements are measurable variables. The design features of elements such as a spool valve or an isodrome do not allow the displacement of their stroke to be measured, since the motion of the pistons of these elements occurs inside a closed volume.
53
For this class of objects, the identification problem reduces to the simultaneous estimation of parameters and unmeasured coordinates. Analysis of known identification methods showed that there are several approaches to solving this problem. Of greatest interest is joint parameter and state estimation, which involves augmenting the state vector with the unmeasured coordinates.
Let the object be described by the system of equations
dX (t)
= F
(
X
t
,U
t
, A ,
dt
( )
( )
)
(1.76)
(
(
)
( )
)
Y
(
)
= G
X
,U
t
t
t
, C ,
where X (t) is the state vector; U (t)
is the vector of input variables;
Y (t) is the vector of output variables;
A is the coefficient matrix; C is the
output matrix.
If the system parameters are constant in time, then the vectorized matrix of unknown parameters can be represented as follows:
dA(t)
= 0.
(1.77)
dt
Then the system of equations (1.76) can be written as
dX (t) = F(X (t),U (t), A),
dt
(1.78)
dA
= 0.
dt
It is easy to see that even if the dynamics of the object are linear, this system of equations cannot be transformed into the form
dX (t)
dt
=
dA
dt
X
D . (1.79)
A
Consequently, the problem of joint parameter and state estimation is nonlinear with respect to the parameter and state vector. In general, only iterative methods exist for solving nonlinear equations, and it is on these methods that the algorithms for simultaneous parameter and state estimation are based:
54
they are:
•quasilinearization;
•the invariant embedding method;
•Kalman–Bucy filters.
1. Quasilinearization.
Quasilinearization is an iterative method based on the Newton–Raphson method.
The joint parameter and state estimation problem is formulated as follows.
It is required to solve the system of nonlinear differential equations:
dX
F(X (t),U (t),A)
dt
,
(1.80)
=
dA
0
dt
subject to the following boundary conditions:
Y (t) = CX (t).
(1.81)
Thus, a multipoint boundary value problem has been formulated. The quasilinearization method proposes the following algorithm for solving this problem:
dX (i +1) = F(X (i),A(i)) + F
(i)(X
(i +1) − X (i)) +
dt
1
+F2 (i)(A(i +1) − A(i)),
(1.82)
dA(i +1)
= 0
dt
subject to the condition
CX (i +1) =Y (i +1),
(1.83)
where F1(x),F2 (x) are Jacobians,
F (i) =
dF
,
1
dX
X (i),A(i)
55
F (i) =
dF
.
2
dA
X (i),A(i)
When applying iterative methods, including the quasilinearization method, it is necessary to assess their convergence. Studies have shown that convergence of the method is guaranteed if the function F(X ,A) is convex, and the elements of the Jacobians F1(x), F2 (x), located
above and below the main diagonal, must be positive. A convergence assessment for a 3rd-order element — and such elements are
among the most common in ACS-AE — showed the following: the presence of negative feedback loops leads to negative elements in both the matrix F(X ,A) and the Jacobians
F1(x),F2 (x).
Thus, convergence of this method is not guaranteed, and therefore the quasilinearization method cannot be applied for joint parameter and state estimation of ACS-AE.
2. The invariant embedding method.
The invariant embedding method is based on transforming a multipoint boundary value problem into a single-point one.
Unlike the previous methods, this method is less dependent on initial values. Convergence of the identification algorithm to the actual parameter values is ensured over a fairly wide range of initial approximations. The drawbacks of the invariant embedding method include the significant time required for transforming the model, which makes it impossible to apply for real-time identification of control objects, as well as the dependence of the algorithm's convergence on the choice of auxiliary matrices.
3. The Kalman–Bucy filter.
The Kalman filter algorithm makes it possible to construct, in real time, an optimal estimate of the system state X (t), based
on measurements Y (t) that contain errors. The measurement vector X (t) is treated as the multidimensional output signal of the system, noisy in this case, while the state vector Y (t) is the unknown
multidimensional signal to be determined. The optimality condition for the constructed estimate is the minimum of the mean-square deviation of the estimated value from the true one.
56
In general, the Kalman filtering problem is formulated for a nonlinear nonstationary system under certain assumptions. A modification of the method — the Kalman–Bucy filter — is applied to a linear stationary object represented in state space under white noise action [6, 8]:
dX (t)
= AX (t) + BU (t) + N (t) ,
(1.84)
dt
where X (t) is a random Markov n-dimensional process specified by the necessary a priori information: M{X} = X0 is the expected
value of X (t); cov(X (t)) = M{(X (t) − M{X (t)})(X (t) − M{X (t)})T}= P0 is the co-
variance matrix of X (t); U (t) is the measured vector input action, which can be either a deterministic or a random quantity; N(t) is a k-dimensional vector of random disturbances,
assumed to be white-noise-type processes: M{N(t)} = 0,
cov(N(t)) =
= M{N(t)N(t)T} =Vδ(t −τ), where δ(t) is the delta function; V
is the symmetric
non-negative definite intensity matrix of the white noise N(t).
The process X (t) is observed by means of a measuring instrument, and the vector of
measured coordinates by the relation
Y (t) = CX (t) +ε(t),
(1.85)
where ε(t) is a p-dimensional vector of measurement noise, assumed to be a random
process
in the form
of white
noise: M{ε(t)} = 0, cov(ε(t)) =
= M ε(t)ε(t)T
}
= Rδ(t −τ),
where R
is the symmetric non-negative
{
definite intensity matrix of the white noise ε(t).
The processes N(t) and ε(t), as well as X (t) and N(t), X (t) and ε(t), are assumed to be uncorrelated: M{N(t)ε(t)T}= 0, M{X (t)N (t)T}= 0,
M{X (t)ε(t)T}= 0.
It is required to construct a linear dynamic system that ensures
obtaining an optimal estimate of the vector, given
X (t) X (t)
that the estimation error is defined as:
57
E(t) = X (t) − X (t).
(1.86)
The optimality criterion is the condition that its quadratic norm be minimal:
J (t) = E(t)E(t)T → min.
(1.87)
The Kalman–Bucy filtering problem corresponds to the following block diagram, shown in Fig. 1.19.
Fig. 1.19. Block diagram of the forming filter (FF) and the measuring instrument under the action of random disturbances and in the presence of noise
Constructing the Kalman filter uses the idea of an n-dimensional observer (a full-order observer), in which the mathematical model of the system is taken as the state identifier.
The Kalman filter performs a recursive estimation procedure based on observations of the input and output signals of the object, where a corrective feedback loop on the system output Y (t) is introduced into the identifier algorithm to reduce the variance of the estimates.
Based on the requirement of obtaining an unbiased estimate of
X (t),
the filter equation has the form [6, 8]
d X (t)
= AX (t) + BU (t) + L(Y (t) −C X (t)),
(1.88)
dt
where L is the filter gain matrix, providing an estimate of the state that is optimal in the sense of minimum variance, and is determined by the expression
L = PCT R−1,
(1.89)
58
where P is the covariance matrix of estimation errors, which in the case of stationary processes is determined by solving the algebraic matrix Riccati equation:
AP + PAT − PCR−1CP +CVCT = 0.
(1.90)
In general, constructing an optimal observer amounts to solving an optimal stochastic control problem under conditions of incomplete information about the state variable vector. Finding the filter gain matrix can be accomplished using the analytical regulator design method [6, 23, 24].
Filtering in the case of sampled measurements of the function of time
k=1, 2,K has the following algorithm:
1. From the a priori values of the signal characteristics, the a priori value of the signal covariance matrix, based on k observations, is found:
Q(k +1) = AP(k)AT +V.
(1.91)
2. The a posteriori value of the signal covariance matrix, based on k +1 observations, is calculated:
P(k +1) = Q(k +1) −Q(k +1)CT [CQ(k +1)CT + R]−1CQ(k +1). (1.92)
3. The gain matrix L(k +1) is determined, setting the weight of the corrections to the initial conditions based on the covariance matrices of the state estimate
L(k +1) = Q(k +1)CT [CQ(k +1)CT + R]−1 = P(k +1)CT R−1. (1.93)
The quantities P(k +1), Q(k +1), L(k +1) are fully determined by the a priori information. Calculations continue until the filter reaches stationarity, the condition for which is the equality P(k +1) − P(k) ≤ ξ or, correspondingly, L(k +1) − L(k) ≤ ξ, where ξ is the
specified accuracy.
Once the matrix L(k) is determined, the operating algorithm of the Kalman–Bucy filter reduces to sequentially processing the incoming input and output data, whereby the current signal estimate is obtained by correcting the previously made estimate using the information arriving at the filter input during observation at each step: this requires additional a priori information, namely
59
the probabilistic characteristics of the disturbances and noise. These requirements are excessive for estimating the parameters and state of a deterministic object.
Thus, the studies described showed that none of the above parameter and state estimation methods is universal. Each of the methods considered has its own drawbacks that limit its applicability.
Therefore, when identifying ACS-AE models with measurement constraints, it is necessary to develop an applied identification method that is simple to implement, meets the requirements of convergence and speed of solution, and is oriented toward a specific class of models.
To solve the problem at hand, the following identification algorithm is proposed [22, 25].
The object is described by the system of equations
V (k +1) = ФV (k),
(1.94)
U (k)
where V (k) is the generalized state vector, V (k) = X (k); U (k) is the
Y (k)
vector of input actions; X (k) is the vector of unmeasured coordinates; Y (k) is the vector of measured coordinates; Ф is the transition matrix.
Let us consider a limited class of models characteristic of the hydromechanical units of ACS-AE. According to the classification given in [25], objects of this class can be represented by the following types of elements:
•a first-order lag element;
•an ideal integrating element;
•a real integrating element;
•a real integrating element with a summed input;
•a series connection of a 2nd-order lag element and an integrating element.
The first two types of elements are fully measurable and therefore can be estimated by any identification method for fully measurable objects. The real integrating element and the real integrating element with a summed input have one unmeasured internal variable. The object comprising a series
connection of a 2nd-order lag element and an integrating element has two unmeasured variables.
The transition matrices of these objects have the following form:
• real integrating element
1
0
0
a11
;
(1.95)
Ф = b1
0
0
a21
0
• real integrating element with a summed input:
1
0
0
a11
;
(1.96)
Ф = b1
0
a21
b2
0
• series connection of a 2nd-order lag element and an integrating element:
1
0
0
0
a11
0
Ф =
b1
0
(1.97)
.
0
a21
a22
0
0
0
a32
1
For 2nd-order elements, the generalized transition matrix has the
form:
1
0
0
a11
0
Ф = b1
.
a21
a22
b2
Given the structure of the transition matrix, the dynamics equation (1.94) can be written as follows:
u(k +1) = u(k),
x(k +1) = b1u(k) +a11x(k),
(1.98)
y(k +1) = b2u(k) + a21x(k) + a22 y(k).
After eliminating x(k) — the unmeasured coordinate — from expression (1.98), it takes the form
61
1
( y(k
+2) −(a11 +a22 ) y(k +1) +
b1a21 +b2 (1−a11)
+a11a22 y(k)) = u(k)
or
c1 y(k +2) +c2 y(k +1) +c3 y(k) = u(k),
where
1
c1
=
,
b1a21
+b2 (1−a11)
c2
=
−(a11 +a22 )
,
b1a21
+b2 (1−a11)
c3
=
a11a22
.
b1a21
+b2 (1−a11)
(1.99)
(1.100)
(1.101)
(1.102)
(1.103)
After observing the object, resulting in the formation of the measurement matrices of the input U (k) and output Y (k) variables,
equation (1.100) is solved with respect to the coefficients c1, c2, c3.
From the known coefficients c1, c2, c3, the coefficients of the transition matrix are determined:
• real integrating element:
b1a21 = 1 ,
c1
a
+a
22
= −
c2
,
(1.104)
11
c1
a a
=
c3
;
11
22
c1
• real integrating element with a summed input:
b a
21
+b (1
−a ) =
1
,
1
2
11
c1
a
+a
22
= −
c2
,
(1.105)
11
c1
a a
=
c3
.
11
22
c1
62
The coefficient b2 represents the response of the system to the input
action over the first discretization interval, and is therefore not difficult to compute.
Thus, identification determines the co-
efficients of the transition matrix: a11, a22, b1a11 or b2, a11, a22, b1a11. This means that the proposed identification algorithm allows
determination of the time constants a
=
1
,
a
=
1
and the overall
11
T1
22
T2
gain coefficient K = K K
2
=
b1a21
.
Σ
1
a11
For 3rd-order elements, in which the two unmeasured internal variables are x1,x2, given the structure of the transition matrix, the dynamics equation (1.94) can be written as follows:
u(k +1) = u(k),
x1(k +1) = b1u(k) + a11x1(k),
(1.106)
x2 (k +1) = a21u(k) + a22 x2 (k),
y(k +1) = a32 x2 (k) + y(k).
After eliminating the unmeasured co-
ordinates from it, expression (1.106) takes the form
c1 y(k +3) +c2 y(k + 2) +c3 y(k +1) +c4 y(k) = u(k),
(1.107)
where
1
c
=
,
(1.108)
1
b1a21a32
c
= −
a11 +a22 +1
,
(1.109)
2
b1a21a32
c
=
a11a22 +a11 +a22
,
(1.110)
3
b1a21a32
c
= −
a11a22
.
(1.111)
4
b1a21a32
63
After observing the object, resulting in the formation of the measurement matrices of the input U (k) and output Y (k) variables, equation (1.107) is solved with respect to the coefficients c1, c2, c3, c4.
From the known coefficients c1, c2, c3, c4, the following are determined:
• the coefficients of the transition matrix:
b a
a
=
1
,
1
21
32
c1
a
+a
22
+1 = −
c2
,
11
c1
(1.112)
c3
a
a
+a +a
22
=
,
11
22
11
c1
a
a
= −
c4
;
11
22
c1
• real integrating element with a summed input:
b1a21 +b2 (1−a11) = 1 ,
c1
a
+a
22
= −
c2
,
(1.113)
11
c1
a
a
=
c3
.
11
22
c1
The system of equations (1.113) is solved with respect to
a11, a22,
b1a21, a32, i.e., just as for 2nd-order elements, the time con-
stants a =
1
, a
22
=
1
and the overall gain coefficient are determined
11
T1
T2
KΣ = K1K2 K3 = b1a21a32 . a11a22
Thus, the proposed approach to solving the identification problem with unmeasured coordinates makes it possible to obtain a simple and practically realizable algorithm for identifying objects of this class.
Comments