IDENTIFICATION OF AN UNMANNED AERIAL VEHICLE PILOT-OPERATOR DYNAMICS MODEL UNDER STOCHASTIC CONDITIONS

The article is devoted to investigation of the new rules for identifying the pilot and his remnants dynamics models in the conditions of work as a part of a multidimensional closed-loop control system for an unmanned aerial vehicle. The article subject is the analysis of existing methods for identifying the dynamics model of a multidimensional object (pilot) as a closed-loop system part. An additive mixture of a vector random useful signals and random noises is allowed at the object output. The research purpose is to substantiate the object (pilot) frequency characteristics and noises’ (remnants) spectral densities matrices determining rules. All signals are assumed to be centered stationary random processes. The research method is based on the use of the control systems statistical dynamics rules. The application of the rules is illustrated by an example in the Matlab programming environment.


Introduction
Problem formulation.Unmanned aerial vehicles (UAVs) have become more and more widespread in recent years.As a rule, their flight is controlled within the framework of a cyber-physical system (Fig. 1), which central element is the pilot operator [1].
Operators control actions are formed on the basis of comparing visual information about the flight plan r formed during preparation for the flight and the actual flight parameters y received from the sensors installed on board the UAV.Control actions are transmitted to the aircraft rudders drives using a digital data transmission system [2].Thus, the quality of solving the object flight control under consideration problems is largely determined by the pilot-operator characteristics, as an element of automation.Studies of the complex multidimensional objects operators psychophysiological and motor characteristics [3] have shown the effects of self-learning and self-organization in a person who is a part of a close loop control system as a regulator.A characteristic feature of the UAV flight control systems operating conditions is the random nature of atmospheric disturbances, measurement noise and transmission of navigation parameters acting at the inputs of this system.As a result, the structure of the pilot-operator model, as an element of automation (Fig. 1), can be represented [3] as consisting of two parts.
The first element, the operator () n Ws, characteriz-© Osadchy S., Ivliiev A., Kolisnichenko S., Tymoshenko G., 2023 es dynamics of the interrelation between the change in visual information u received by the pilot and UAV controls position change x , made by him.The second element, the operator () p Ws characterizes controls' additional movements carried out by the pilot during the control ρ.These motions are called preemption or remnants.In a general case, additional motions are considered to be stationary random processes [4], which can be formed from a unit intensity white noise vector δ.
Analysis of recent research and publications.Depending on the visual information composition and the controls number, these operators can be either transfer functions or transfer function matrices.As the studies of various authors [5][6][7][8][9][10][11] have shown, presented above operators structure and parameters depend on many factors, which can be conditionally divided into two parts.The first part is determined by control system design, flight control task and execution conditions.The second part of the set is determined by pilot state and his training degree.
Pilot dynamics operators models knowledge allows not only to assess his work quality, psychophysiological state and training degree, but also to form the initial data for synthesizing the UAV flight automated remote control system, which has the highest possible quality under certain conditions [5].In addition to this, the presence of these models will allow one to successfully solve the problem of building an adaptive humanmachine cyber-physical avionics system, posed in the article [7].
Given the operators n W and p W structure and parameters depend on many factors, the most appropriate method for obtaining pilot's models is a structural identification [5].An analysis of such identification various methods [6][7][8][9] made it possible to identify a number of limitations that do not allow using them to determine the pilot's and remnants' dynamics models as multidimensional dynamic objects operating in stochastic conditions.The main of these limitations are: the requirement that the disturbances acting on the UAV in flight, the errors in measuring (setting) the vectors r , y (Fig. 1) represent random "white noise" type processes [7][8][9] ; representation of the pilot-operator as a dynamic object with one input and output without taking into account the remnants [10; 11]; requirement of vectors u and r mutual influ- ence absence [7][8][9]; sensors represent ideal measuring instrument [10; 11]; presentation of the identification process as a black box problem solution procedure [5; 6].
An article objective.This article is devoted to the development of a pilots' and remnants' dynamics models structural identification method, which would be free from the above limitations.

Statement of basic materials 1. Purpose and task of research
The predicted spatiotemporal flight path of any aircraft can be represented as a program signals vector r (Fig. 1), which carries information about the flight plan.Information about the flight parameters actual values is contained in the vector y of the results of UAV output signal vector 1 x components measuring using the on board sensors system.The interrelation between these vectors equation has the form [4]: where 1 W is the sensors transfer functions matrix, which dimension is determined by the vector 1 x components number; - is a centered vector of stationary random measurement noise.
The UAV dynamics model is known, for example, from the monograph [2], and can be represented as a system of ordinary differential equations with constant coefficients: where 1 P , 1 M are appropriate size polynomial matrices from the differentiation operator; -1 u is the UAV rudder control signals vector; - is the vector of disturbances acting on the in- struments in flight.
Taking into account the accepted designations (1), (2) allows to bring the UAV flight control system architecture shown in Fig. 1 to a block diagram (Fig. 2).Let us assume that all components of the program signals vector r , the results of measurements vector y and changes in the position of the UAV controls vector x are measured and recorded, and the disturbances and noise of measurements represent vector stationary random processes.Since components of the program signals, disturbances and measurement noise vectors have different sources, the hypothesis of the relationship between them absence is accepted.Then the task of developing a structural identification method is formulated as follows.Using the known polynomial matrices 1 P , 1 M , the sensors transfer functions matrix 1 W , the transposed spectral densities matrices of disturbing influences / S  and of measurement noises / S  , as well as the vectors r , u , x , y components records, find the transfer func- tions matrices n W and p W .

Basic research
The analysis of the block diagram (Fig. 2) made it possible to obtain the following equations system that connects vectors of different signals: where .
The vector 1 x exclusion from the equations system (3) and taking into account expression (4) lead to the result: in which matrices 0 K and 10 K are obtained as a result of the Multifractional Decomposition described by the formula: Bringing into consideration two vectors: the recorded output signals of the system (Fig. 2)  and the input signals  of the form: Makes it possible to represent the system of equations (5) in a vector-matrix form: ; , where n E is an nn  identity matrix; n O is a zero matrix of the same dimension; n is a number of the vector x components.The equation ( 8) implies that: where 1 G is the matrix of closed-loop UAV flight control system transfer functions from the vector  to the vector  , which is equal to: Using the Frobenius formula [12] for inverting the block matrix 0 G and performing the multiplication prescribed by expression (11), it is obtained the rule for determining the transfer functions matrix 1 G : where 0 10 0 1 .
As it is known from the multidimensional control systems' statistical dynamics [4], if two vector stationary random processes are interrelated by an equation of the form (10), then there is a correlation between the transposed matrices of the spectral and cross spectral densities of the vectors  and  .superscript / denotes transposition [12].

This relationship is characterized by an expression
According to the Wiener-Khinchin theorem (relation) for vector random processes, the transposed matrix of the vector  spectral densities is defined as Its elements are transposed matrices of mutual spectral densities between vectors indicated in subscripts.
If matrices ( 14), (15) elements are known by the problem statement, then the matrix 1 G is easily determined from equation (13) as Comparison of the elements of this matrix and the matrix (12) makes it easy to find an identification problem solution.
For example, the pilot transfer function matrix n W is given by the equation ( ) The second method is based on the application of the Wiener-Khinchin theorem to the vector (10) with taking into account the matrix (12)

AS A F S F W AS A F S F W AS A F S F W AS A F S F W F K M W S A F S F F K M W S A F S F W AS
S is a fractional rational matrix which elements depend on the unknown matrix of transfer functions p W so that To solve the equation ( 19) with respect to the matrix p W , the matrix (23) was transformed to the form 2 0 3 0* , S T S T = (25) Thus, the additional communication equation is represented as ( ) where the transposed matrix / S  is found from the records of the vector  components, obtained experimentally, in accordance with expression (18), and the subscript * denotes the Hermitian conjugation [8].Equations ( 18), ( 27) and (28) comparison allows one to find six independent equations for finding the pilot remnant spectral densities transposed matrix / S  .
In accordance with the Wiener-Khinchin theorem, taking into account (Fig. 2), this matrix is equal to To search for this matrix, the first block from 3 S of dimension nn  was used in the equation ( 28).As a result, the following rule for determining the matrix (29) was found:

W S W M K S K M S S S S S S S S S S S S
It is advisable [4] to find the remnant transfer functions matrix p W as a result of the Wiener factorizing the matrix on the right hand of the equation (30), since  is a stationary white noise of unit intensity according to the problem statement.
To implement this method, before factorization, it is necessary to approximate the matrix elements on the right side of the equation (30) on the class of complex argument's sj = fractional-rational functions.
Thus, the identification problem has been solved and the rules ( 16), (30) for calculating the pilot transfer functions matrix and the remnant (leading) transfer functions matrix from experimental data have been obtained.

Results
In order to analyze the effect of using the proposed new method for identifying the pilot and remnant dynamics models the UAV pitch angle stabilization system operation simulating problem was posed and solved.The stabilization system corresponds to the block diagram in the Fig. 2 with the following symbols and data.
1.The program signal r is equal to the given value of the pitch angle, has a zero mean, and its dynamics is characterized by the spectral density 2. At the pilot's input, there is a control signal u , which is equal to the deviation of the measured pitch angle y from its given value.
It is necessary: to develop a simulation model of the stabilization system; to get records of signals r , u , x , y ; to identify estimates of the pilot transfer function n W and remnant spectral density / S  ; to evaluate the degree of correspondence between estimates and true values.
To solve the problem, a simulation Simulink model of the stabilization system was developed in the Matlab 2022a environment (Fig. 3).
As can be seen, the model includes four unit intensity white noise generators WNG1 -WNG4.WNG1 together with the shaping filter Grr serve for create a program signal r corresponding to the spectral density (31).The generator WNG2 together with the transfer function p W form a remnant ρ.WNG3 together with the shaping filter G11 serve to create a UAV disturbance, the spectral density of which is determined by expression (34).The remaining WNG4 generates the measurement noise.Additive signals shown at the outputs of the generators are designed to compensate generators systematic zero drift.The pilot model that implements the transfer function (32) consists of two blocks (Fig. 3) p W and a pure delay block connected in series.A four-beam oscilloscope O1 is used as a recorder.As a result of stabilization system (Fig. 3) operation modeling, realizations of four random processes were obtained for 700 seconds, fragments of which are shown in Fig. 4-7.
These processes are stored in an out object created in the workspace of the Matlab environment according to the O1 oscilloscope data    Suu=frd(Suu0,2*pi*F); where dz is the identifier of a variable that stores the program signal; ks is the identifier of a variable that stores the vector ς; Srr, Suu, Srx, Sry, Sru are identifiers of variables that store frd-objects [14] corresponding to spectral densities rr S , uu S , rx S , r S  , ru S .
Since the system consists of elements with one input and one output, the spectral and mutual spectral densities of the signals coincide with the transposed matrices in the notation from expressions (16) and (30).
The substitution of initial data into equation ( 16) was carried out as a result of the executing the command: Wn=Srx/Sru; where Wn is the identifier that saves the operator transfer function, obtained as a result of identification in the frd format [14].Comparative analysis of the Bode diagrams (Fig. 8) of the identified transfer function and the origin transfer function OW1 constructed from the data according to the expression (32) shows their good agreement.
The identified remnants spectral density's correctness verification was carried out as a result of the two commands groups' execution.The first command made it possible to find the first term on the right side of the equation (30), namely: The second commands group allows one to find the following spectral densities sum: S22=minreal(S22); where P1, M1, K0, K10, WFR, Spsipsi, Sfifi, S22 are identifiers of variables that store the original data (31-34) in zpk format [15].As it is followed from the equation (30), a balance must be maintained between the spectral densities S11 and S22 (Fig. 9).

Conclusion
An obtained results analysis allows one to draw a number of conclusions and formulate restrictions on the scope of justified identification rules: justified identification rules make it possible to obtain correct results under conditions when the disturbances acting on the UAV in flight, the measurement errors (setting) of the vectors r , y (Fig. 1) represent stationary random processes with arbitrary dynamics; the pilot-operator can be represented as a multidimensional dynamic object; it is allowed to have mutual influence between vectors of program signals at the system input and signals at the operator's inputs; sensors are non-ideal measuring devices with known metrological characteristics; the identification process is reduced to solving the gray box problem; to use the proposed identification method, it is necessary to provide an active identification experiment with synchronous recording a certain signals set; random signals in the control loop should be close to stationary.
Directions for further use of the research results: planning experiments to identify models of the pilot dynamics; obtaining pilot and remnants dynamics models; identification of the multidimensional object dynamics as a tracking system part.

Fig. 1 .
Fig. 1.Unmanned aerial vehicles flight control system Source: developed by the authors.

Fig. 2 .
Fig. 2. Structural diagram of the UAV flight control system Source: developed by the authors.
S  is the transposed spectral densities matrix of the vector  ; / S  is the transposed cross spectral densities between vectors  and  matrix;

Fig. 6 .
Fig. 6.An elevator position change graph Source: developed by the authors.

Fig. 8 .Fig. 9 .
Fig. 8. Bode diagram of a pilot transfer function Source: developed by the authors.
 is the transposed matrix of white noise spectral densities.The cross spectral densities transposed matrix  , found on the basis of the Wiener-Khinchin theorem, has the form / S / S