Modeling of Electriﬁed Transportation Systems Featuring Multiple Vehicles and Complex Power Supply Layout

: The paper proposes a novel approach to modeling electriﬁed transportation systems. The proposed solution reﬂects the mechanical dynamics of vehicles as well as the distribution and losses of electric supply. Moreover, energy conversion losses between the mechanical and electrical subsystems and their bilateral inﬂuences are included. Such a complete model makes it possible to replicate, e.g., the impact of voltage drops on vehicle acceleration or the necessity of partial disposal of regenerative braking energy due to temporary lack of power transmission capability. The modeling methodology uses a ﬂexible twin data-bus structure, which poses no limitation on the number of vehicles and enables modeling complex traction power supply structures. The proposed solution is suitable for various electriﬁed transportation systems including suburban and urban systems. The modeling methodology is applicable i.a. to Matlab/Simulink, which makes it broadly available and customizable, and provides short computation time. The applicability and accuracy of the method were veriﬁed by comparing simulation and measurement results on an exemplary trolleybus system operating in Pilsen, Czech Republic. Simulation of daily operation of an area including four supply sections and maximal simultaneous number of nine vehicles showed a good conformance with the measured data, with the difference in the total consumed energy not exceeding 5%.


Introduction
Electrified transportation systems such as trams, trolleybuses, or metros feature complex power supply structures consisting of a set of interconnected traction substations, feeders, and catenary sections. In urban systems, there is typically a large number of vehicles covering various routes and being in different drive modes at a particular time. This makes the distribution of electrical power in the supply system complex and variable. Due to the transmission losses, the voltages on vehicle current collectors change. Upon an excessive voltage drop, the vehicle control system needs to limit the collected current to assure stable operation of the whole system. In turn, during regenerative braking, the voltage on the vehicle collector rises and the level of this voltage is used to control the portion of braking power that cannot be transferred to other vehicles and needs to be dissipated in the braking resistor. Hence, energy flow in electrified transportation systems depends on the electric power supply structure and parameters, as well as on mechanical variables associated with driving mode and position of the vehicles. Moreover, these mechanical and electrical subsystems influence each other, and this influence may also be altered by the vehicle control systems. All these features make comprehensive modeling of electrified transportation systems very difficult. uniform velocity profiles, and they cannot include a certain degree of randomness required to model realistic road traffic [18,28].
This paper presents an approach to modeling electrified transportation systems with multiple vehicles and complex power supply layout, typical for trolleybus or tram systems. The novel aspects of this work include modeling every part of the transport network as a separate entity, allowing for unrestricted power supply layout (including branch lines), and timetable implementation. Similarly, the vehicles may be of different types and follow different routes, as each one of them has its own schedule. Furthermore, velocity profiles are derived from a large set of measurement data and are selected semi-randomly. High flexibility was achieved by the adoption of a modular model architecture and introduction of a twin data bus structure. The modeling methodology is explained using Matlab/Simulink (MathWorks, Inc., Natica, MA, USA) as exemplary software implementation. The applicability and accuracy of the method were verified by comparing simulation and measurement results for an example trolleybus system in Pilsen, Czech Republic.

Proposed Model Structure and Software Setup
Various transportation systems have different features. For instance, in railway systems, a power section is most often supplied by two traction substations located near the ends of the section (bilateral supply). In turn, urban systems use a single traction substation to supply the power section, but multiple feeders (cables) may be used to connect the substation to different points of the section. Moreover, the spatial layout of an urban power supply system is typically much more complex than that of a railway system. Railway vehicles are able to strictly follow the timetable and their speed profiles are very repetitive due to the lack of influence of traffic. In turn, trolleybuses run in congestion, which makes their speed profiles unrepeatable. The proposed modeling approach is feasible to include features of various electrified transportation systems. However, in order to keep the paper concise and clear, its further content refers to the trolleybus transportation system as an example. Modeling such a system is considered challenging due to the complex spatial and electrical structure of the power supply system, large number of vehicles running in particular power sections, and random speed profiles due to varying congestion.
The general structure of the proposed model is shown in Figure 1. It consists of a group of subsystems related to vehicles (green) and a group related to the electric power supply (blue). These groups exchange data throughout data buses that work similarly to industrial communication networks (ICNs) [29,30]. This paper presents an approach to modeling electrified transportation systems with multiple vehicles and complex power supply layout, typical for trolleybus or tram systems. The novel aspects of this work include modeling every part of the transport network as a separate entity, allowing for unrestricted power supply layout (including branch lines), and timetable implementation. Similarly, the vehicles may be of different types and follow different routes, as each one of them has its own schedule. Furthermore, velocity profiles are derived from a large set of measurement data and are selected semi-randomly. High flexibility was achieved by the adoption of a modular model architecture and introduction of a twin data bus structure. The modeling methodology is explained using Matlab/Simulink (MathWorks, Inc., Natica, MA, USA) as exemplary software implementation. The applicability and accuracy of the method were verified by comparing simulation and measurement results for an example trolleybus system in Pilsen, Czech Republic.

Proposed Model Structure and Software Setup
Various transportation systems have different features. For instance, in railway systems, a power section is most often supplied by two traction substations located near the ends of the section (bilateral supply). In turn, urban systems use a single traction substation to supply the power section, but multiple feeders (cables) may be used to connect the substation to different points of the section. Moreover, the spatial layout of an urban power supply system is typically much more complex than that of a railway system. Railway vehicles are able to strictly follow the timetable and their speed profiles are very repetitive due to the lack of influence of traffic. In turn, trolleybuses run in congestion, which makes their speed profiles unrepeatable. The proposed modeling approach is feasible to include features of various electrified transportation systems. However, in order to keep the paper concise and clear, its further content refers to the trolleybus transportation system as an example. Modeling such a system is considered challenging due to the complex spatial and electrical structure of the power supply system, large number of vehicles running in particular power sections, and random speed profiles due to varying congestion.
The general structure of the proposed model is shown in Figure 1. It consists of a group of subsystems related to vehicles (green) and a group related to the electric power supply (blue). These groups exchange data throughout data buses that work similarly to industrial communication networks (ICNs) [29,30].   The vehicle subsystems, acting as transmitters, send their output data vectors to the Vehicle Data Bus, which broadcasts these vectors to receivers. The receivers, which are power section subsystems, use selectors to extract the data of vehicles that run along a particular power section at a considered time. Hence, the selectors perform a task similar to that of frame ID-based filtering used in ICNs. Particular substations are modeled as independent subsystems that exchange data with power section subsystems. This provides the possibility to model various interconnections between the substations and the power sections, and the bilateral influence of their state variables. The power section subsystems send the results of computations to the Supply Data Bus, which broadcasts them to vehicle subsystems. Similarly to the Vehicle Data Bus and the power supply subsystems, the vehicle subsystems use deselectors to extract the results related to particular vehicles. The proposed modeling approach allows for comprehensive analysis of electrified transportation systems, without typical limitations from other solutions. High versatility has been achieved through making the model structure flexible and thanks to local data handling by multiple matrices instead of merging the data into the global matrix.
More details on the data exchange are shown in Figure 2, on the model of a single vehicle and single power section selected as an example. The vehicle subsystem models the vehicle run according to a pre-programmed velocity vs. time profile defined in the look-up table. In order to model the impact of congestion on the trolleybus service, the vehicle subsystem may consist of several velocity profiles (e.g., recorded on the considered route or designed based on the analysis of the route). In such a case, the profiles are selected randomly by a random number generator (RNG) seed assigned to each vehicle. The computations of motion dynamics provide a set of mechanical variables such as motive or braking force, velocity, and covered distance. These computations exchange data with the uploaded route profile, which allows for reflecting additional resistance forces related to riding up or downhill. The vehicle subsection also comprises a look-up table that links the covered distance with the identifier of corresponding power supply section and computes the position of the vehicle within this section.
The current and energy of the vehicle collector are calculated based on mechanical variables, with consideration of the variable efficiency factor and the current collector voltage. The computed energy is split into consumed, regenerated, and dissipated in the braking resistor. This, along with the drivetrain losses calculation, allows for the complete vehicle energy efficiency analysis. Additionally, a feedback path from electricalto-mechanical calculations is provided, in order to reflect the power limitation caused by excessively low voltage on the current collector. Each vehicle is represented by an independent subsystem, therefore the number of modeled vehicles may be freely adjusted. Moreover, particular vehicles may have different parameters and may run on different routes.
The vehicle subsystems send the output data vectors to the Vehicle Data Bus. Each of these vectors consists of vehicle identifier k veh , supply section identifier k sec , vehicle collector current i veh , and position of the vehicle in section x veh . The selectors isolate vectors corresponding to particular power sections and store them in local matrices. It is worth noting that the number of vehicles within the section may vary in time, so the local matrix must be declared with respect to the maximum number of vehicles (because variable-size arrays are inefficient or not supported, depending on the programming environment used), and filled with dummy vehicle data in case of a lower vehicle number. The vehicles can run on different routes, and they will not necessarily enter the section always in the same order. Hence, an important role of the selector is to sort the vehicle vectors in the local matrix in order of their location. This is required by the main algorithm of the power section subsystem that merges power section and vehicles data and builds the electric circuit equivalent. The model of power section synthesizes the electric circuit based on the parameters prepared by the selectors and on the spatial layout and parameters of the power section. The layout of the electric power supply does not change in time, so the data connections between power section subsystems and substation subsystems are fixed. The substation subsystem calculates the feeder output voltage that depends on the load current, internal resistance, and feeder resistance. Since most substations in DC systems are equipped with diode rectifiers, the reverse current is not allowed. Interacting with the linked substations, the power section subsystem calculates the catenary voltage for vehicle nodes (collector voltages) and compiles the data into the result matrices that are delivered to the Supply Data Bus. These matrices consist of vehicle identifiers kveh and current collector voltages uveh. The voltage subsystems use deselectors to extract the collector voltage of particular vehicles. This completes the loop of model computations.
The basic parameters of the analyzed transport network, as well as the initial conditions for the simulation, are loaded from external files. The results are obtained directly from the corresponding models. It is possible to save them into a file or display them as waveforms using commands from outside the simulation.
The models of vehicles, supply sections, and substations are independent-they can be extracted from the main model and used for individual analysis if relevant parameters are provided. The number of vehicles is limited only by the available computer resources, and not by the model design. The modeled vehicles may be of different type and follow individual schedules and routes. Similarly, there is virtually no limit for the number of power sections or connections between them. This allows for simulating complex transport networks because the sectors need only be connected with the data buses and the substation, while the connection between them is defined in the initial parameters file. The model of power section synthesizes the electric circuit based on the parameters prepared by the selectors and on the spatial layout and parameters of the power section. The layout of the electric power supply does not change in time, so the data connections between power section subsystems and substation subsystems are fixed. The substation subsystem calculates the feeder output voltage that depends on the load current, internal resistance, and feeder resistance. Since most substations in DC systems are equipped with diode rectifiers, the reverse current is not allowed. Interacting with the linked substations, the power section subsystem calculates the catenary voltage for vehicle nodes (collector voltages) and compiles the data into the result matrices that are delivered to the Supply Data Bus. These matrices consist of vehicle identifiers k veh and current collector voltages u veh . The voltage subsystems use deselectors to extract the collector voltage of particular vehicles. This completes the loop of model computations.
The basic parameters of the analyzed transport network, as well as the initial conditions for the simulation, are loaded from external files. The results are obtained directly from the corresponding models. It is possible to save them into a file or display them as waveforms using commands from outside the simulation.
The models of vehicles, supply sections, and substations are independent-they can be extracted from the main model and used for individual analysis if relevant parameters are provided. The number of vehicles is limited only by the available computer resources, and not by the model design. The modeled vehicles may be of different type and follow individual schedules and routes. Similarly, there is virtually no limit for the number of power sections or connections between them. This allows for simulating complex transport networks because the sectors need only be connected with the data buses and the substation, while the connection between them is defined in the initial parameters file. Since the power sections also have their identifiers and localization parameters, the implementation of branches can be achieved easily-the power section can be connected in the same manner as the vehicles (but with constant location).
Using Simulink as the programming environment has some limitations. The most notable one is the lack of dynamic matrix size. Moreover, the input and output ports of the program blocks have fixed width [31]. To deal with these limitations, countermeasures in the form of data buses, selectors, and deselectors were implemented. It is also worth noting that all variables in the matrices must be of uniform numeric type. The program is designed to use the ode3 solver and a fixed time step equal to 1 s, which provides sufficient accuracy while retaining high computing performance [32]. The selection of time step has been validated by comparing simulation results obtained for execution with the chosen and decreased time step (0.01 s). For the simulation of a 10-min trolleybus run, the difference in its consumed energy was less than 0.5%. The time-step size of 1 s has been commonly used also in similar models [3,6,8,18].

Vehicle Subsystem
Each modeled vehicle has its individual subsystem that includes vehicle parameters, timetable, route inclination profile, and reference speed profiles. The ability to follow the reference speed profile is constrained by several factors of both mechanical and electrical nature. The vehicle motion is modeled mostly using kinematic equations, but the vehicle subsystem also includes algorithms that compute losses of conversion between electric and mechanical power. Some control-related factors are also included, e.g., the power limitation due to excessive voltage drop at the current collector.

Reference Speed Profiles and Timetabling
The reference speed profiles must be prepared before the execution of the model and stored in the look-up tables that are uploaded to the vehicle subsystem. Trolleybuses run along public roads, so modeling their ride with simple velocity profiles that include a sequence of acceleration, cruising, coasting, and braking would not reliably reflect their complex operating conditions. Therefore, the reference velocity profiles for the analyzed case were prepared with the use of onboard vehicle data recorders. For each considered trolleybus line and each type of vehicle, numerous recordings were carried out to obtain a representative set of reference velocity profiles. The recordings include runs that took place at different times of the day, in order to reflect time-specific factors such as traffic, stop-atlights patterns, etc. This allows for aligning the recorded runs into groups representing morning rush, midday, afternoon rush, and evening operation. The recordings also include dwelling time, so analysis of their content allows for the determination of minimal and maximal times that the trolleybuses spent at particular stops.
The aim of the proposed model is to reflect both deterministic features of the transportation system, such as the operation according to the timetable, and stochastic features such as different run and dwelling times. This is achieved by a set of signals that control the time when the vehicle enters the modeled area of the transportation system and when it starts the run from the terminal (end stop). These signals are also used for the random selection of velocity profiles related to rides between stops (featured by different run times) as well as the random selection of dwelling times at stops.
Exemplary waveforms of vehicle control signals are presented in Figure 3. The permission signal, when active, allows for the run of a given vehicle. The rising slope of this signal determines the time when the vehicle enters the considered area of the transportation system or resumes the operation after a service pause at the end stop (terminal). Each run between stops is initiated by the triggering signal, whose impulses are counted by the run counter. The initial value of this counter must be set before model execution to define the initial position of the vehicle. The value of the run counter determines between which stops the trolleybus is running. Each run starts with selecting the reference velocity profile from the set related to the considered route, stored in the vehicle subsystem. This selection is based on the random number generator (RNG). The probability of selecting Energies 2021, 14, 8196 7 of 20 each profile recorded for a considered route between stops is not equal. The RNG uses global simulation time to set weights for selecting profiles related to the different parts of the day. The highest probability is attributed to the group of velocity profiles that were recorded in the time similar to the actual global simulation time. The RNG is also used to set the dwelling time at the next stop (not shown in Figure 3). In this process, the range of selected numbers is defined by minimal and maximal time derived based on the recordings. 2021, 14, x FOR PEER REVIEW 7 of 20 from the set related to the considered route, stored in the vehicle subsystem. This selection is based on the random number generator (RNG). The probability of selecting each profile recorded for a considered route between stops is not equal. The RNG uses global simulation time to set weights for selecting profiles related to the different parts of the day. The highest probability is attributed to the group of velocity profiles that were recorded in the time similar to the actual global simulation time. The RNG is also used to set the dwelling time at the next stop (not shown in Figure 3). In this process, the range of selected numbers is defined by minimal and maximal time derived based on the recordings. The reference velocity profiles feature slightly different times of travelling between stops. As these times, as well as dwelling times at stops, are selected randomly, the travelling time between terminals (end stops) may vary. Hence, the maximal time for travelling between the terminals must be assumed, and if the vehicle reaches the end stop sooner than this time expires, the time that it spends at this stop is prolonged according to the recorded difference. This travelling time is reflected by the direction change signal, whose rising slope represents the latest assumed arrival at the end stop, while the falling slope determines the departure. The rising slope in the direction change signal also reverses the value of the direction signal, which can be set to −1 or 1. The value of the direction signal is multiplied by velocity in the algorithm computing the vehicle motion dynamics. This, in turn, influences the computations of the covered distance, which increases or decreases depending on the direction of movement.

Vehicle Motion Dynamics
The vehicle can be considered as a material point (vehicle mass is concentrated at a single point), which greatly reduces the complexity of the model. This assumption provides satisfactory modeling accuracy for most types of vehicles, except for freight trains with a large number of wagons, where the mass is distributed along a distance of hundreds of meters [33].
The equation for vehicle dynamics is given by: where a is the vehicle acceleration, F is the total motive force, W is the resistance force, m is the vehicle mass, and k is the coefficient of rotating mass. The reference velocity profiles feature slightly different times of travelling between stops. As these times, as well as dwelling times at stops, are selected randomly, the travelling time between terminals (end stops) may vary. Hence, the maximal time for travelling between the terminals must be assumed, and if the vehicle reaches the end stop sooner than this time expires, the time that it spends at this stop is prolonged according to the recorded difference. This travelling time is reflected by the direction change signal, whose rising slope represents the latest assumed arrival at the end stop, while the falling slope determines the departure. The rising slope in the direction change signal also reverses the value of the direction signal, which can be set to −1 or 1. The value of the direction signal is multiplied by velocity in the algorithm computing the vehicle motion dynamics. This, in turn, influences the computations of the covered distance, which increases or decreases depending on the direction of movement.

Vehicle Motion Dynamics
The vehicle can be considered as a material point (vehicle mass is concentrated at a single point), which greatly reduces the complexity of the model. This assumption provides satisfactory modeling accuracy for most types of vehicles, except for freight trains with a large number of wagons, where the mass is distributed along a distance of hundreds of meters [33].
The equation for vehicle dynamics is given by: where a is the vehicle acceleration, F is the total motive force, W is the resistance force, m is the vehicle mass, and k is the coefficient of rotating mass. By integrating the result of (1), the model computes the vehicle speed. Then, double integration provides the covered distance. In (1), the vehicle mass m is increased artificially by multiplying it by the coefficient k, whose value is slightly greater than 1 (for trolleybuses: k = 1.2 [34]). This mass increase reflects the impact of the moment of inertia of drivetrain rotating parts on the vehicle acceleration.
The total motive force F is either propelling (positive) or braking (negative), and its value is controlled by the driver. The positive motive force is generated solely by the electric drive. Hence, its value is limited by the nominal power and torque of the drive, see further comments in Section 3.3. The negative motive force that occurs during braking may be generated by both electric drive and friction brakes. The distribution of the braking force (set by the driver) between the electric drive and the friction brakes is performed by the drivetrain control algorithm.
The motion resistance W is the sum of resistance forces that can be attributed to the vehicle (fundamental resistance W f ) and the route (additional resistance W i ). The fundamental motion resistance of the vehicle is calculated as the sum of aerodynamic drag, transmission drag, and rolling resistance. The following second-order polynomial is used to approximate the value of the fundamental resistance force: where v-the vehicle velocity; A, B, C-the coefficients of the fitting function (e.g., for Škoda 26Tr trolleybus: The additional motion resistance W i , related to route inclination, is calculated as: where g is the gravitational acceleration and i is the longitudinal inclination of the route. The route profile that defines the route inclination as a function of covered distance can be derived from GPS recordings or extracted from a map supported by Shuttle Radar Topography Mission (SRTM)-based altitude data.

Control of Motive Force
The vehicle subsystem is designed to regulate the motive force F in such a way as to follow the reference velocity profile. This is carried out by the use of the proportionalintegral (PI) controller that sets the value of total motive force F according to the error computed as the difference between reference and actual speed of the vehicle (Figure 4). The gains for the proportional and integral terms were tuned to ensure control dynamics similar to those observed during the experiments. Since the controller output depends on the velocity error and its gains, the resulting force value can exceed the maximum motive force available in the traction drive. This cause of possible errors was solved by saturating the PI controller output (value of motive force) with respect to factors described below. By integrating the result of (1), the model computes the vehicle speed. Then, double integration provides the covered distance. In (1), the vehicle mass m is increased artificially by multiplying it by the coefficient k, whose value is slightly greater than 1 (for trolleybuses: k = 1.2 [34]). This mass increase reflects the impact of the moment of inertia of drivetrain rotating parts on the vehicle acceleration.
The total motive force F is either propelling (positive) or braking (negative), and its value is controlled by the driver. The positive motive force is generated solely by the electric drive. Hence, its value is limited by the nominal power and torque of the drive, see further comments in Section 3.3. The negative motive force that occurs during braking may be generated by both electric drive and friction brakes. The distribution of the braking force (set by the driver) between the electric drive and the friction brakes is performed by the drivetrain control algorithm.
The motion resistance W is the sum of resistance forces that can be attributed to the vehicle (fundamental resistance Wf) and the route (additional resistance Wi). The fundamental motion resistance of the vehicle is calculated as the sum of aerodynamic drag, transmission drag, and rolling resistance. The following second-order polynomial is used to approximate the value of the fundamental resistance force: where v-the vehicle velocity; A, B, C-the coefficients of the fitting function (e.g., for Škoda 26Tr trolleybus: The additional motion resistance Wi, related to route inclination, is calculated as: where g is the gravitational acceleration and i is the longitudinal inclination of the route. The route profile that defines the route inclination as a function of covered distance can be derived from GPS recordings or extracted from a map supported by Shuttle Radar Topography Mission (SRTM)-based altitude data.

Control of Motive Force
The vehicle subsystem is designed to regulate the motive force F in such a way as to follow the reference velocity profile. This is carried out by the use of the proportionalintegral (PI) controller that sets the value of total motive force F according to the error computed as the difference between reference and actual speed of the vehicle (Figure 4). The gains for the proportional and integral terms were tuned to ensure control dynamics similar to those observed during the experiments. Since the controller output depends on the velocity error and its gains, the resulting force value can exceed the maximum motive force available in the traction drive. This cause of possible errors was solved by saturating the PI controller output (value of motive force) with respect to factors described below.  The maximal output torque of an electric drive (hence also the motive force) is most commonly defined as a function of velocity (see Figure 5). Upon low velocities, the maximal Energies 2021, 14, 8196 9 of 20 output torque is constant. However, after reaching the base speed at which the nominal value of output power is achieved, the torque is limited according to the constant (nominal) power curve. In certain operating conditions, the available motive force may be further decreased. For instance, if the voltage on the trolleybus current collector drops below a predefined value, the output power of the electric drive is decreased linearly with respect to the catenary voltage in order to prevent further voltage drop and to sustain the stability of the supply system. Such emergency power limitation was implemented into the model [35]. When modeling some transportation systems (e.g., railway), it is sometimes desired to reflect limited adhesion of driving wheels [36]. However, in the trolleybus system, it is assumed that wheel slip or skid do not occur. The above-mentioned limitations of the propelling force generated by the electric drive were implemented as the upper saturation threshold for PI controller output (see Figure 4). cordings were carried out, heating and air conditioning were not used.
The efficiency of the mechanical part of the drivetrain was set at 95%. It was assumed that there is no need to model the variable efficiency, as the operating conditions related to significant efficiency drops (very low torque or speed) do not occur often during trolleybus runs [37].
Computing power losses Pem_loss in the electric motor and Pinv_loss in the power electronic inverter is a complex problem. As the operating conditions of vehicle drives change in a wide range, assuming a constant efficiency would lead to considerable errors [38]. Trolleybuses in the considered transportation system use induction motor drive. The induction motor converts electric to mechanical power using the electromagnetic field. The efficiency of the induction motor is highly dependent on angular velocity and torque. While for most operating conditions the efficiency exceeds 90%, there are also low-efficiency regions corresponding to the low-speed and low-load operation. In order to accurately reflect motor power losses in the model, the efficiency map η = f(T,ω) is required. Since experimentally derived efficiency values were unavailable, the authors used parameters of the electric motor to set up a mathematical model, which was used to compute losses as a function of velocity and torque. The obtained dataset was loaded into the vehicle subsystem as a 2-dimensional lookup table, where the efficiency values for arguments between the defined points are computed using linear interpolation. A graphical representation of the derived efficiency map is shown in Figure 5. A detailed analysis of the power electronic inverter would require simulating transient states related to the switching of its transistors. This, in turn, would lead to setting simulation time step at the order of microseconds, which is unjustified for other parts of the considered model. Therefore, the inverter losses are estimated based on the actual operating conditions and parameters of power electronic components specified in their datasheets. The total inverter losses Pinv_loss are computed as the sum of switching losses of transistors Pt_sw, reverse-recovery diode losses Pd_rec, and conduction losses for both diodes Pd_cond and transistors Pt_cond [39,40]. All these losses need to be included because the typical voltage drop of 1.2-2.5 V on a semiconductor switch translates to kilowatts of losses. The negative motive force, which corresponds to braking, was assumed as a blended force generated by friction brakes and electric drive. The latter part used for braking is commonly referred to as an electrodynamic brake. The available braking force related to blended braking has large values, allowing even for emergency braking, but this is assumed not to occur in the simulated case. However, for formal reasons, the maximal braking force was defined and used to set the lower saturation threshold for the PI controller output (see Figure 4). This value is negative as it corresponds to braking.
To compute the vehicle dynamics (1) for the braking regime, the total motive force F generated by both electrodynamic and friction brakes should be considered. However, the vehicle subsystem also includes computations related to the electric drive. For these computations, it is necessary to isolate the braking force F em generated solely by the electric motor. In most blended braking systems, the vehicle controller maximizes the use of electrodynamic brake to save friction brake pads and allow for generating the electric power that can be delivered to other vehicles. Therefore, it was assumed that friction brakes are engaged only if the total braking force exceeds the value defined by the electric drive relation between the motive force and velocity ( Figure 5). Consequently, the motive force generated by the electric drive is computed as a total motive force with a lower saturation threshold set accordingly to the inverted (negative) maximal electric drive motive force at actual speed.

Electric Drive Efficiency and Electric Current Computation
Electric current I veh collected by the vehicle is determined from mechanical power on electric drive output but also includes other important factors. The electric current is computed from the following formula: (F em ·v) + P br_res + P aux + P mech_loss + P em_loss + P inv_loss U veh , where (F em ·v) is the mechanical power on electric motor output (may be either positive or negative); P br_res is the power dissipated in the braking resistor (positive or zero); P aux is the power of auxiliary loads (positive); P mech_loss , P em_loss , P inv_loss , are the power losses in electric motor, power inverter, and the mechanical part of the drivetrain, respectively (positive or zero); and U veh is the voltage on vehicle pantograph. In order to exclude the operation of friction brakes, which has no impact on the vehicle current, the mechanical power computed using Equation (4) is not based on total motive force F but the part of this force F em generated by the electric motor.
The power balance in the numerator of (4) includes power P br_res dissipated into the heat in the vehicle braking resistor. This power is equal to zero during propelling and coasting. On the other hand, P br_res may be positive during braking if other vehicles on the same power supply section are unable to consume the electric power generated by regenerative braking in the considered vehicle. The power P br_res is controlled using the vehicle pantograph voltage U veh and voltage threshold levels. Typically, the braking resistor is engaged when the pantograph voltage exceeds 760 V. The power P br_res increases linearly with the further rise of U veh , reaching the full value of vehicle power balance at U veh = 790 V. Conversely, when the vehicle pantograph voltage falls below 580 V, the vehicle power is reduced linearly down to zero at U veh = 520 V, which is the minimum allowed voltage for traction drive of the analyzed vehicle.
The auxiliary power was set at a constant value, P aux = 2.5 kW, derived from onboard recordings carried out in trolleybuses during their standstill. In the period when the recordings were carried out, heating and air conditioning were not used.
The efficiency of the mechanical part of the drivetrain was set at 95%. It was assumed that there is no need to model the variable efficiency, as the operating conditions related to significant efficiency drops (very low torque or speed) do not occur often during trolleybus runs [37].
Computing power losses P em_loss in the electric motor and P inv_loss in the power electronic inverter is a complex problem. As the operating conditions of vehicle drives change in a wide range, assuming a constant efficiency would lead to considerable errors [38]. Trolleybuses in the considered transportation system use induction motor drive. The induction motor converts electric to mechanical power using the electromagnetic field. The efficiency of the induction motor is highly dependent on angular velocity and torque. While for most operating conditions the efficiency exceeds 90%, there are also low-efficiency regions corresponding to the low-speed and low-load operation. In order to accurately reflect motor power losses in the model, the efficiency map η = f(T,ω) is required. Since experimentally derived efficiency values were unavailable, the authors used parameters of the electric motor to set up a mathematical model, which was used to compute losses as a function of velocity and torque. The obtained dataset was loaded into the vehicle subsystem as a 2-dimensional lookup table, where the efficiency values for arguments between the defined points are computed using linear interpolation. A graphical representation of the derived efficiency map is shown in Figure 5.
A detailed analysis of the power electronic inverter would require simulating transient states related to the switching of its transistors. This, in turn, would lead to setting simulation time step at the order of microseconds, which is unjustified for other parts of the considered model. Therefore, the inverter losses are estimated based on the actual operating conditions and parameters of power electronic components specified in their datasheets. The total inverter losses P inv_loss are computed as the sum of switching losses of transistors P t_sw , reverse-recovery diode losses P d_rec , and conduction losses for both diodes P d_cond and transistors P t_cond [39,40]. All these losses need to be included because the typical voltage drop of 1.2-2.5 V on a semiconductor switch translates to kilowatts of losses.
A typical traction power inverter consists of 6 transistors and 6 freewheeling diodes, so the total power losses in the inverter can be approximated as: For a single transistor, the switching losses are approximated as [39,40]: where U veh_n -the nominal voltage at vehicle collector (on inverter DC bus); Es w_on and E sw_off -energy losses for switching the transistor on and off, respectively; f -the switching frequency.
The reverse recovery losses of a diode are calculated as [39]: where E rec -energy dissipated due to the reverse recovery charge current. The conduction losses for a transistor and a freewheeling diode are computed according to (8) and (9), respectively [38,39]: where I C -collector current; U CE -collector-emitter saturation voltage; D-inverter duty cycle; cosϕ-power factor of the motor; and U f -diode forward voltage drop. The collector current I C of the transistor and diode was assumed equal to the vehicle motor current I veh -averaged over the duty cycle. The motor power factor was implemented in the second term of Equations (8) and (9).

Power Supply Subsystem
The power supply of electrified transportation systems may have a complex structure, especially in the case of an urban system. The trolleybus supply system consists of traction substations that transform the electric power from medium-voltage AC to DC voltage at a level of 600-700 V, which is suitable for vehicles. Typically, a trolleybus traction substation has a rated power of a few hundred kilowatts and supplies the catenary at the area of several square kilometers. The catenary is divided into power sections of a typical length of a few hundred meters. The substation supplies several neighboring power sections through separate feeders (supplying cables). The neighboring power sections of the catenary are electrically separated from each other by section insulators fixed to the catenary. However, they are electrically linked by the feeders if they are supplied from the same substation. Some power sections may include intersections and then the catenary has some side branches. Unlike railways, power sections of urban systems are supplied from a single traction substation. However, several feeders may connect the substation to different points of the power section.
The power supply subsystem includes models of power sections and traction substations. Such a model composition allows for flexible representation of complex supply structures. Moreover, such a division splits the process of numerical solving into smaller portions, which is favorable in terms of computation performance and stability. Following the above-discussed features of the trolleybus supply system, typically a single substation model is linked to models of several associated power sections.

Traction Substation and Feeders
The model of substation and feeders computes the output voltages for the feeders. The results are delivered to the linked power section model. The traction substation is modeled as voltage source E and internal series resistance R s . Such a model reflects the voltage drop at substation output (at its DC switchgear) that is linear to the output current.
The best method to derive the model parameters is to analyze the recordings collected during the operation of the considered substation. Approximating the relation between the substation output voltage U sub and output current I sub with a linear function allows for deriving both the idle voltage E (as the voltage for I sub = 0) and the series resistance (as the slope coefficient of the linear function).
The traction substation model is supplemented by the models of all feeders associated with this substation. The feeders are modeled as series resistances R fdr whose values depend on the feeder length, cross-section, and the restiveness of the material used (typically aluminum). The trolleybus catenary is symmetrical, i.e., the feeders and the catenary contact wires are identical for positive and negative poles. For simplicity, each feeder is modeled as a single resistance R fdr , whose value is the sum of resistances of both negative and positive pole feeders. The feeders of the substation have different lengths and different load currents, so the output voltage for each feeder needs to be computed individually. The output voltage for the k-th feeder in the substation supplying n feeders is computed as follows: The resistances of the feeders are also used for computing the power transmission losses in the feeders: The substations in the DC supply systems are usually equipped with diode rectifiers, so they do not allow reverse current. To reflect this property, the value of the modeled internal resistance R s is a function of the substation output current I sub . The resistance R s increases when the output current approaches zero and saturates at 1 GΩ for negative currents. Such a continuous change of resistance synergizes well with the voltage-dependent reduction of regenerative braking power implemented in the vehicle model and ensures stable operation of the numeric solver.

Power Section
The power section model uses a circuit equivalent of the catenary to compute feeder currents and voltages at vehicle pantographs. The feeders are modeled as voltage sources with values obtained from the model of substation and feeders. The vehicles are modeled as current sources, whose values are set according to the vehicle data obtained from the Vehicle Data Bus through the selector (see Figure 1). The catenary is modeled as series resistance R cat . The catenary is assumed to have a unified value of unit resistance R cat expressed in Ω/m. Consequently, the resistances between the mentioned current and voltage sources are computed based on their distances. The layout of the catenary is fixed, which also refers to the positions where the feeders are connected to the catenary. The variable topology of the equivalent circuit is related to the number and positions of vehicles that run within the considered power section. This problem is resolved by assuming a maximal number of vehicles in the equivalent circuit. If the actual number of vehicles is lower, the current sources related to absent vehicles are set to zero (while their set positions may be arbitrary). Consequently, the absent vehicles do not influence the current distribution in the circuit. This approach leads to a fixed structure of the equivalent circuit, where the values of voltage sources, current sources, and resistances are not stationary. Such a circuit can be solved with relatively simple and fast methods.
Another means that simplifies the equivalent circuit synthesis is sorting the vehicles related to the considered power section according to their position, which is done by the selector. Due to the sorting, it is easy to automatically set the catenary resistances and current source values according to the received data at each step of numeric computations.
Typical structures of power sections and their equivalent circuits are shown in Figure 6.
Another means that simplifies the equivalent circuit synthesis is sorting the vehicles related to the considered power section according to their position, which is done by the selector. Due to the sorting, it is easy to automatically set the catenary resistances and current source values according to the received data at each step of numeric computations.
Typical structures of power sections and their equivalent circuits are shown in Figure  6. The most common case of power supply section layout is shown in Figure 6a, where a single set of feeders is connected at the boundary of the power section. This case can be solved using self-reliant mathematical formulas for each derived parameter. A much less frequent case, where the power section is supplied by two sets of feeders is presented in Figure 6b. In this case, both feeders are supplied from the same substation. However, due to their different length and load, the output voltages of the feeders (represented by voltage sources Ufdr1 and Ufdr2) are different. This case can be solved automatically, e.g., by nodal analysis, which leads to solving the following set of equations: The most common case of power supply section layout is shown in Figure 6a, where a single set of feeders is connected at the boundary of the power section. This case can be solved using self-reliant mathematical formulas for each derived parameter. A much less frequent case, where the power section is supplied by two sets of feeders is presented in Figure 6b. In this case, both feeders are supplied from the same substation. However, due to their different length and load, the output voltages of the feeders (represented by voltage sources U fdr1 and U fdr2 ) are different. This case can be solved automatically, e.g., by nodal analysis, which leads to solving the following set of equations: If there are two or more sets of feeders connected to the power section, but not at the boundaries (the case not shown in Figure 6), such a section can be modeled as a combination of several self-reliant circuits from the set shown in Figure 6a,b. Figure 6c shows a more complex case, where a part of supply sections splits into two branches. The side branch of the catenary may be considered as a subcircuit represented by a current source whose value is equal to the total current of the vehicles running at the branch.
Since all vehicle models communicate with the power supply model using the Vehicle Data Bus (see Figure 1), there is no explicit limitation to the system layout. When part of the power section branches off, this case can be implemented as a separate stationary section in the place where it is connected to the main line.
The output values calculated by the power section model are loaded into the data vectors along with the corresponding vehicle ID-the order of results is consistent with that of vehicles. This allows omitting the use of search functions or multidimensional matrices, thus improving the simulation performance. The data vectors are then sent to the Supply Data Bus, where they are awaited by the deselectors, which send the line voltage value to individual vehicles.
Aside from the above-discussed results, the power section model computes the power transfer losses related to the catenary. This is done in a similar manner to computing power transfer losses in feeders (11).

Experimental Verification
The proposed model was verified within the framework of the Interreg project "Effi-cienCE", where the improvement of energy efficiency in urban electrified transportation is the main objective. The model was parametrized according to a part of the trolleybus network in the City of Pilsen, Czech Republic. The Pilsen supply system for electrified transportation systems includes 10 traction substations that supply both trolleybuses and trams. The verification was carried out for the supply area of substation MR5 (Figure 7), which was selected because all vehicles running in this area are equipped with AC-motor drive that enables regenerative braking. The verification covered a single workday operation (from 4 a.m. to midnight). The substation MR5 powers four trolleybus supply sections (marked from "A" to "D" in Figure 7) and three tram supply sections (not included in Figure 7). Two of the trolleybus sections are supplied by feeders connected at the boundaries, i.e., sections B and C. These sections are represented by the equivalent circuit shown in Figure 6a. Section D is supplied by two sets of feeders marked as 52A and 52B, and these sets are connected to the catenary at a 230-m distance. This section is represented by a set of equivalent circuits composed of one circuit shown in Figure 6b and two circuits from Figure 6a. The south section, marked as A, is powered by twin feeders connected to the catenary at the same point. Hence, they may be modeled as a single set of feeders with a double crosssection. However, this section includes a branch that leads to the Nova Hospodá loop, which is the end stop for trolleybuses. This branch is represented in the model as the equivalent circuit shown in Figure 6c.
The  The substation MR5 supplies both trolleybus and tram supply sections. The verification was focused on the trolleybus system, but the trams have a notable influence on the operation of the substation and the trolleybus supply section voltages. Therefore, in order to include this influence, the load caused to the substation by the trams was recorded during a selected mid-week day and this recording was used to set current waveforms of the tram feeders in the model. The trolleybus supply sections powered from substation MR5 are shown in Figure 7.
The substation MR5 powers four trolleybus supply sections (marked from "A" to "D" in Figure 7) and three tram supply sections (not included in Figure 7). Two of the trolleybus sections are supplied by feeders connected at the boundaries, i.e., sections B and C. These sections are represented by the equivalent circuit shown in Figure 6a. Section D is supplied by two sets of feeders marked as 52A and 52B, and these sets are connected to the catenary at a 230-m distance. This section is represented by a set of equivalent circuits composed of one circuit shown in Figure 6b and two circuits from Figure 6a. The south section, marked as A, is powered by twin feeders connected to the catenary at the same point. Hence, they may be modeled as a single set of feeders with a double cross-section. However, this section includes a branch that leads to the Nova Hospodá loop, which is the end stop for trolleybuses. This branch is represented in the model as the equivalent circuit shown in Figure 6c.
The area lies within the bounds of trolleybus lines 12, 15, and 17 (as well as tram line 2). The locations of the trolleybus stops and the route lengths are shown in Figure 8.  The substation MR5 powers four trolleybus supply sections (marked from "A" to "D" in Figure 7) and three tram supply sections (not included in Figure 7). Two of the trolleybus sections are supplied by feeders connected at the boundaries, i.e., sections B and C. These sections are represented by the equivalent circuit shown in Figure 6a. Section D is supplied by two sets of feeders marked as 52A and 52B, and these sets are connected to the catenary at a 230-m distance. This section is represented by a set of equivalent circuits composed of one circuit shown in Figure 6b and two circuits from Figure 6a. The south section, marked as A, is powered by twin feeders connected to the catenary at the same point. Hence, they may be modeled as a single set of feeders with a double crosssection. However, this section includes a branch that leads to the Nova Hospodá loop, which is the end stop for trolleybuses. This branch is represented in the model as the equivalent circuit shown in Figure 6c.
The area lies within the bounds of trolleybus lines 12, 15, and 17 (as well as tram line 2). The locations of the trolleybus stops and the route lengths are shown in Figure 8   Lines 12 and 15 are operated by standard two-axle trolleybuses type 26Tr. Line 17 is serviced by three-axle articulated vehicles type 27Tr. The parameters of these vehicles were provided by the manufacturer, Škoda Electric. The most important parameters used for the simulation are specified in Table 1. The model was set up according to a workday timetable. The schedule for line 12 begins with a 15-min service interval between subsequent runs, which is decreased to 5 min during the morning peak hours, i.e., between 5 a.m. and 8 a.m. Then, between 8 a.m. and 1 p.m., this interval is increased to 10 min. From 1 p.m. to 3 p.m., there is another peak with a 5-min interval, which from 5 p.m. starts increasing gradually to 10 min. The service interval is further increased to 15 min in the evening, starting from 7 p.m. Trolleybuses on line 15 operate with a 30-min service interval during the whole day, starting from 5:12 a.m. Vehicles on line 17 run along the analyzed area regularly between 5:20 and 6:00 a.m., with a 5-min interval. During the rest of the day, trolleybuses on line 17 appear at the analyzed area rarely, during only a few services. The maximal number of vehicles running simultaneously in each supply section is 6 for section A and 4 for each of the sections B, C, and D. These values were used to set the number of vehicle nodes in the equivalent electrical circuits included in the power supply section models.
The model was verified by comparing the simulation outcomes with real measurements from substation MR5 in Pilsen. The time of the measurements was the same as in simulation, i.e., workday from 4 a.m. until midnight. The collected measurements cover the voltage at the substation output bars and the currents of the substation feeders.
The currents and voltages, both in the measurements and simulation, were postprocessed to obtain the following quantities that were used for verification: • substation total output current: sum of feeder currents in substation MR5; • substation output voltage: voltage at substation MR5 output bars; and • substation output energy: integral of the total output power of substation MR5.
The model executes according to the timetable. As, in fact, the real operation times of the trolleybuses differ slightly from the scheduled time, the instantaneous voltage, current, and power values may differ between the model and the measurements. In order to disregard these differences, the waveforms of substation output voltage and current were averaged using a standard 15-min averaging window. The averaged variables still reflect daily load changes due to varying frequency of services well. The output energy does not require averaging as its computation involves integration.
The results of the verification are shown in Figures 9-11. Figure 9 compares the 15-min averages of substation output current. The results show good conformance between the measurement and simulation, considering the variability of the trolleybus network. The greatest current values were observed during morning peak hours when the service frequency was the highest. The current differences between 5 a.m. and 7:30 a.m. can be attributed to the differences between the assumed and the actual heating power and vehicle mass, which was set constant in the simulation. During the evening hours, the simulated current values slightly exceeded the measured ones, most likely due to traffic conditions and differences in the number of carried passengers. Figure 10 shows the 15-min averages of substation output voltage. The voltage values were very similar for both the measurement and simulation, with only small differences reflecting the deviations in substation output current that were commented in the previous paragraph. The 15-min averages of current ( Figure 9) showed good agreement between the measurement and simulation, and the 15-min averages of voltage ( Figure 10) were almost exactly equal. This confirms the good accuracy of the proposed model and therefore, the good feasibility of the simulation results. hicle mass, which was set constant in the simulation. During the evening hours, the simulated current values slightly exceeded the measured ones, most likely due to traffic conditions and differences in the number of carried passengers.  Figure 10 shows the 15-min averages of substation output voltage. The voltage values were very similar for both the measurement and simulation, with only small differences reflecting the deviations in substation output current that were commented in the previous paragraph. The 15-min averages of current ( Figure 9) showed good agreement between the measurement and simulation, and the 15-min averages of voltage ( Figure 10) were almost exactly equal. This confirms the good accuracy of the proposed model and were almost exactly equal. This confirms the good accuracy of the proposed model and therefore, the good feasibility of the simulation results.
In order to quantify the differences between measured and simulated waveforms of substation current and voltage, Root Mean Square Errors (RMSEs) were computed. The RMSE for the output voltage was 8.8 V, which constitutes less than 2% of the idle-state voltage. The RMSE for the output current was 49.2 A, i.e., 7% of the maximal measured current.  Figure 11 compares the output energy waveforms. Both waveforms are rising monotonously, as no energy can be delivered back to the substation. The end value of the energy obtained in the simulation was relatively close to that computed from the measurements. The notable and increasing difference between the waveforms was observed in the evening. The obtained results are consistent with the values of current and voltage. The observed differences stem from the assumed input parameters and the set schedule. As the deviation was less than 4.5%, the differences are within the expected margin. The performed verification has shown that good conformance between the measurement and the simulation was achieved. Slight differences can be attributed to unavailable parameters, which could not be set accurately in the model. For instance, the number of passengers was set constant in the simulation, whereas the real number varies in different services. Moreover, accurate modeling of auxiliary loads, such as heating and air conditioning, would require an additional thermal analysis [41] and precise knowledge about the weather and vehicle climate control settings. This knowledge was not available. Finally, the variability of traffic congestion was dealt with in the simulation by randomly selecting the recorded velocity profiles. However, the traffic conditions on particular days may differ due to specific occurrences.

Conclusions
The proposed model includes both mechanical and electrical subsystems of a transportation system, which makes it a comprehensive and versatile tool for various analyses. Adopting a structure similar to the industrial communication networks enables local data processing within the subsystems and overcomes common restrictions found in similar models. The verification carried out for an exemplary trolleybus system showed good accuracy of the model. Both the error in total energy consumption (less than 5%) and the RMSEs of substation current and voltage (7% and 2%, respectively) are considered satisfactorily small, taking into account numerous factors that could not be precisely reflected. The errors would have probably been even lower if the variable number of passengers and the variable auxiliary power were considered in the model. The model is ready for such a parametrization, but the mentioned factors were unknown in the analyzed case.
The model may be used for analyzing the behavior of both vehicles and power supply working under realistic operating conditions. The implementation of RNG-based velocity profile selection improves the accuracy of modeling variable traffic conditions. The use of Matlab as the programming basis allows for model interaction with numerous toolboxes and scripts. This, in turn, enables multi-objective optimization that is not limited In order to quantify the differences between measured and simulated waveforms of substation current and voltage, Root Mean Square Errors (RMSEs) were computed. The RMSE for the output voltage was 8.8 V, which constitutes less than 2% of the idle-state voltage. The RMSE for the output current was 49.2 A, i.e., 7% of the maximal measured current. Figure 11 compares the output energy waveforms. Both waveforms are rising monotonously, as no energy can be delivered back to the substation. The end value of the energy obtained in the simulation was relatively close to that computed from the measurements. The notable and increasing difference between the waveforms was observed in the evening. The obtained results are consistent with the values of current and voltage. The observed differences stem from the assumed input parameters and the set schedule. As the deviation was less than 4.5%, the differences are within the expected margin.
The performed verification has shown that good conformance between the measurement and the simulation was achieved. Slight differences can be attributed to unavailable parameters, which could not be set accurately in the model. For instance, the number of passengers was set constant in the simulation, whereas the real number varies in different services. Moreover, accurate modeling of auxiliary loads, such as heating and air conditioning, would require an additional thermal analysis [41] and precise knowledge about the weather and vehicle climate control settings. This knowledge was not available. Finally, the variability of traffic congestion was dealt with in the simulation by randomly selecting the recorded velocity profiles. However, the traffic conditions on particular days may differ due to specific occurrences.

Conclusions
The proposed model includes both mechanical and electrical subsystems of a transportation system, which makes it a comprehensive and versatile tool for various analyses. Adopting a structure similar to the industrial communication networks enables local data processing within the subsystems and overcomes common restrictions found in similar models. The verification carried out for an exemplary trolleybus system showed good accuracy of the model. Both the error in total energy consumption (less than 5%) and the RMSEs of substation current and voltage (7% and 2%, respectively) are considered satisfactorily small, taking into account numerous factors that could not be precisely reflected. The errors would have probably been even lower if the variable number of passengers and the variable auxiliary power were considered in the model. The model is ready for such a parametrization, but the mentioned factors were unknown in the analyzed case.
The model may be used for analyzing the behavior of both vehicles and power supply working under realistic operating conditions. The implementation of RNG-based velocity profile selection improves the accuracy of modeling variable traffic conditions. The use of Matlab as the programming basis allows for model interaction with numerous toolboxes and scripts. This, in turn, enables multi-objective optimization that is not limited by any program code and is compatible with other works carried out in this environment.
In future work, the authors are going to use this model to seek the optimal idle voltage of substations. Increasing the idle voltage minimizes the transmission losses but decreases the utilization of regenerative braking energy. The presented model will help to find the optimal value of the idle voltage in terms of total energy consumption.