A methodology for developing the digital twin mathematical model for a cyber-physical system in the air heating process form using an electric heater is proposed. At the passive identification expense, the technique gives possibility of synthesis and adaptation of the digital twin discrete mathematical model by the air heating measured data in the operating electric heater. The influence of analytical model uncertain parameters of the electric heater on the calculation's accuracy is considered. The quality criterion for assessing the mathematical model adequacy to a physical dynamic process is proposed. The algorithm of passive identification the mathematical model uncertain parameters is developed, in which the deviations of calculated model variables from measured values in the process of air heating are minimized. A numerical study of the considered method has been carried out. It is shown that the passive identification of the uncertain model parameters belongs to the oneextreme optimization problem. The considered simulation examples prove that the proposed methodology for the digital twin development using MatLAB software is an effective and convenient tool for the rational creation of a digital twin for cyber-physical systems.
Modem IT technology has advanced a great deal and is making it possible to increase the efficiency of industrial production to a great extent. Experts around the world predict the final arrival of the digital age in the near future. The prospect of digitalization looks attractive, progressive and efficient. The basic technological unit of digital production is the cyber-physical system (CPS) [1]. There are different interpretations of CPS concept because these systems are simultaneously located at the different spheres intersection of human activity. The main characteristic of these systems is the interaction between physical and computational processes. The interaction mechanisms are provided by a network structure known as the Industrial Internet of Things (IIoT).
Fig. 1 (a) shows a classic production control hierarchy. Here PLC are program logic controllers, SCADA is Supervisory Control and Data Acquisition (dispatching system and management), MES is Manufacturing Execution System (production management system), ERP is Enterprise Resource Planning (enterprise resource management system). However, this structure does not meet the modem production requirements. Today, production management requires Industry 4.0 competences, which are relatively new and largely IT-related. Thus, the classic hierarchical pyramid of production management is being transformed into a system with an IIoT (see Fig. 1, b).
Changing the concept of cyber production management dictates the need to develop new methods for the analysis and synthesis of such systems. The main characteristics of CPS are its heterogeneity, uncertainty of different nature, multidisciplinary nature, and provision of functioning guaranteed strategy throughout system life cycle [2]. CPS should be able to analyze multidimensional data, taking into account the latent factors of production. Based on this data, it can autonomously solve optimization problems and make the right decisions. A key unit in the CPS maintenance and management process is the digital twins.
The digital twin concept as a basic prerequisite for managing a physical object throughout its life cycle (design, implementation, operation maintenance and disposal) has been proposed by Michael Grieves [3]. The digital twin concept is part of the fourth industrial revolution and is designed to help enterprises detect physical problems faster, predict their results more accurately, and produce better products. The use of a digital twin during the creation phase allows the entire life cycle of the intended object to be modelled and simulated. The digital twin use makes it possible to simulate the entire life cycle of the intended object [4]. Most often, digital twins are created to simulate objects directly related to industrial production or being an important element of technical systems [5]. The first practical definition of a digital twin was given by NASA in an attempt to improve the modelling of physical spacecraft models in 2010 [6].
According to the Industrial Internet Consortium (IIC) [7], a digital twin is a formal digital representation of some asset, process or system that captures attributes and behaviors of that entity suitable for communication, storage, interpretation or processing within a certain context. Researchers point out that digital twin’s design is based on the use of simulation methods that provide the most realistic representation of a physical object in the virtual world [8]. Mathematical representation of digital twins can be obtained applying methods: statistical modeling, using machine learning or analytical modeling, which uses a mathematical description the physical processes laws.
Statistical analysis methods can be divided into three groups: regression analysis models, classification models and anomaly detection models [9]. The machine learning method choice depends on the size, quality and data nature, as well as the problems type to be solved. In [10], the authors point out that the analytical model’s determinism in engineering is a valuable property that should be exploited in CPS. In [11], CPS development using deterministic models, which have proven to be extremely useful, is discussed. Deterministic mathematical models of CPS are based on differential equations and include synchronous numerical logic and single threaded imperative programs. However, CPS combine these models in such a way that determinism is not preserved.
The development of digital twins is taking place within the concept of IIoT, which encourages developers in the standardization direction. One of the fundamental works on digital twins’ standardization is the Industrial Internet Reference Architecture (IIRA) reference model proposed by IIC. The document describes guidelines for the development of systems, solutions and applications using the Internet of Things in industry and infrastructure solutions. This architecture is abstract, and provides general and consistent definitions for different stakeholders, system decomposition, design patterns, and terms glossary.
The IIRA model for the IIoT identifies at least four stakeholder perspectives: business, use, function, implementation. Each perspective focuses on the functional components of the IIoT, their structure and relationships, the interfaces and interactions between them, and the relationship and interaction of the system with external elements of the environment to support the use and CPS operation. These features are of particular interest to CPS architects, developers and integrators. Based on the IIRA model, digital twin information includes (but is not limited to) a combination of the following categories [7]: physical model and data; analytical model and data; temporal variable archives; transactional data; master data; visual models and calculations. Based on the IIRA model, digital twins are mainly created for commercial use using simulation and optimization techniques, databases, etc. to implement rational real-time production management and support strategic and operational decision making. The concept of a digital twin has a multifaceted architecture and correspondingly complex mathematical support for its implementation. The development of a digital twin mathematical model, which after machine learning or identification gives the ability to follow and predict the physical process behavior, is the subject of many scientific discussions. In this paper the approach to develop a mathematical model of a digital twin, taking into account its fractionalrational uncertainty is proposed.
The aim of the publication is to develop digital twin mathematical model of CPS for manufacturing plants, taking into account the parametric uncertainties of the physical process mathematical description. The application of the mathematical model synthesis methodology for the electric heater digital twin is given.
Let's consider a mathematical model of the air heating process on an electric heater TE d dt E E k0 N E k1 A , TA d dt A A k2 E k3 A0 k4 GA , Td d dtd A d A k5 d A0 k6 GA ; here
TE cE M E ,
KE
KE 0 F0 , k0 1 KE , k1 1 ; TA cAKMA A ,
K A cA GA 0 F0 , k2
(
KA k3 1 k2 , k4 cA A0 A ; Td VA , k5 1, k6 dA0 dA .
KA GA GA
Mathematical model (
However, the heat transfer coefficient 0 depends on many factors. Its numerical value can vary
several times depending on: the heating element temperature E ; temperature A , humidity d A and
air expense GA ; heat exchange surface design features and other factors. This parameter is the subject
of thermal engineering research. Numerous papers have been published on methods for calculating
the heat transfer coefficient, which are based on experimental studies and similarity theory [13]. The
parameters TE , TA , k0 , k2 k4 of model (
The linear properties of the mathematical model must also be taken into account. The linearization
of model (
It is convenient to use Laplace transforms to find a solution to (
1 d A k5 d A0 k6 GA ; Td p 1 here a2
TE TA 1 k1 k2 , a1 TE TA ; b0 1 k1 k2
k3 1 k1 k2 , b1 k T
3 E 1 k1 k2 , b2
k4 1 k1 k2 , b3 k T
4 E 1 k1 k2 b4 k0 k2 ; p is the Laplace operator.
1 k1 k2
The mathematical model (
Y( p) W ( p)X( p) ; (3) A , W
W11 here Y d A 0 0 W22
W13 W23
W014 , X A0 d A0 GA N E T ;
W11 a2 pb12 p a1bp0 1 , W13 a2 pb23 p a1bp2 1 , W14 a2 p2 b4a1 p 1 , W22 Td kp5 1 , W23 Td kp6 1 .
Applying the inverse Laplace transform, it is easy to find an analytical solution (
Let's analyze variable parameters influence of mathematical model (3) on electric heater HE 36/2 dynamic properties, simulation of which is considered in [12]. Let's assume that the electric heater basic static mode parameters θA0, θA, dA0, dA are determined using CFS sensors readings. Suppose, for model (3), due to reduction of heated air flow rate from GA 0.43 kg s to GA 0.2 kg s (the new function parameters will change
5.6 p 1 from
W p 2.4 p2 6.7 p 1 0 nominal operating mode of the heater, all other conditions being equal) the heat exchange coefficient is changed from 0 161 mW2оС to 0 100 mW2оС . At the same time, the model (3) matrix transfer to W p 8.2 p2 11.3 p 1 8.2 p2 11.3 p 1 8.2 p2 11.3 p 1 . (5) 0 0.911p 1 0 0
The numerical values of the transfer matrices (4) and (5) differ significantly. Fig. 3 shows the transient simulation results. In this figure and further on the symbols are used: the output vector Y (4) c) d)
HE 36/2 through the action channels X Y (the graphs shows output variables received increments as the result of changing 0 ); а) A0 Y , A0 1 оС; b) d A0 Y , d A0 1 g/kg; c) GA0 Y , GA0 0.1 kg/s; d) N E Y , N E 1kW for the reference mathematical model with the transfer matrix (4) where the initial values of parameters 0 and GA were used; the output vector Y for the model with the transfer matrix (5) where the new values of parameters 0 and G A were used. The output vectors Y and Y contain the variables: A and A are heated air temperature; dA and d A are heated air moisture content. In transients’ simulation for HE 36/2 electric heater a control action X was used by successive stepwise changes in the input variables. The input vector X contains the variables: A0 , d A0 are temperature and humidity of electric heater input air; GA is air flow through the electric heater; NE is electric power supplied to the electric heater elements.
The heat transfer coefficient 0 affects the numerical values of the matrix transfer function of the model (3). The dynamic characteristics of the air heating process are significantly influenced by this parameter, as demonstrated in the simulation graphs (see Fig. 3). For these reasons, specialists in the field of heat exchange processes simulation should determine this parameter quite accurately, since its numerical value significantly affects the investigated process properties.
In paper [14] mathematical models uncertainties classification in the control of organizational and technological objects is considered. In general terms, the electric heater mathematical model can be represented by the input-output relation
Y p, q W p, q X p , (6) here W p, q A p A p, q B p B p, q is transfer matrix, A p, q and B p, q are numerator and denominator of perturbation, p is Laplace operator, q is vector of uncertain parameters. According to the proposed classification, the electric heater mathematical model has fractional-rational uncertainty, which is proposed to be identified by numerical methods.
The electric heater mathematical model (3) is developed on the basis of theoretical laws for a specific physical process. Therefore, the analytical model requires specification of some parameters. Formally for model (3) there are six changing parameters , GA, θA0, θA, dA0, dA (see Fig. 2). In the 0 process of identification, it is necessary to determine two parameters 0 , GA, parameters of the basic static mode θA0, θA, dA0, dA can be estimated, using CFS sensors readings or with sufficient accuracy these parameters can be determined from the data sheet of the electric heater.
Model (3) parameters identification is proposed to be performed in the passive experiment operation by measured signals values of electric heater input vectors X and output vectors Y. As the identification criterion we use the least squares criterion, which is defined as the measured values error square of the physical process Y and the estimates of the identified model (3) output vector Y at the same input action X
I M t0 t f Y Y T Q Y Y dt min , (7)
t0 here t0 is the initial trend time, tf is the trend duration, Q is the unit square matrix, T is the matrix transpose operator, M is the mathematical expectation operator, which takes into account industrial perturbations when measuring variables by sensors. The general structural diagram of parameters identification of the electric heater mathematical model is shown in Fig. 4.
It is recommended to use numerical zero-order optimization methods to identify the mathematical model parameters [15], since the search function is not set analytically and is calculated directly during the search algorithm implementation. Given the substantial computational resources to implement numerical methods, passive identification can be implemented at the enterprise management middle level within a decision support system (DSS) [16].
Electric heater
Model (3) search min, , (t)
The algorithm for identifying the parameters of a mathematical model consists of steps. 1. The CFS sensors monitor the input vector X(t) and the output state Y(t) of the heater in real time and the data is processed in the DSS, thus forming the time base of the physical process trends. 2. The output state Y t of the identifiable model (3) with the parameters initial values 0 , GA is estimated from the time trends of the input action X(t) on the electric heater. 3. Based on the values of the output Y(t) vectors and the identifiable Y t states, the identification quality criterion (7) is determined. 4. According to criterion (7), the parameters 0 , GA of the identified model (3) are optimized using numerical optimization methods. 5. If the minimum of criterion (7) is found, move on to the development of the digital twin mathematical model, else move on to step 2 and continue the process of identifying model parameters.
During the identification process, the optimization method (7) can have its own specific features that need to be taken into account. For qualitative identification, the time interval tf must be several times longer than the duration of transients in the heater. Also, it is necessary to set the parameter variation limits 0 and GA, based on the physical feasibility of the duct heater model.
The proposed identification algorithm was implemented using MatLAB software package. The reference model with transfer matrix (4) was used to form time trends X(t) and Y(t) of the functioning electric heater. The random signal with amplitude ±0.2 was added to the reference output variables in order to simulate production disturbances in DSS measurement channels. In the identified model the initial values of parameters 0 , GA differed significantly from the reference values 0 , GA. The MatLAB function fminsearch(...) was used for parametric identification of 0 and GA with the NelderMeade simplex optimization method. The main results of the numerical study are presented in Fig. 5. and Fig. 6.
Fig. 5 shows the case where the initial values of the identification parameters are as follows: 0 161, GA 0.43; GA 0.65. With a step change in the input action vector 0 120, Xt 1, 1, 0.2, 1T , produces the transients for the output of the reference Y t and the identifiable Y t model, which are shown in Fig. 5 (a).
The difference between the output values of vectors Y t and Y t is significant. In Fig. 5 (b) shows the surface isolines of criterion (7) and its minimization trajectory. From Fig. 5 (b) shows that the quality criterion has one extremum in the area of identification parameters. Identification parameters 0 143.7, GA 0.428 were determined using the above algorithm. The numerical value 0 was found with low accuracy, while GA was determined quite accurately. Identification results obtained can be explained by weak sensitivity of criterion (7) to parameter 0 , that can be seen from curves of isolines in Fig. 5 (b). Fig. 5 (c) shows the output time characteristics of the reference Y t and the identified Y t model after identification. In the graph, the variables of the identified model output vector Y t are in the region of the reference model vector Y t with random perturbation. It can be concluded that the proposed identification algorithm has good convergence in the stepped input influences case.
GA with the input Xt 1, 1, 0.2, 1T ; а) Transients modelling before identification;
to criterion (7);
and GA
However, under real-world conditions, the input can be of any shape. Consider the more complex case where the input signal has a harmonic component. In this case, the input action Xt 1, 1 0.2sin(0.2t), 0.2 0.1sin(0.3t), 0.5 0.5sin(0.1t)T is applied to both models. Fig. 6 (a) shows the simulation case when the parameters of the identified model 0 300, G 0.25 differ A significantly from the reference ones 0 161, GA 0.43. Fig. 6 (b) shows surface isolines of criterion (7) and its minimization trajectory. In the process of identification parameters values 0 161.9, G 0.426 are optimized. Numerical values of identification parameters are sufficiently close to
A reference ones. Fig. 6 (c) shows temporal characteristics of the reference vector Y t and identifiable Y t model after identification under harmonic change X(t). It can be concluded from the simulation results that the proposed passive identification algorithm has good convergence in the harmonic component presence the vector Xt and the random noise presence.
a)
under harmonic variation of the input vector Xt ; а) Transients modelling before identification;
to criterion (7);
and GA
The development of the digital twin mathematical model will be based on the continuous model (3), which previously passed the uncertain parameters identification stage using control samples of trends X(t) and Y(t) for the operating electric heater. The digital twin model must reflect the physical process dynamics in real time. To ensure this condition, the mathematical model must be discrete, with the sampling time ensuring the information distribution over CFS network.
The transition from the continuous model (3) to the discrete analogue can be obtained by using ztransformations [17]
Y(z) W z X(z) , (8) here W ( z)
z 1 the multiplier
z quantization moments TKV .
z 1 z
W ( p) Z is the discrete transfer matrix of the multidimensional system, where p indicates the presence of zero-order extrapolator for fixing the signal between Thus, the application of mathematical model for s digital twin for electric heater consists of steps. 1. The uncertain parameters (0 and GA) identification of electric heater mathematical model (3) by the considered algorithm.
2. Transition from continuous model (3) to discrete model (8), which is digital twin. 3. If during operation, the digital twin accuracy has deteriorated due to the non-stationarity of the physical process, it is necessary to go to step 1 to determine the model parameters.
In accordance with the proposed methodology, the digital twin mathematical model simulation of the electric heater was carried out using the MatLAB software package. The reference model (3) with transfer matrix (4) was used as a basic model. The function c2d(...) of the MatLAB software package was used to calculate the matrix W ( z ) of digital twin (8). The simulation results are shown in Fig. 7.
the
In Fig. 7 (a) simulates the case with initial conditions for the reference model 0 161, G 0.43 and A identifiable model 0 120, GA 0.22 and at the input step action Xt 1 0.5sin(0.15t) 1 0.2sin(0.15t) 0.2 0.1sin(0.3t) 0.2 0.8sin(0.05t)T . In Fig. 7 (b) shows the surface isolines of criterion (7) and parameters minimization trajectory , GA, resulting in 0 154, 0 G 0.428. From these parameters the model transfer matrix being identified is calculated A Y . Based on the simulation results, it can be concluded that the proposed methodology makes it possible to implement sufficiently exact the digital twin of air heating.
A distinctive feature of the proposed digital synthesis technique is the model parameters identification that are imprecisely defined at the model development stage and are refined in the passive identification process. As a rule, in modern identification methods, the model structure is specified and all parameters are determined from its dynamics [17, 18]. In the classical case, more software resources are required for identification, and the identification result may be unsatisfactory if the model structure is incorrect. In the proposed methodology, the model structure is known, for which only uncertain parameters are identified and not the model as a whole.
The introduction of real-time industrial control systems leads to the complex mathematical models use with uncertain parameters and the complexity of their identification tasks. As an example, the methodology of analytical model passive identification of electric heater with subsequent digital twin synthesis for CFS is considered. For the considered model changing parameters analysis is made and their influence on modelling results is demonstrated. An algorithm for passive identification of the mathematical model uncertain parameters is proposed that uses criterion (7) to numerically optimize the uncertain model coefficients values using the accumulated trend base of the physical process.
Several examples of parameter identification 0 and GA for the analytical model (3) are demonstrated. The examples show that the uncertain model coefficients identification should be classified as a one-extremes optimization problem. Application of the mathematical model for the electric heater digital twin CPS is proposed and numerically investigated. The simulation results confirmed the effectiveness of the considered digital twin design methodology using MatLAB application package. In the future, using available data trends of reference model (3), it is supposed to synthesize electric heater digital twin based on ARCS model and analyze the obtained results.
The use of digital twins in CPS makes it possible to identify process bottlenecks, improve product quality, and reduce the risks of abnormal operation throughout its equipment lifecycle. Digital twins can be used to optimize equipment modes, predict and detect faults, and find-modifications to the structure of a real physical system based on its observed effects in the future. This approach provides a highly accurate assessment of a plant's production capacity when drawing up a production program.
[3] M. Grieves, Completing the Cycle: Using PLM Information in the Sales and Service Functions [Slides], SME Management Forum, Troy, MI. (2002). [4] M. Grieves, Can the digital twin transform manufacturing, World Economic Forum, 2015. https://www.weforum.org/agenda/2015/10/can-the-digital-twin-transformmanufacturing/?DAG=3&gclid=Cj0KCQjwwtWgBhDhARIsAEMcxeCFyfAQIbgE1WIoIQUT1 wAeywKBphlxnxPUaNBTN5mtDPykV7GLwlcaAho2EALw_wcB [5] V. Slyusar, The concept of networked distributed engine control system of future air vehicles, Proceedings of AVT-357 STO NATO Workshop on Technologies for future distributed engine control systems (2021) 1-12. doi:10.14339/STO-MP-AVT-357 [6] E. Negria, L. Fumagalli, M. Macchia, A review of the roles of Digital Twin in CPS-based production systems, 27th International Conference on Flexible Automation and Intelligent Manufacturing (2017) 939-948. doi:10.1016/j.promfg.2017.07.198 [7] Digital Twins for Industrial Applications, An Industrial Internet Consortium White Paper, 2020. https://www.iiconsortium.org/pdf/IIC_Digital_Twins_Industrial_Apps_White_Paper_2020-0218.pdf [8] B.A. Talkhestani, T. Jung, B. Lindemann, N. Sahlab, N. Jazdi, W. Schloegl, M. Weyrich, An architecture of an intelligent digital twin in a cyber-physical production system.
Automatisierungstechnik, (2019) 762–782. doi:10.1515/auto-2019-0039. [9] P.I. Bidyuk, O.L. Tymoshchuk, A.E. Kovalenko, L.O. Korshevniuk, Decision support systems and methods. Kyiv, Igor Sikorsky Kyiv Polytechnic Institute, 2022. [10] E.A. Lee, Fundamental Limits of Cyber-Physical Systems Modeling. ACM Transactions on
Cyber-Physical Systems. (2016) 1-26. doi:10.1145/2912149 [11] E.A. Lee, The Past, Present and Future of Cyber-Physical Systems: A Focus on Models, Sensors (2015) 4837–4869. doi:10.3390/s150304837 [12] N. Pankratova, I. Golinko, Electric heater mathematical model for cyber-physical systems, System research and information technologies (2021) 7–17. doi:10.20535/SRIT.23088893.2021.2.01 [13] Y.M. Kornienko, Yu.Yu. Lukach, I.O. Mikulonok, V.L. Rakytskyi, G.L. Ryabtsev, Processes and equipment of chemical technology, NTUU "KPI", Kyiv, 2012. [14] T.O. Prokopenko, Classification of uncertainties in the management of organizational and technological objects, Technological audit and production reserves (2014) 23–25. doi:10.15587/2312-8372.2014.30336 [15] S.D. Shtovba, Optimization methods in the MatLab environment, VDTU, Vinnytsia, 2001. [16] N. Pankratova, P. Bidyuk, I. Golinko, Decision support system for microclimate control at large industrial enterprises, CMIS-2020: Computer Modeling and Intelligent Systems (2020) 489–498. https://ceur-ws.org/Vol-2608/paper37.pdf [17] V.M. Dubovoi, Identification and modeling of technological objects and control systems, VNTU,
Vinnytsia, 2012. [18] P.I. Bidyuk, O.M. Trofymchuk, A.V. Fedorov, Information system of decision-making support for forecasting financial-economic processes based on structural-parameter adaptation of models, KPI Science News (2011) 42-53.