Multirate DAE-simulation and its application in system simulation software for the development of electric vehicles

This work is devoted to the eﬃcient simulation of large multi-physical networks stemming from automated modelling processes in system simulation software. The simulation of hybrid, battery and fuel cell electric vehicle applications requires the coupling of electrical, mechanical, ﬂuid, and thermal networks. Each network is established by combining the connection structure of a graph with physical equations of elementary components and resulting in a DAE (diﬀerential algebraic equation). To speed up the simulation a non-iterative multirate time integration co-simulation method for the system of coupled DAEs is established. The power of the multirate method is shown via two representative examples of battery powered electric vehicle with a cooling system for the battery pack and a three phase inverter with a cooling system. This work is an extended version of (Kolmbauer et al. in Scientiﬁc Computing in Electrical Engineering (SCEE 2020), Springer, Cham, pp. 231–240, 2021).

lubrication system of the ICE and battery packs, the electrical propulsion system including the engine and a battery pack, the air conditioning and passenger cabin models, cf. [5], waste heat recovery and finally according control systems. Due to the complex interaction of the subsystems, the challenges in the development of future power trains do not only lie in the design of individual components but in the assessment of the power train as a whole. On a system engineering level, it is required to optimize individual components globally and to balance the interaction of different subsystems. Due to the increasing complexity of the models, the systems exhibit largely varying time scales and are difficult in the numerical simulation. A mainly automatized multirate approach is a promising way to decrease the computational effort. One possible implementation for this is the time integration via a co-simulation approach, cf. [6][7][8].
The structure of the work is the following: In Sect. 1 the individual physical networks are introduced, and the coupling conditions are stated in order to obtain a fully coupled system of network DAEs. The multirate time integration technique for the coupled system of network DAEs is described in Sect. 2 and the corresponding numerical results are presented in Sect. 3 and Sect. 4. Finally we conclude in Sect. 5.

Problem formulation
In the following, a network is discussed that consists of several multi-physical elements. The network elements describing the electric contribution are provided by current sources, voltage sources, nodes, ground, resistors, capacitors, and inductors. The fluid network consists of pipes, pumps, demands, junctions, and reservoirs. The electrothermal coupling is established by lumped mass elements representing the pipe wall and the masses from the battery and heat transfer connections. The individual components are assembled to a network N , which is represented by a linear directed graph. The graph structure is described by an incidence matrix A, which can be used for the model descriptions, cf. [9]. In the following we state the DAEs for the three main involved physical networks.
Electrical network The electrical network N E = {R, C, L, V , I, N, G, B} is composed of resistors R, capacitors C, inductors L, voltage sources V , current sources I, nodes N , grounds G and batteries B. The DAE for the network in N E in input-output form is given by: For T and the outputs y = y R , such that for given boundary conditions e G = 0 and predefined resistance r, capacitance c and inductance l as well as prescribed currentsj I and prescribed voltagesv V andv B . The coupling variables are expressed as temperature of the resistor u R , the capacitor u C and the battery u B as well as the energy flux of the resistor y R , where • denotes the Hadamard product, i.e. the elementwise vector product, and | · | the elementwise absolute value.

Multi-physical model
The multi-physical model is derived by combining (1), (2) and (3) with appropriate coupling conditions. The coupling conditions describe the relation between the inputs and outputs of the individual models. For the model used in Sect. 3 and Sect. 4, the following coupling conditions are used, see e.g. [2].
The connectivity equation (4) represents the electro-thermal coupling of the electrical network and the cooling systems. Combining all subsystems and their connectivity equations (4) yields a DAE: Find DAEs resulting from automated modelling software typically obtain a structure with dindex (differentiation index) greater than 1, cf. [3,4,10] and hence are not suitable for a direct simulation with standard solvers. In the setup of multiple physical networks, it is not sufficient, that the full DAE (5) can be reduced to a d-index 1. Additionally, each subsystem, to which a solver is applied, has to fulfill d-index 1 conditions as well, cf. [7,11]. In our applications an automatic index reduction is performed if the electric or the fluid system are detected to be of d-index 2. This is achieved by providing surrogate models according the corresponding literature, cf. [3,4], which enables additionally a consistent definition of initial values. Similarly, it is possible to describe the coupling of thermal solid systems with gas systems. An approach for coupling fluid, gas and thermal solid systems can be found in [5].

Multirate integration for coupled network DAEs
In our multirate approach the full DAE (5) is partitioned due to the physical background to n ∈ N subsystems (typically n 2). Each subsystem is index reduced according to the available literature, cf. [3,4,10]. Since in the global network the individual subsystems are interacting with each other, i.e. inputs and outputs are connected according to the connectivity equation (4), it is necessary to put it into an input-output form. For this purpose, each subsystem i = 1, . . . , n classifies its inputs u i , state variables x i , algebraic variables a i and outputs y i . To conclude, this approach yields a coupled system of n semi-explicit DAEs in input-output form of d-index 1. For inputs u i specified by equation (4), find x i ,ẋ i , a i and y i , such thaṫ for i = 1, . . . , n. A careful choice of the connectivity matrix given in (4) guarantees that the coupled system obtains d-index 1 as well, cf. [2,7]. E.g. one possible choice is the usage of differential states, which are not involved in any index reduction, as coupling variables. For each subsystem (6) an arbitrary Runge-Kutta method with micro-step sizes h i is used, cf. Fig. 1, in an indirect approach, cf. [12]. This means that just the differential part is solved using a Runge-Kutta scheme, while the algebraic part is solved for each given time and differential state. The choice of the actual integration technique depends on the properties of the underlying system and can be explicit, implicit, fixed, or adaptive. The whole system (5) is integrated via a non-iterative co-simulation technique with macro-step size H = max(h i ). All systems are updated at the end of each macro-step. This principle relates to synchronous communication and we refer to these points in time as synchronization times, cf. Fig. 1. In the general case, the individual time steps vary significantly, which means that a global synchronization time quickly pushes the multirate approach to its limits. A hierarchical modelling approach allows for defining individual synchronization times for coupled subsystems. A partitioning algorithm detects co-simulation components that can exchange information on smaller synchronization step sizes. By Figure 1 Macro-step of the ith system from synchronization time t k to t k+1 grouping on multiple levels, multiple synchronization times can be configured. This allows higher frequency subsystems to be exchanged on lower time scales. The evaluation of each macro-step of the subsystems is done in a sequential Gauss-Seidel-approach. The values u i are handled with appropriate interpolation or extrapolation techniques, depending on the slow or active characteristic of the interacting subsystems. To support energy conservation across circuit borders, fluxes are integrated and interpolated accordingly such that incoming and outgoing energies are identical. It is important to mention that for our multi-physical model (5), we always have cyclic dependencies between the coupled physical domains. For example, the temperature and the heat flux are exchanged when coupling electric and solid systems. Therefore, extrapolation is always needed at some point. Nevertheless, it is possible to state convergence order results of the method for slow-first, fast-first and fully decoupled co-simulation approaches for n = 2, cf. [7]. Due to n 2 this approach is generalized to a more general one. This is done by establishing an evaluation order for all subsystems. By defining dependencies within and between the subsystems, it is possible to analyze the evaluation order based on a graph on multiple levels. According to this graph, an algorithm computes a valid evaluation order. In most relevant cases, this satisfies the conditions, cf. [7], such that the co-simulation approach has convergence order p if the Runge-Kutta method has convergence order p and the interpolation is of order p -1. This applies to the examples presented in Sect. 3 and 4.

Simulation of a BEV with cooling system
The given BEV example demonstrates the modelling of an electrical system linked with the required cooling system, cf. Fig. 2. The model consists of an electrical propulsion and two cooling circuits. An oil circuit is used for cooling of the electric machine and a coolant  Table 1 for BEV with cooling system simulation circuit is used for cooling of the battery pack, inverter, and DC-DC converter. The involved subsystem of the coupled electro-thermal model can be reduced to DAEs of d-index 1.
The multirate approach presented in Sect. 2 is put into context with the reference solution of a single solver approach (both sequential/single CPU). In this example eight thermal circuits of the form (2), three mechanical circuits, an electric circuit of the form (1), 14 gas circuits given in [5] and two fluid (oil and cooling and lubrication) circuits of the form (3) are present which represent in total 461 equations. The preferred choice of solvers for both, the single solver approach and the multirate approach are adaptive explicit solvers [13] but are compared for completeness with a fixed step singlerate approach as benchmark solution. Hence the step size of the single solver is limited to the minimum step size of all subdomains, while the multirate approach is limited to the synchronization time and the characteristic of the global domain. Four different solver parametrizations were simulated and compared. The first parametrization is a fixed step singlerate solver with a time step of 1 ms. We refer to this in the charts and tables as Singlerate (Fixed: 1 ms). The second one is an adaptive singlerate solver with a maximum time step of 50 ms, referenced by Singlerate (Adaptive: 50 ms). In this example no singlerate parametrization was found that achieves a better performance. On the other hand, we have two multirate solvers which use the benefit of adaptive solvers. In the first approach, which is called Multirate (Synctime: 20 ms), the electric, gas and oil circuits use a maximum step size of 10 ms while cooling and lubrication, mechanical and thermal domains use a maximum step size of 20 ms. Here the information is exchanged at a global synchronization time of 20 ms. The last parametrization (Multirate (Synctime: multiple)) uses a maximum step size of 10 ms again for the electric, gas and oil circuits. For mechanical and thermal it uses 20 ms and for the cooling and lubrication it uses 40 ms as a maximum step size. In this constellation two interacting subsystems derive their synchronization time from the highest maximum step size of both parametrizations. Thus, synchronization times of 10 ms, 20 ms and 40 ms are used. Since many physical subsystems are involved, the overall evaluation sequence of the individual physical domains is complex to visualize. We give a simplified list of the evaluation order with the main contributions: 1. Cooling and lubrication circuit 2. Oil circuit 3. Electrical circuit 4. Mechanical circuits 5. Gas and thermal solid circuits This evaluation order is valid for both the Multirate (Synctime: 20 ms) and Multirate (Synctime: multiple) approaches. In the Multirate (Synctime: multiple) approach the oil and electric circuit are exchanging data at synchronization time 10 ms and form a new group. This itself is coupled to the mechanical, gas and thermal solid circuits with a synchronization time of 20 ms and again creates a new group. At the end for the coupling with the cooling and lubrication circuit a synchronization time of 40 ms is applied.
A comparison of the results in the heat flow between the battery and the cooling system shows that the two singlerate solutions differ from the multirate solutions, cf. Fig. 3. Because of the information exchange on macro step size, multirate approaches are more inert in the simulation. For real-time simulations, such deviations in results are usually not of major importance. It is more important that the overall model is calculated in an energy conserving manner to avoid unphysical effects and unrealistic controller behaviours. This Figure 3 Comparison of heat flow between battery pack and cooling system of multirate cases against single solver cases for BEV with cooling system simulation Figure 4 Comparison of the energy loss of the electric system of multirate cases against single solver cases for BEV with cooling system simulation is achieved by using flux averaging techniques similar to the more advanced approach in [14] for mechanical couplings. The verification of the energy conservation is mainly done empirically and with a-posteriori consistency checks. The energy loss in the electrical system shows that there is no difference between singlerate and multirate solvers, cf. Fig. 4. The multirate solvers provide sufficiently accurate results.
Besides accuracy, performance is an important criterion for real-time applications. Both, simulation time and real-time factor (RTF) are interesting characteristics to assess the performance. The real-time factor is calculated based on the computation time This value can be used to evaluate the real-time capability of a simulation. In Table 1 the simulation time and average RTF of the singlerate cases is compared with those of the multirate case, cf. Fig. 5. A significant speed up in the calculation time can be achieved and one can see that this is achieved per time step not only in some critical areas. Using multiple synchronization times, another 10% of the computing time can be saved compared to the standard multirate approach. In total, the computation time can be reduced to a fifth of the best singlerate case.   Table 2 for BEV with cooling system simulation

Simulation of a three phase inverter with cooling system
We consider a detailed physical model of an inverter with transitors, IGBT (insulated-gate bipolar transistor), switches, an RC (resistor-capacitor) filter as well as a 3 phase ohmic load. The inverter is used to convert a DC (direct current) voltage through timed switching of the six transistors into a PWM (pulse width modulation) signal. The RC filter then averages the PWM and thus creates a 3 phase AC (alternating current) voltage, cf. Fig. 6.
In total this example consists of 178 equations which are spread over 20 solvers. A fluid circuit described by (3), seven gas circuits modeled like in [5] and eleven solid thermal circuits of the form (2) are responsible for modelling the cooling of the electric system (1).
In this example three solver constellations are compared. The reference solution (Singlerate (Fixed: 1 μs)) of the system is using an explicit fixed step solver with 1 μs. In the first multirate scheme (Multirate (Synctime: 1 ms)) each physical domain is solved using again an explicit fixed step method. A step size of 1 μs is applied for the electric system. All other domains take 1 ms as step size. Again, an explicit fixed step method is used, whereby the chosen step size is now 1 μs. The information exchange takes place after each macrostep of 1 ms. The second multirate scheme (Multirate (Synctime: Multiple)) uses, again as in the previous example, multiple synchronization times. The electric system is still solved with 1 μs. But now, only the solid thermal circuits apply steps of 1 ms. The gas and fluid circuit use a step size of 5 ms. Therefore, within the system two synchronization times are used, 1 ms and 5 ms. Again, like in the previous example a simplified evaluation order can be given by:  O(1000)). The results of an IGBT show that the temperature is the same in all three constellations, cf. Fig. 7. In general, the results of the multirate solvers are sufficiently accurate in this example.
Again, significant speed up in the calculation time can be achieved, cf. Fig. 8. The multirate approaches are more than 6 times faster than the singlerate approach, cf. Table 2.
Comparing the two multirate solutions with each other, it can be seen that the use of different synchronization times does not provide any further significant improvement. This is because the thermal networks of fluid, gas and solid circuits require to be resolved on almost the same time scale compared to the electric system. In general, the RTFs for this example are very high as the electrical system must be resolved very accurately. Thus, despite enormous performance improvements, the simulation is still far from being realtime capable, which means that further improvements both on modelling side as well as the solver side are necessary.

Figure 7
Comparison of temperature in an IGBT of multirate cases against a single solver case for three-phase inverter with cooling system simulation

Conclusion
As shown, the multirate approach offers a possibility to reduce computation time considerably. Its can be adapted to the requirements of the model to be simulated. Any physical system can be accompanied by a solver tailored to its own needs. If the solver step sizes in the model vary strongly between the different solvers (n 2), a stable solution can still be guaranteed by a hierarchical synchronization approach. To ensure a stable simulation, automatic index reduction of the physical networks, appropriate solver settings for each subsystem and an adequate coupling procedure, are important. For the correct choice a significant speed up can be achieved, while satisfying the accuracy requirements. Despite the considerable improvements using the presented co-simulation approach, this still can be improved. Thus, there is the possibility to save additional computation time by parallelisation and intelligent inter-/extrapolation methods. However, further analysis of the considered systems and the applied methods are necessary to achieve this properly.