Development of an Operational Digital Twin of a Locomotive Parking Brake for Fault Diagnosis

In recent years, a growing role in digital technologies has been filled by model-based digital twinning. A digital twin produces a one-to-one mapping of a physical structure, operating in the digital domain. Combined with sensor technology and analytics, a digital twin can provide enhanced monitoring, diagnostic, and optimization capabilities. This research harnesses the significant capabilities of digital twining for the unmitigated challenge of fault type classification of a locomotive parking brake. We develop a digital twin of the locomotive parking brake and suggest a method for fault type classification based on the digital twin. The diagnostic ability of the method is demonstrated on a large experimental dataset.

mathematical models using sensor data 28 .These algorithms can address the differences between simulation and reality 29,30 and, hence, provide DT improvements.
In this study, a model-based DT approach for classifying various faults that may occur in a locomotive parking brake is proposed and validated by experiments.The DT is based on a physical model of the parking brake and is optimized using a machine learning approach.The DT of the locomotive parking brake is continuously updated to generate possible fault conditions that approximate those of the actual system.The study has three major contributions: 1. Development of a DT of a locomotive parking brake.2. Development of a robust diagnosis of various faults in a locomotive parking brake by estimating the internal latent physical variables within a DT and training of a learning model on the residuals signals and model based estimated internal variables.3. Show how deep physical understanding and model-based DT can improve the diagnostic capabilities of complex systems, especially in scenarios where machine learning algorithms alone may not be sufficient.
The study is divided into six sections.Section "Theoretical background" provides a theoretical background and introduces the locomotive parking brake and the physical model DT, Section "Experimental validation of the proposed parking brake model" validates the physical model DT using experimental data, Section "The proposed algorithm" presents the new algorithm based on DT, and Section "Fault diagnosis" demonstrates the algorithm using experimental data.The study is summarized in Section "Summary and conclusions".

Parking brake working principle
The JT42 locomotive is equipped with a spring-loaded brake.Four of the eight brake blocks are equipped with load springs, as shown schematically in Fig. 1.These four brake blocks contain in the rear part of their housing the following main components: a piston, a load spring, and a manual release device.The cylinder part equipped with the load spring has an independent compressed air connection; the piston chamber is hermetically isolated from the service brake chamber and is additionally separated by an atmospheric sluice.In the released state, the loading spring cylinder is filled with compressed air, and the piston and connecting rod are retracted against the resistance of two springs until no force is applied against the piston base of the brake shoe (in the front part of the cylinder).In this position, the springs of the parking brake are loaded and energy is stored in them.In the applied (brake operative) state, the piston of the spring storage cylinder chamber is empty of compressed air.The force of the springs press on the piston of the service brake cylinder via the connecting rod, which presses the brake pads onto the wheel.
To apply or release the parking brakes, an electrical switch is located in each driving cab.A relay controls the one-way solenoid valve and allows compressed air from the main reservoir line to flow into the spring chamber through a pressure-reducing valve set at 7.7 bar and through an orifice in the parking brake valve.This compresses the loading spring in the spring chamber and the parking brakes remain in the non-operated mode.A healthy parking brake should be completely released after 25 s.By switching off the solenoid valve, the compressed air is released from the spring chamber so that the brake shoes are driven against the wheel by the force of the spring.For more detailed information on the braking system, see 31,32 .

Parking brake model
Physical systems must first be represented in terms of a mathematical equation for each machine before any particular analysis can be performed for those particular systems.Mathematical modeling is required before any control over the systems can be explored.In general, the modeling of each system can be created either from www.nature.com/scientificreports/physical derivatives or experimental data measured on the real system.This paper presents nonlinear mathematical modeling of a JT42 locomotive parking brake based on readily accessible physical derivatives.The mass flow rate, pressure dynamics, and equation of motion were derived from previous research 30 .
The compressible mass flow rate, ṁ through a valve orifice, can be described as 1 ; where P u , P d are the upstream and downstream pressure respectively, C f is the non-dimensional discharge, T temperature and k is the specific heat ratio, and ṁ(P u , P d ) is the mass flow rate.The constants C 1 and C 2 are given by; where R is the gas constant, P cr is a critical pressure related to the ratio between downstream pressure P d and upstream pressure P u .The value of P cr can be expressed as; The upstream and downstream pressures of the cylinder are different, depending on whether the cylinder chamber is charging or discharging, according to the following functions: where ṁin , ṁout are the mass flow rate into and out of cylinder chamber A respectively, P chamber is the pressure inside the cylinder chamber, P s is the supply pressure and P a is the ambient pressure.
Equation (5) represents the charging process, in which the pressure in the reservoir is considered to be upstream and the pressure in the cylinder chamber is downstream.For the discharging process, represented by Eq. ( 6), the pressure in the chamber is the upstream pressure and the ambient pressure is the downstream pressure.The flow condition can be classified as either choked flow or under-choked flow, depending on the downstream pressure P d and the upstream pressure P u of the orifice.The flow condition is considered to be choked flow when P d P u is less than the critical pressure ratio P cr .Assuming that the gas is perfect and the process is considered adiabatic, the rate of pressure change in a pneumatic chamber can be expressed as follows; where P(t) is the pressure inside the chamber, ∂ ∂t P(t) rate of change in pressure inside chamber, ∂ ∂t m in (t) inlet mass flow rate to chamber, ∂ ∂t m out (t) outlet mass flow rate from the chamber,V (t) is the volume of the chamber, and ∂ ∂t V (t) rate of change in volume of the chamber.The origin of the piston placement can be chosen either from one side at the end of the stroke or in the middle of the stroke.By choosing the origin on one side at the end of the stroke, the volume of the chamber can be expressed as follows: where x(t) is the piston displacement and A p is the area of the piston.
By differentiating Eq. ( 8), the rate of change of the volume of the chamber is obtained as follows: Substituting Eqs. ( 8), ( 9) into Eq.( 7) gives the rate of pressure change in the chamber of the pneumatic cylinder: where x(t) is the piston position, ∂ ∂t x(t) is the rate of change of piston position, and V 0 is the inactive volume at the end of stroke and admission port.
(1) The equation of motion for the piston rod, including the mass and friction effects of the pneumatic cylinder, can be expressed by Newton's second law as follows: where M p is the mass of the piston rod, F f is the friction force, P a is the atmospheric pressure, A r is the cross- sectional area of the piston rod and K is a constant factor characteristic of the spring.
The LuGre friction model proposed by Liu et al. 33 is considered in this study.The frictional force, F f , the dynamics of the internal state ( z ), and the Stribeck effect function g(v) are given in Eqs. ( 12), (13), and ( 14), respectively.From these equations, the relationship between velocity and friction force for the steady-state of motion can be derived as Eq. ( 15): where σ 0 is a stiffness coefficient, σ 1 is a damping coefficient, B is the viscous friction, F c is the Coulomb friction, and F s is the static friction.

Experimental validation of the proposed parking brake model
The parking brake DT is configured to simulate three typical fault types, including cylinder leakage, blocked upstream pressure, and low upstream pressure.A total of five working states of the parking brake system are simulated using a DT, including one health state, two types of single fault states, and two types of composite fault states as shown in Table1.The model parameters were calculated assuming the following; the total heat energy of the gas does not change during compression.The specific heat coefficient for air is k = 1.4 .The value of the gas constant R , temperature T and ambient pressure Po at a relative humidity of 65% are 288 J/kgK, 293.15K and 100 kPa , respectively (Standard ISO 6358).This gives constants C 1 and C 2 of 0.040418 and 0.156174 , respectively, which were calculated using Eq.(3).Equation ( 4) is used to set the critical pressure ratio p cr to 0.528 .The pressure after the pressure-reducing valve was set at 7.7 bar .The cylinder bore diameter is 0.1 m.
To verify the proposed mathematical model, a series of experiments was conducted in which the pressure in the cylinder chamber and after the pressure-reducing valve was were measured using piezoelectric absolute pressure transducers, as can be seen in Fig. 2.
The experiments were conducted as follows: First, the parking brake was applied (i.e. the pressure in the cylinder chamber is 0 bar), and then the parking brake was released using the electric switch in the driver's cab. Figure 3 shows both the model and experimental results for the pressure in the cylinder chamber.There is close agreement between the theoretical and experimental curves, with very good amplitude agreement.Figure 3 also shows a time delay between the valve command (0 s) and the corresponding pressure output.This phenomenon could be due to the connecting pipes and the dynamics of the actuator valve.In the first stage the pressure in the cylinder chamber rises rapidly to 4.07 bar.Then, as the piston begins to move, the rate of change of pressure decreases as the volume of the cylinder chamber increases.Finally, when the piston reaches its end position (25 s), the pressure rate of change is zero and the pressure in the cylinder chamber is equal to the pressure of the pressure-reducing valve, i.e. 7.7 bar.The results show that the DT simulates the parking brake with great accuracy, which reflects the quality of the mathematical procedures used.
Based on the verified DT model, different failure states of the parking brake system can be simulated by adjusting some parameters of the model.Figure 4a shows the results obtained when an experiment was performed with an upstream pressure of 7.7 bar illustrates a healthy parking brake system (red curve).The figure Working conditions of parking brake.

Working states of the parking brake Total number of samples
Normal state 500,000 Air leakage 500,000 Blocking due to an incorrectly set orifice 500,000 Cylinder leakage & blocking 500,000 Low upstream pressure 500,000 also shows the simulated curves of the pressure in the cylinder chamber obtained using the DT for a scenario where the outlet (upstream) pressure of the pressure reducer is lower than 7.7 bar.In all scenarios excluding the case where the upstream pressure is lower than the set pressure (i.e.4.027 bar, yellow), there is close agreement between shape of the curves.The results in Fig. 4a shows that in the case where the upstream pressure is lower than the set cylinder pressure, there is a rapid increase in the pressure inside the cylinder up to the set pressure.From this point it remains constant (4.027 bar).This means that the chamber volume does not change because the force acting on the piston is too small (the piston remains in place, i.e. the parking brakes do not release).
Figure 4b shows the results of an experiment in which the upstream pressure for both DT and the real parking brake system (red curve, illustrates a healthy parking brake system) was set to 7.7 bar.In this scenario, blocking of the upstream pressure due to an incorrectly set orifice was simulated using the DT. Figure 4b shows that the pressure in the cylinder chamber quickly rises (3 s) to the upstream value (7.7 bar) when no orifice was inserted into the parking brake valve outlet (yellow).In this case, the piston moves violently towards its end position.This situation can cause severe damage to various parts of the parking brake or even to the entire mechanism.In contrast, too small and orifice can cause a significant delay in releasing the parking brake.The results in Fig. 4b  show that the rate of change of pressure in the cylinder chamber using an orifice of 0.4 mm (dotted curve), for example, is small compared to the healthy condition (red curve, illustrates a healthy parking brake system).In this case, the pressure in the cylinder chamber reached the set pressure value (4.027 bar) only after 50 s, causing a significant delay in the parking brake release application.Figure 4c shows the results of an experiment in which the upstream pressure drops due to leakage in the cylinder or in the connecting pipes.The figure shows the different rates of pressure drop in the cylinder chamber found experimentally and simulated using DT compare to a healthy parking brake system (red curve).Some curves represent a rapid pressure drop, simulating a severe leak typical of a torn pipe, while other curves (yellow) represent a pressure drop reminiscent of a small leak due to improper sealing.Figure 4c shows that a severe leak due to a torn pipe causes the air pressure in the cylinder chamber to drop rapidly, causing the piston to move uncontrollably toward its initial position.Such a movement can severely damage the brake pad.On the other hand, a leak resulting from improper sealing can be overcome by the supplying pressure to the cylinder.However, this results in a considerable delay in releasing the parking brake.
To evaluate the capacities of the DT, the relative Normalized Mean Squared Error (NMSE) between the measured and simulated pressures in the cylinder chamber at each time point was calculated as follows: where N is the total number of time points, P A is the measured signal, P A is the DT simulated signal and i the index represents the time points.
To calculate the NMSE, the same experiments simulating the different scenarios (see Table 1) were performed with six locomotives.Figure 5 shows the NMSE comparing each DT simulated scenario and the corresponding measured scenario.For example, the leftmost column shows the NMSE index between the measured and simulated air pressure in the cylinder chamber at normal state.The highest error was calculated for the experiment simulating an air leak in the cylinder together with an air blockage.The relatively high error can be explained by the relative complexity of the experiment and the sensitivity of the different parking brake valves.However, Fig. 5 shows that the mean square error is relatively small (highest value, approx.5%) for all simulated scenarios, reflecting the quality of the DT.

The proposed algorithm
The proposed algorithm is integrated in the locomotive DT.It consists of four steps, illustrated in Fig. 6.The pressure in the parking brake cylinder (marked by P A ) was routinely measured on the tested locomotives, as can be seen in Fig. 6.In the current study, the internal latent physical variables ( K, B, σ 0 and σ 1 ) are estimated by solving a least squares problem, where the variables minimize the constraints presented in Eq. ( 17).This is achieved in Eq. ( 18), also known as a least squares estimation: where x k [n] represents the corresponding values of the internal parameters θ k in the coordinate n e.g., for θ k in time t the corresponding value is P A t �t , y[n] represents the measured parameters in the coordinate n and X is the matrix of x k [n] with N rows and j − i + 1 columns.
The internal variables ( θ 1 , θ 2 , θ 3 and θ 4 ) can be used to calculate the process coefficients presented in Eq. ( 19): For each Real Twin (RT) (i.e., physical measurement P A ) in the training set, an individualized DT is gener- ated by estimating the internal variables ( θ 1 − θ 4 ) .These variables are estimated by the least squares method presented in Eq. (17).
Based on the internal estimated parameters, a model of Deep Neural Network (DNN) is trained on the estimated parameters where, at first, the training set is divided into 80% training and 20% validation.The DNN model uses the sigmoid activation function in the input and hidden layers, while the Softmax function is used in the output layer.The selected optimization algorithm was Adam with an initial learning rate of 0.001 and the loss was categorical cross entropy.The DNN model included an input layer, eight hidden layers with eight neurons each, and one output layer.The number of epochs was regulated using early stopping based on three followed unimproved loss on the validation set.
For each Real Twin (RT) (i.e., physical measurement P A ) in the test set, an individualized DT is generated by estimating the internal variables ( θ 1 − θ 4 ) .These variables are estimated by the least squares method presented in Eq. (17).The trained model is used to predict the classes of the test set.

Fault diagnosis
The new algorithm described in Section "Experimental validation of the proposed parking brake model" is tested and compared with two other algorithms: 1.A regular machine-learning algorithm using a DNN model consists of Steps 3, 4, and 5 described in Fig. 7.
The regular machine-learning algorithm was trained with 21 different possible features: Mean, Variance, Maximum, Kurtosis, standard deviation, skewness, Absolute Sum etc., of each RT's features extracted directly from the measured signals, i.e., P A .The DNN model is trained on these extracted features, as described in Fig. 7.This algorithm does not use the DT. 2. A simple residual algorithm which uses Steps 1, 2, 3, 4 and 5 (including the DT) described in Fig. 7.For each RT in the training set, a DT is generated by estimating the internal variables, as explained in Fig. 7. Based on the internal estimated parameters, the DT calculates the residuals.From each residual, 21 features are extracted: mean, variance, maximal value, kurtosis, absolute sum etc.A model of deep neural network is trained on these extracted features, as described in Fig. 7.
The comparison process between the new algorithm and the regular machine-learning algorithm demonstrates the contribution of the DT concept, while the comparison to the simple residual algorithm demonstrates the contribution of the internal variables estimation.These three algorithms, i.e., the new algorithm described in Figure 6.The new suggested algorithm for machine fault diagnosis using a DT described in Section "The proposed algorithm".
Section "Experimental validation of the proposed parking brake model", the regular machine-learning algorithm, and the simple algorithm, were tested on an experimental dataset containing a 5500 RT training set and a 250 RT test set.An example of the measurements of an RT in the training set is depicted in Fig. 3.
The results of the three tested algorithms are presented in Fig. 8 on 250 test examples, 50 from each condition.Each table in Fig. 8a-c presents the confusion matrix after applying the tested algorithm on the test set, which contains 50 examples of each category (overall, 250 examples).As can be seen in Fig. 8d, the new algorithm achieved a significant improvement from an accuracy of 72-96.Furthermore, the use of the residuals helps to improve the accuracy from 72 to 81.2, as can be seen in Fig. 8d.This result demonstrates the ability of the DT to improve diagnosis by incorporating physical knowledge of the system.

Summary and conclusions
In this study, a DT algorithm for diagnosing faults in a locomotive parking brakes is presented.The algorithm comprises four steps, including estimation of the internal DT variables, training of a deep neural network, and prediction.The algorithm was tested with a dataset of 5500 training RTs and 250 test RTs.The results show a significant improvement in accuracy from 72 to 96 compared to a traditional machine-learning algorithm.It is shown that the DT improves diagnosis by incorporating physical knowledge about the system.The described method for model-based fault detection and diagnosis for locomotive parking brakes is based on standard measurement of the pressure in the parking brake cylinder.A key advantage of the described DT based fault detection method is that it uses a physically-derived parking brake model.Therefore, it is applicable to a wide range of operating points and directly transferable to other parking brakes.In addition, its symptoms are usually easy to interpret and understand.The model-based DT approach presented here provides an application of engineering performance that goes beyond traditional engineering design.In addition to impacting train operation, it can also play a role in improving maintenance activities in which maintenance failures are effectively diagnosed and corrected.Ultimately, the most important benefit of the DT is its impact on customer experience and operating costs.This paper demonstrates the value of performance-focused engineering where real-time performance provides inputs to an adaptable system packaged in a digital layer-the DT.

Figure 2 .
Figure 2. Schematic diagram of the parking brake system and experimental setup.

Figure 3 .
Figure 3. Model (blue curve) and experimental (red curve) results for the pressure in the cylinder chamber of healthy parking brake.

Figure 4 .
Figure 4. Model results show the response of the pressure in the cylinder chamber for different faults compared with healthy experimental results (red curve).(a) low up-stream pressures, (b) different sizes of the orifice, (c) sudden air leak in up-stream pressure.

Figure 5 .
Figure 5. NMSE index showing the mean square error between the experimental and simulated pressures in the cylinder chamber for different scenarios for six locomotives.

Figure 7 .
Figure 7.The two other algorithms use as comparison to the new suggested algorithm described in Section "Experimental validation of the proposed parking brake model".

Figure 8 .
Figure 8. Results of the new suggested algorithm.(a) simple residual (the new algorithm with the DT), (b) regular machine-learning algorithm without the DT and (c) the new suggested algorithm.(d) The summarized accuracies of (a-c).