Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

Modelling a Viennese ballroom: agent-based simulation to investigate complex behaviour

Modelling a Viennese ballroom: agent-based simulation to investigate complex behaviour MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 2020, VOL. 26, NO. 2, 169–192 https://doi.org/10.1080/13873954.2020.1727930 ARTICLE Modelling a Viennese ballroom: agent-based simulation to investigate complex behaviour a,b c c M. Bicher , S. Winkler and A. Körner a b Institute of Information Systems Engineering, TU Wien, Vienna, Austria; dwh Simulation Services, dwh GmbH, Vienna, Austria; Institute of Analysis and Scientific Computing, TU Wien, Vienna, Austria ABSTRACT ARTICLE HISTORY Received 4 March 2019 Dancing Viennese Waltz in one of the great historic ballrooms is an Accepted 6 February 2020 important and indispensable part of Austrian culture. This dance, while being tradition, is quite difficult to perform, especially if the KEYWORDS dance-floor is crowded. There, it is additionally challenging to avoid Agent-based simulation; collisions with other dancers, as they pace through the ballroom at ballroom dancing; complex a high velocity. Dependent on the dancer’s skill level, spinning speed system; Viennese Waltz can be adjusted to succeed. This paper presents an agent-based waltz model which makes it possible to investigate the influence of heterogeneously skilled dancers on the movement smoothness of the dancing crowd. Herein, each agent represents one dancing cou- ple in reality and it moves on the dance-floor by a rotatory motion with periodically switching rotation axes. Interaction between agents occurs via inelastic collisions. By performing a couple of case studies, we analyse and quantify the widespread rumour that the presence of only a few unskilled dancers disturbs the flow of the dancing crowd. 1. Introduction The waltz was first recorded in writing by Salomo Jakob Wolf in 1797 [1]. He described a novel dance that stands out by its dynamic spinning motion (old German walzen translates to to spin). This dance has evolved and became part of Austrian culture in the form of the famous Viennese waltz. It is part of every ball in Vienna and is traditionally danced at midnight on New Years Eve. On first glance, rhythm and step sequence seem fairly simple, but quick pace and spin make it more difficult than it looks. In contrast to other ballroom dances like Boogie or Cha-Cha-Cha, which are part of the Latin American dances, couples do not remain in a small dancing area, but explore the entire dance-floor. Moreover, in comparison to dances like Quickstep, Rhumba or Tango, dancing couples are much faster and con- stantly spinning. Hence, beside correct execution of the steps themselves, dancers need to be careful to avoid collisions with other dancers and furnishing. This is often challenging as typical ballrooms are usually quite crowded. All potential obstacles have to be identified and collision avoiding/handling strategies have to be made in a fraction of seconds while constantly spinning at more than 180 per second. Dizziness CONTACT M. Bicher martin.bicher@tuwien.ac.at Institute of Information Systems Engineering, TU Wien, Vienna, Austria © 2020 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any med- ium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way. 170 M. BICHER ET AL. is a typical side effect that makes the dance even more challenging [1,2]. This is the reason why the basic steps and simple movements of Viennese waltz can lead to complete chaos on the dance-floor. Collisions between dancers, bottlenecks and clusters of several couples are common phenomena. Yet, many effects are only understood on the scope of individual dancers so far. Detailed analysis of the behaviour of the crowd would require different research techniques, for example modelling and simulation. 1.1. Research objective In our study, we aim to model and simulate a crowded ballroom with waltz dancers to gain deeper insights and identify influential parameters that lead to clusters and bottle- necks on the dance-floor. The outcomes of the study should lead to general advises, that help dancers to avoid getting stuck on the dance-floor. In order to do so, we need to research a typical example of a complex system [3]: We observe the macro-behaviour of a system, the dancing crowd, which is the direct result of the micro-behaviour of many interacting elements, the dancers. Hereby, we need to model research problems on different scales: Micro-Scale: How does a dancing couple perform a technically correct waltz-spin maintaining the right distance between each other? Meso-Scale: How does a dancing couple adapt its movements in response to other dancers, spatial limitation and its targeted dance direction? Macro-Scale: How does the crowd-density change on the dance-floor? Where and when can we detect clusters of dancers? We will present a simulation model that combines all three levels of abstraction and present our attempts of validation. Based on chosen case studies, we will draw conclu- sions on how clusters and collisions on the dance-floor can be avoided. 1.2. Novelty of the research Apart from medically oriented ones, there have been only a few studies on this topic. Most of them only deal with one of the three previously mentioned scales but do not investigate the entire system. Some studies investigate the physics of one single couple dancing (the Micro- Scale) in order to get insights into competitive ballroom dancing [4,5]. TheMeso-Scaleis analysed in [6], where trajectories of competitive ballroom dancers are tracked, yet the interactions between dancing couples are not analysed. In [7] the Macro-Scale was studied but only investigated from a psychological and not a spatial perspective. In summary, literature that considers all three scales and their effects on each other could not be found. 1.3. Modelling approach We decided to tackle the research question using a model, which is best interpreted as an agent-based modelling approach (ABM) with continuous time update. Creating a model combining all three levels of complexity, as explained in the previous section, excluded any other alternatives: MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 171 First of all, the system can hardly be modelled by one (or a few) aggregated ordinary differential equations, as the model needs to depict temporal and spatial changes. The corresponding extension, partial differential equations, are sometimes used for modelling of systems with interacting, moving particles such as traffic flow [8]or flocking [9]. Unfortunately, the macroscopic description form leads to a loss of information on the Micro- and the Meso-scale. This also holds true for all other modelling strategies in which the individuality of the dancers is lost, for example cellular automata. Consequently, a microscopic modelling concept is required. With respect to microscopic modelling techni- ques, classical discrete event simulation is not handy as well because the system involves continuous movement which cannot be neglected. Therefore, the strategy needs to be continuous as well. Continuous microscopic models are used exemplary for simulating particles, so-called molecular dynamics models. These models would dedicate the whole set of necessary (ordinary) differential equations for each particle (see [10]). Other approaches that model the individuals in a similar way are so-called gravity models [11] or the classical car-following models [12]. These approaches have the same foundation, and the behaviour and interaction are based on (differential) equations for each entity. Unfortunately, modelling every dancer or dancing couple solely via a continuous (system of) kinematic differential equations is not enough: the beat of the music poses for a discontinuity at which dancers need to perform a complicated decision-making process. Consequently, we chose to extend the ideas of the mentioned three modelling approaches by complicated time and state events, that involve classical decision-making processes. In [13] a pedestrian dynamics model is presented, that is based on the latter idea. As the author denotes the presented model by the term ‘agent-based’, we decidedtodosoaswell. The generic concept of an agent-based modelling [14] is a suitable framework to describe such highly complex models. It does neither involve any restrictions with respect to time dynamics nor with respect to the decision-making processes involved. Therefore, it is easily capable of combining continuous and discrete processes on a microscopic scale. Moreover, this flexible strategy is already famous for depicting group dynamic effects [15,16] and stands out with a very natural and picturesque description of the modelled system. One could call the model presented in this work a hybrid agent-based model: time-continuously updated with time- and state-events. 2. Model description To increase comprehensibility of the agent-based model (ABM) its description follows the ODD (Overview, Design Concepts, Details) protocol [17,18]. To keep the description of this quite complicated model as short and readable as possible, certain details are found in Appendix A for reproducibility purposes. 2.1. Overview 2.1.1. Purpose The ABM model simulates an arbitrary number of couples dancing in a ballroom to several beats of waltz music. The model allows to perform scenario tests like mixing dancers with different skill levels, comparing different distributions on the dance-floor as 172 M. BICHER ET AL. well as test different shapes of the ballroom. The results of these tests can furthermore be used to draw qualitative and partially quantitative conclusions about the dynamics of a waltz-dancing crowd. 2.1.2. Entities, state variables and scales In this model, each agent represents one dancing couple in reality. As we are interested in spatial interactions it is important to consider the physical properties of the dancing couples. Therefore, the physical appearance of a dancing couple is modelled by a solid cuboid with specific measurements and mass M, representing the two dancers common mass. The agent’s cuboid has two possible rotational axes R and R which represent the A B physical body centres of the two dancers and a distance d in between, as shown in Figure 1. The rotational axes are defined as two-dimensional coordinates relative to a global coordi- nate system. Note, that there is currently no specific purpose for the vertical dimension of the cuboid as all movement takes place in the plane. Yet, this feature allows 3D visualiza- tions and possible model extensions in the future. The movement of agents is, of course, inspired by waltz spins in reality: Agents may lock the position of one of the axes and apply a kinetic energy E to perform a spin around it. Periodically switching the rotation axis, the agent moves forward (see Figure 2). Hence, the coordinates of the rotational axes R ¼ R ðtÞ and R ¼ R ðtÞ change over time. A A B B The agents’ spinning behaviour depends on their individual parameters. These include, most importantly, a maximal kinetic energy E an agent can ‘produce’ for max spinning, as well as a default spinning angle α , which purpose will be explained later. def By default, agents apply their energy to spin in the clockwise direction, yet some of the agents are capable of spinning counterclockwise as well. We will describe this ability by the value of a boolean variable left. The parameters M; E ; α , left and d as well as the max def dimensions of the cuboid stay constant during simulation, while the kinetic energy applied for each spin will vary. The dance-floor represents the environment of the agents. As its appearance also influences the crowd’sbehaviour,wewilltestdifferent shapes of it. It is fully con- tinuous and every position on the dance-floor is reachable. In order to make agents ‘dance’ a music piece, represented as a sequence of beat-times (i.e. time instances of the first beat of every measure), is chosen, defining the length of the simulation run. As mentioned, the physical positions R ðtÞ and R ðtÞ of the agents are updated A B R R A B . . Figure 1. Abstraction of a dancing couple via a cuboid (dashed line) with two rotational axes R and R A B seen from above. The black and white circles should visualize the leading and following dance partner. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 173 R R A B . . Figure 2. An Illustration of the periodically switching rotation axis which makes the dancing couple both spin and move forward. The shape of the arc described by the dancing couple is solely defined by the achieved spinning-angle α. continuously (see Section 2.1.3) which makes the model both time- and spatially continuous. A full collection of model parameters and parameter values is found in Table 1. 2.1.3. Process overview and scheduling All agents in this model are treated simultaneously and update their state/position time- continuously according to classical frictionless equations of motion. Within this process, two types of events interfere with the movement of the agents: Table 1. Full list of all model parameters as well as used values or value-ranges, respectively. A value- range indicates, that this parameter is either varied in a case study or is defined by an agent-collective (see last column). Symbol Parameter Value/value-range Unit N total number of agents 2f1; .. . ; 120g 1 Case Study 1 N total number of agents per collective 2f1; .. . ; 120g 1 Case Study 3 t ; i ¼ 1; .. . ; time instants of the music beats see Section 4 s a outermost super-ellipse, first axle 2½17:7; 25:0 m Case Study 4 b outermost super-ellipse, second axle 2½12:5; 17:7 m Case Study 4 p initial distribution power-parameter 2½0:6; 1:4 1 Case Study 2 E maximum kinetic energy 2½228; 240 J see Table 2 max left counterclockwise ability ftrue; falseg boolean see Table 2 d distance between R and R 0:5m A B ð, ; , ; , Þ cube dimensions ð1:0; 2d; 1:7Þ m 1 2 3 m mass of a couple 150 kg Emax E default kinetic energy J α default angle 70 deg. st 1 measure nd 2 measure 174 M. BICHER ET AL. T T 1 2 a a 1 2 Figure 3. Hybrid automaton representation of the simulation model for the dancing couple, where the two events for the scheduling procedure are illustrated. The time event T represents the beat of the music and the state event T represents collisions between agents. Collisions (state event). Whenever the cuboids of two agents touch, a collision occurs. Inelastic collision handling is applied and both agents are deflected. This process is handled by classic collision equations and the event’s outcome can not be influenced by agents. Music-beat (time event). On the first beat of every three-four- or six-eight-time of the waltz music, the movement of all agents is entirely stopped and a new spin is initiated. Agents decide the amount of energy and the direction used in performing their next spin. During these events, agents are actively influencing the simulation via decision processes. Following the dynamic hybrid system (DHS) approach (see [19]), a graphical diagram of the described process can be depicted. Hereby, each agent, explained above, can be described by an (Alur-Henzinger) hybrid automaton (see Figure 3). Vertex (control mode) A describes the physical dynamics of the spinning movement and uses the state vector x ¼ðR ; R ; φ; ωÞ ; x y where ðR ; R Þ is the position of the current rotation axis of the agent, φ is the angle of the x y movement and ω the angular velocity. The mathematical description of the hybrid automaton requires event formulations, which are given by T ¼ Aa A and T ¼ Aa A: 1 1 2 2 Transition T corresponds to the time event Music-beat and T to the state event Collisions. 1 2 While T can be described as a typical time-event, using the classic DHS formulation, the state event cannot, as agent-agent-interactions are involved. Hereby, different to a classical DHS model, the guard and jump of T need to consider the states of all agents simulta- neously: The guard is represented on the interaction rules of the agents (collision detection) and the jump corresponds to the collision handling of the two involved agents. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 175 All three components, A, T and T are explained in Section 2.3 in more detail. 1 2 Note, that for implementation of such a system, time is usually discretized again. This is done by an ordinary differential equation solver (see Section 3.1). 2.2. Design concept 2.2.1. Basic principles From the modelling point of view, the most interesting feature of waltz dancing is the unusual way dancing couples move forward on the dance-floor: After spinning for one measure of waltz music, the orientation of the two partners is more-less switched and the couple, seen as a whole, has moved into the direction of their first step. After one more measure the orientation of the partners is almost back to the original setting, yet the couple itself has continued moving into the (more-less) same direction, compare Figure 2. Consequently, the movement of the whole couple emerges from the individual spin- ning movement of the partners and it is quite challenging to influence it: By adapting the speed of the rotation, the partners can not influence the velocity of the whole couple’s movement, but their orientation after the spin – they can ‘steer’ this way, but can neither break nor accelerate. Also, changing the spinning direction from the easier clockwise ‘Rechtswalzer’ to the much more sophisticated counterclockwise ‘Linkswalzer’ will not influence it. The actual velocity can only be changed by modifying the length of the individual steps, a parameter which is very difficult to change without getting problems with the synchronicity of both partners, e.g. physiological limits, stepping on partners feet, problems with keeping the angular momentum, etc. This observation founds the main concept of the simulation model: As agents are unable to change the distance of their rotation axes on run-time, they cannot influence how fast they actually proceed on the dance-floor. Yet, modifying their angular momen- tum and the direction of their spin, they can decide (up to some extents) where the whole couple is facing after the next measure. This type of steering is extremely unusual and becomes even more strange considering that aiming to the right or to the left with a clockwise spin makes a huge difference. 2.2.2. Emergence There are two elements in the present model which can be interpreted as sources of so- called emergent behaviour. On the lower level, the aforementioned trajectory of the agents can be considered as emerging from the individual spins of the agent. Even though the mathematical relation between individual angle and orientation of the agent after one spin is quite easy, the achieved individual angle depends on the kinetic energy applied and, in particular, the length of the music beat. For a music piece, especially the latter is usually not constant and the agents need to adjust the kinetic energy to compensate it (see Section 2.2.3). On the higher level, the interaction between the agents via collisions is a second source for emergent behaviour. 2.2.3. Objectives In the present model, agents only follow one objective: to stay on track. The tracks of the agents are defined as a set of steady curves that lead around the centre of the ballroom. In the real system, dancers are trying to follow more or less concentric paths around the 176 M. BICHER ET AL. centre of the ballroom in the counterclockwise direction. Based on observations, we found that neither ellipses nor rectangles sufficiently describe the shapes of these paths. The prior leaves the corners of the rectangular ballroom completely unused while the latter unrealistically over-exploits them. Hence, we decided to model the tracks as a family of fourth-order curves, so-called super-ellipses 4 4 x y þ ¼ c; c 2½0; 1: (1) 4 4 a b Hereby 2a and 2b correspond to the length and width of the ballroom. Moreover, any point ðx ; y Þ2 R on the dance-floor corresponds to one unique p p curve of this family. Hence, whenever an agent makes a decision about its spinning behaviour, it determines the super-ellipse that contains its current position, and aims to remain on this curve after the next spin. For detailed information on how this is done, the reader is referred to Section 2.3.4. Finally, the modelled agents do not aim to avoid collisions with other agents, though we would find such behaviour with dancers in reality. We found it feasible to neglect this feature because (1) modelling would be too complex, and (2) it can be substituted by a correct interpretation of the results: A collision of two agents in the model can be interpreted in the real system by either (a) at least one of the two couples is unable to finish its spin, (b) at least one of the two couples is performing some collision-avoidance move, or finally (c) an actual collision of two couples. 2.2.4. Sensing As mentioned in Section 2.1.3 an agent can only actively influence its motion at the music- beat events. Hereby the agent senses its velocity and its current position in the ballroom as well as thegraph of thecorrespondingsuper-ellipse(seeabove). The agent uses this information to calculate a target angle which it aims to spin in the course of the music- beat. Yet, in the course of this decision process, the agent is not aware of any other agents around. In case of a collision event, the agent is aware of the other agent and its motion, but cannot use this information for any future decisions. 2.2.5. Interaction In this model, the only form of interaction between agents is the collisions mentioned above (for more details on collisions, see 2.3.3). 2.2.6. Learning Due to occasional tempo changes in classical music, it is difficult for an agent to stay on track. In particular, the mapping α : energy applied ¼ E 7! αðE Þ¼ angle spinned kin kin depends on many additional parameters, namely (1) length of the current measure, (2) mass of the couple, (3) whether the spin is clockwise or counterclockwise, (4) distance between the rotation-axis and (5) whether a collision occurred or not. Although the physical relations between (1), (2), (3) and (4) are well known, (1) and (5) are difficult to predict. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 177 Table 2. The different skill level characterized through maximal kinetic energy and ability to switch spinning direction. The rightmost column contains the estimated maximal spinning angle per second. Collective name E (maximum energy) left (counterclockwise spin) est. angle range per second max perfect agents E ¼ 240 J true \ ½188; 188 max good agents E ¼ 240 J false \½0; 188 max slow agents E ¼ 228 J false \ ½0; 183 max So we decided to model this function as a black box combined with a learning process: At the beginning of the simulation, agents are unaware of the parameters of the black-box and perform their initial spin with a default energy E . Based on the outcome of each passed spin, they learn how much energy they have to apply to achieve a certain target angle. This process is described in 2.3.4 in more detail. 2.2.7. Stochasticity The model is a perfect example of an initial-value-stochastic microscopic model (com- pare with the classification proposed in [20]). While the initial positions and orientations of the agents are stochastic, the model dynamics, i.e. the decision rules of the agents as well as the kinematic movement, are deterministic. Hence, the same initial configuration of agents lead to the same simulation outcome. 2.2.8. Collectives In order to test the heterogeneity of the crowd more easily, we decided to summarize agents with equivalent traits to collectives. Hereby, we specifically refer to two para- meters, that define the agent’s ‘dancing skills’: The maximum kinetic energy E max applicable for spinning and the agent’s ability to perform counter-clockwise spins left (see Table 2). Note, that the estimated angle-range in this table is given for spins of precisely 1 s. Hence, a slow agent might get troubles performing 180 degree spins, if the measure of the waltz music is shorter than 1 s. To avoid misunderstanding, we want to clarify that these groupings are designed for analysis purposes only and do not have any impact on the model logic. 2.2.9. Observation First, the time-dependent positions R ðtÞ and R ðtÞ are recorded at each beat-event. A B After the simulation, this information is investigated for distribution of dancers with different skill levels. Second, the total number of agent-contacts, an agent has within the course of one simulation run, is tracked. We analyse the latter with respect to average collisions per agent and per skill-level. 2.3. Details 2.3.1. Initialization We initialize the model with a total number of N ; i 2f1; 2; 3g agents of each of the three collectives, explained in Section 2.2.8 and their skill parameters set as specified. We define N ¼ N þ N þ N as the total number of agents in the model. Furthermore, initial 1 2 3 positions for the agents are sampled. Given a uniformly drawn orientation in ½0; 2π of 178 M. BICHER ET AL. the agent, the centre point of the agent, i.e. ðR ð0Þþ R ð0ÞÞ, is set according to A B a specific distribution: Let X ¼ðX ; X Þ be a uniformly distributed point in 1 2 4 4 x x 1 2 ðx ; x Þ : þ  1 ; (2) 1 2 4 4 a b then the couple’s initial position is given by TðXÞ with p1 4 4 4 * * X X 1 1 TðXÞ¼ X  þ ; (3) 4 4 a a with a power-parameter p 2 R . Consequently, agents are distributed fully uniformly on the super-elliptical region for p ¼ 1. For p > 1 agents are more likely to be close to the centre of the room and for p< 1 vice versa. We will investigate the impact of this parameter in one of the case studies. 2.3.2. Submodel a – physical description of agent-movement The solely physical process is modelled frictionless and with full respect to energy con- servation. By applying angular velocities to the agents’ collision bodies, all angular and translational accelerations are and remain equivalent to zero at all times. In between any two events, the angular and translational velocity is constant and all positions change uniformly. Directly after the execution of a music-beat event (see 2.3.3) only the angular velocity around the rotation axis is non-zero, leading to a rotation with constant speed. In contrast, after the execution of a collision event also translational velocities can be non- zero (see 2.3.4) and the movement becomes a mixture of non-accelerated translation and rotation around a non-centred axis. 2.3.3. Submodel T – music beats Whenever a first beat of the three-four or six-eight time of the music occurs, i.e. at every time instant t , every agent starts a new spin around the switched spinning axis. Hereby, it decides on the energy it applies for the next spin and, in case the agent has an advanced skill level, whether it performs the spin clockwise or counterclockwise. This process can be split into the following steps: (1) The agent stops moving and all velocities are set to zero. (2) The rotation axis is switched. (3) Agents perceive their current position in the environment and find their target path. (4) Agents decide on how much energy they want/can apply for their next spin to reach the target and if they want/need to change spinning direction to succeed. (5) The kinetic energy is transformed to a new angular velocity and applied to the cuboid. A detailed description of steps 3–5 can be found in Appendix A. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 179 2.3.4. Submodel T – collisions In case the cuboids of two agents get in contact with each other, i.e. the surfaces of any two cuboids of agents intersect, inelastic collision-equations, i.e. collisions with restitu- tion coefficient equal to zero, are evaluated to solve it. This is feasible as collisions between humans usually lead to bending of body parts that absorb almost all the energy of the push. Equations for inelastic collisions in 2d and 3d are well known (e.g. see [21]) therefore they will not be described any further. 3. Simulation model analysis 3.1. Implementation In order to simulate the time-continuous elements of the model, the open-source 3D multi-body physics engine Bullet [22] was used. This powerful engine makes it possible to evaluate the collisions of hundreds of agents in reasonable time. As visual output is very important for face-validation purposes, we used Bullet’s interface provided by the graphics engine Panda3D [23]. Panda3D is primarily intended for the development of PC games (e.g. see Pirates of the Caribbean Online [24]) therefore the output of the simulation is nicely visualized in 3D and can be observed on run-time. In order to get a more realistic impression, the agents’ cube-shapes are replaced by more realistic 3D models of dancing couples, as seen in Figure 5. Both Panda3D and Bullet are freely available and are used as importable Python3 modules. A code for running the base model is available on GitHub https://github.com/MBicher/waltz_simulation. 3.2. Algorithms, complexity and numerics As the solution of the frictionless equations of motions does not involve any problems, Bullet’s internal semi–implicit Euler method is sufficient for integrating the occurring equations with a fixed time-step dt ¼ 0:025 ½s. Moreover, Bullet stands out with its fantastic features of collision detection and handling. Hereby it uses a Sweep and Prune algorithm to detect, and a variant of the projected Gauss–Seidel method to sequentially handle all occurring collisions [25]. Alike most simulations that involve the interaction of individuals, the computational and temporal complexity of the algorithm rises non-linearly when the number of agents is increased. Considering N as the total number of agents in the simulation, this procedure would lead to the computational effort of OðÞ N . Luckily, Bullet combines a Sweep and Prune algorithm and axis-aligned bounded boxes (AABB) to pre-process collisions and therefore reduces the computational effort to OðÞ N [26]. These observations do not take into account what happens in the course of the simulation on the dance-floor. The number of collisions that actually take place increases quadratic as long as the space of the dance-floor is kept constant (see next section). Consequently, the computational effort will increase by square as well. Apart from the long-term behaviour, our experiments revealed that up to an agent- count of about 50  70 the simulation time almost remains constant. Up to this number, the standard processes of Panda3D and Bullet, which occur during time-evolution, 180 M. BICHER ET AL. dominate the efforts for collision detection, evaluation and other iterative procedures. Thereafter, the computational effort increases as expected. 3.3. Value calibration Besides parameters actually measured from observing the real system (mass of a couple, size of the ballroom, . . .), Table 1 also shows one particular parameter, calibrated from simulation runs, namely the skill-dependent maximum kinetic energy E . This was max done by a simple parameter sweep comparing the applied kinetic energy with the maximum achievable angle. In other words, the second column of Table 2 was calibrated using the values in the fifth. 3.4. Verification To ensure that all implemented procedures work as modelled, a couple of verification scenarios has been set up and tested. These include: ● scenarios with only one dancing couple following a given trajectory (circle, ellipse, super-ellipse, sine-wave . . .), to make sure that the goal-seeking behaviour works, ● scenarios with two colliding dancing couples in different settings, to make sure that collisions are treated correctly, and ● scenarios with several dancing couples perfectly aligned in a circle, to investigate the propagation of numerical errors. 3.5. Validation Like most ABMs that include complex human behaviour, our dancing model is very difficult to validate. So far, we restricted our attempts to qualitative face validation together with dancing experts comparing open-source videos and pictures of famous balls in Vienna: Opernball, Zuckerbäckerball, Kaffeesiederball and, of course, the TU Ball, hosted annually by our university. Feedback of the face-validation was fully positive, as many group-dynamic effects, observed in reality, could also be reproduced in the virtual ballroom (see Section 5). In the future, we aim to perform additional validation attempts via trajectory-tracking devices. We plan to use motion capturing software to track a few number of dancers in a laboratory setting filmed from above, as well as GPS-based devices to track single trajectories of volunteers in a real ballroom situation. 4. Simulation scenarios and results In order to draw conclusions for the real system, we specified four case studies in which we varied specific model parameters and track the defined observables: the average number of collisions per agent and their time-dependent positions. Although the model only contains randomness within the initial conditions, the model is not only a stochastic one, but we expect the model dynamics to depend heavily MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 181 on the initial setting as well. Consequently, we need to perform a Monte Carlo simulation to make valid statements about the outcome: Each defined scenario is simulated 80 times with different random-number generation seeds. Around the gathered sample mean of the results of the simulation runs we will always state estimates for a confidence interval that contains the unknown real mean of the stochastic simulation. This estimate is based on the number of samples, the conver- gence speed of the sample mean, the sample variance of the results, and a Gaussian distribution assumption: With Φ being the cumulative probability function of the standard normal distribution and X and s being the sample mean and variance of N independent Monte Carlo samples, it is shown [27] that ! rffiffiffiffiffi 1  p s 1 N P X  μ Φ  p: (4) 2 N By adding these confidence interval estimates, we want to convince the reader that observed trends seen in the scenarios are not coincidental, exemplary due to fluctuations. We perform all scenarios for manually extracted measures of Waltz No. 2 by Dmitri Shostakovitch, which is 3:55 minutes long and has a total of 222 measures. The time- spans between consecutive measures vary between 0:83 s and 1:25 s and average in 0:97 s. Consequently, slow agents are not able to dance 180 degree spins all the time. 4.1. Case study 1: varying the number of dancing couples Our first subject of interest is the agent-density. Though it is clear that increasing the total number of dancers has a negative influence on the overall performance of the waltz, the rate of influence can only be quantified by simulation. Our simulation model reveals that the total number of agents N on the (same) dance- floor influences the average collision count per agent with approximately quadratic order. This can be justified analytically by thinking about the average number of collisions as the product of number of agents and the fraction of occupied space on the dance-floor. As the latter depends linearly on the number of agents N, the product grows approximately with N . Figure 4 shows a graph of the collision count for N from 30 to 130 with 80 simulation runs each. Hereby, only perfect agents have been considered (see Section 2.2.8). That means that all modelled agents are sufficiently skilled to perform spins with at least 180 degrees and are able to spin counterclockwise as well. The initial distribution is uni- form (p ¼ 1). The shaded area around the sample mean, calculated from 80 Monte Carlo runs, shows the mentioned estimated confidence interval that contains the real, unknown mean of the actual distribution of the simulation observable with 95% probability. Moreover, Figure 5 gives a snapshot of the simulation for N ¼ 100. The upper part of the figure shows the picturesque 3D-output of the simulation provided by Panda3D. In the lower part of the figure, the positions of the rotation axis R ðt Þ of all agents are B end depicted as red dots, and the corresponding trajectories R ðt Þ; i 2f0; .. . ; t g are given B i end as lines fading from black (very recent) to white (long time ago). For the regarded simulation setting with N ¼ 100, it shows quite smooth super-elliptic lines which 182 M. BICHER ET AL. empiric mean value of 80 simulation runs 95% confidence interval for the actual mean 30 40 50 60 70 80 90 100 110 120 130 140 150 number of agents N Figure 4. Average number of collisions per agent within 222 measures of waltz-music under variation of the total numbers of agents. The shaded area shows a confidence interval around the sample mean of the Monte Carlo analysis that contains the real, unknown mean with 95% probability (see (4)). Figure 5. Snapshot of the simulation with 100 perfect agents after spinning 222 measures. The upper image shows a snapshot of the visualization from Panda3D for case study 1, the bottom image shows the corresponding traces of the agents on the dance-floor, to be precise, the positions of R ðt Þ; i 2f1; .. . ; 222g. B i number of collisions per agent MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 183 indicate that agents hardly collided. It seems, the crowd circulated quite fluidly around the centre in this simulation run. 4.2. Case study 2: varying the initial distribution In the second case study, we investigate the impact of agents’ initial positions. Hereby the power-parameter p is varied (see Equation (3)) to find out, which initial distribution results in the minimal number of collisions. Parameter p is varied in steps of 0:1 between 0:3 and 1:3 for N ¼ 100 perfect agents. The mean value of collisions is drawn out of a Monte Carlo analysis with 80 simulation runs each. Results are shown in Figure 6. Although we would assume that a fully uniform distribution (p ¼ 1) of dancers would result in the most homogeneous crowd, the least number of collisions occurs with a parameter around p ¼ 0:6. Consequently, being placed slightly more densely to the border is advantageous with respect to avoiding conflicts. Figure 7 shows the difference between the optimal distribution suggested by the simulation compared with a fully uniform one. 4.3. Case study 3: varying the skill level Before performing tests with the heterogeneity among the crowd, the impact of reducing the skill level of all dancers, is investigated. It seems feasible that the skill to switch spin direction and to achieve larger maximum angles, should increase the couple’s ability to react to disturbances. Hence, we would expect the total number of collisions to increase when these skills are reduced. empiric mean value of 80 simulation runs 95% confidence interval for the actual mean 0.30.40.50.60.70.80.91.01.11.21.3 initial distribution coefficient c Figure 6. Impact of the initial distribution on the number of collisions. The average number of collisions per agent (gathered from 80 simulation runs each) is plotted against the power-parameter p. The shaded area shows a confidence interval around the sample mean of the Monte Carlo analysis that contains the real, unknown mean with 95% probability (see (4)). number of collisions per agent 184 M. BICHER ET AL. 10 10 5 5 0 0 −5 −5 −10 −10 −20 −10 010 20 −20 −10 010 20 ball-room length [m] ball-room length [m] Figure 7. Comparison of 3000 points drawn according to a fully uniform distribution (left, p ¼ 1) and the experimentally determined ‘optimal’ distribution for dancing (right, p ¼ 0:6). Table 3. Average number of collisions per agent for different skill levels and total agent counts. The expressions in the square brackets describe an estimated 95% confidence interval around the depicted sample mean of 80 simulation runs (see Equation (4)). Skill-level N ¼ 60 N ¼ 80 N ¼ 100 N ¼ 120 N perfect agents 48.96 [±3.40] 79.49 [±4.02] 121.42 [±6.48] 170.52 [±7.24] N good agents 58.20 [±4.76] 101.37 [±7.14] 178.31 [±9.88] 292.67 [±14.47] N slow agents 152.89 [±11.47] 237.50 [±11.74] 384.65 [±12.13] 562.70 [±16.59] mix of all three 119.81 [±6.23] 184.10 [±7.46] 304.68 [±12.18] 429.70 [±14.89] The simulation setting of Section 4.1 is slightly changed, (1) removing the agents’ ability to switch the spin direction (good agents) and moreover (2) reducing their maximum achievable kinetic energy (slow agents). Afterwards, heterogeneity is tested as we perform simulations with a homogeneous mix of agents with each of the three skill levels: A total of bN=3c perfect, bN=3c good, and bN=3c slow agents was homogeneously distributed (p ¼ 1) on the dance-floor for different total numbers N. The most important model outcomes are summarized in Table 3. The results clearly show that worsening the skill level leads to vastly more collisions occurring in the model. This effect is particularly significant as soon as the total number of agents on the dance-floor increases. Interestingly, the collision counts for the homo- geneously mixed scenarios are much more similar to the slow agent scenarios than the perfect or good agent ones. Apparently, this is a situation in which, proverbially, good apples, put into a bad situation, partially become bad apples. Anyway, in the course of this study we recognized that good agents and in particular slow agents tend to cluster in the centre of the ballroom due to their missing ability to steer to the right, i.e. towards the border of the ballroom. This situation is depicted via a Gaussian kernel-density estimation (KDE)of8000samples (80simula- tion runs with 100 agents each) in Figure 8. The upper heat-maps show the KDE for the agents’ positions at the beginning of each simulation, the lower show the KDE at the end. As a direct application of the previous finding, the influence of the initial positions of 50 slow agents together with 50 perfect agents is analysed. We assumed that, while on their way to the middle, slow agents would drag perfect agents with them leading to unneces- sarily many collisions. ball-room width [m] MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 185 good agents slow agents t =00 : 00 t =03 : 55 Figure 8. Kernel density estimation (KDE) of the positions of 100 agents in 80 simulation runs. The upper plots show the KDE at the start of the simulation while the lower plots depict the situation at the end of the music piece. The two left plots show the KDE for good agents, the two right plots for slow agents. Both collectives tend to cluster in the middle of the ballroom. Table 4. Average number of collisions compared for two different initial settings for a heterogeneous crowd. The expressions in the square brackets describe an estimated 95% confidence interval around the depicted sample mean of 80 simulation runs (see (4)). Scenario Avg. collisions per perfect agent Avg. collisions per slow agent Avg. collisions per agent ‘reference’ 325.05 [±13.68] 333.18 [±11.52] 329.12 [±10.55] ‘slow start at centre’ 74.26 [±4.51] 543.83 [±25.18] 309.04 [±12.82] In our ‘reference’ scenario the perfect and slow agents are mixed homogeneously and their position is initialized uniformly (p ¼ 1) on the dance-floor. In the alternative scenario (‘slow start at centre’), the initial positions are determined the same way, yet, afterwards, agents are resorted so that the initial positions of all slow agents are closer to the centre than those of all perfect agents. The results are depicted in Table 4 and partially confirm the assumptions: Leaving the outside tracks to perfect agents significantly reduces the number of collisions for them, but also vastly increases them for the slow agents. Interestingly, though, with respect to the overall average number of collisions per agent, both scenarios are almost equivalent. As the latter was initially defined as a ‘measure for fluidity’ in the ballroom, the alternative setting is not convincingly better than the reference. 4.4. Case study 4: varying the ballroom Run-time observations of simulation indicate that agents show immense problems staying on their track while spinning along the longer axis of the super-ellipses, as they need to achieve almost 180 degree per spin. Therefore, in this final case study, we analyse the impact of the dimensions of the ballroom on the collision count. 186 M. BICHER ET AL. Table 5. Average number of collisions compared for different ratios of the ballroom. The expressions in the square brackets describe an estimated 95% confidence interval around the depicted sample mean of 80 simulation runs (see (4)). Agents a : b ¼ 2 : 1 a : b ¼ 3 : 2 a : b ¼ 1 : 1 N ¼ 100 perfect agents 121.42 [±6.48] 104.29[±5.52] 101.75 [±5.84] N ¼ 100 slow agents 384.65 [±12.13] 258.70[±10.96] 202.63 [±8.14] a : b =1 :1 a : b =2 :1 Figure 9. Comparison of a quadratic and a rectangular ballroom. The colours of the displayed super- ellipses show how much angle an agent needs to spin per measure in order to stay on track. Keeping the total area of the ballroom constant, the ratio a : b of the two axis is varied and the total number of collisions is observed. Results are displayed in Table 5. The results show that the quadratic (a : b ¼ 1 : 1) ballroom helps especially slow agents as the regions in which they are forced to spin angles close to 180 degree are much shorter, as visualized in Figure 9. 5. Discussion Though especially quantitative validation is difficult, analysis of the model and the performed case studies (henceforth, short CS) leave enough scope for interesting con- clusions for the real system. Some of them confirm already known properties about the real system, while some others are, up to some extent, surprising. 5.1. Learn to spin counterclockwise In CS3 the feature of dynamically switching spinning direction, from clockwise to counterclockwise and vice versa, is an essential advantage compared to other dancers. It significantly decreases the total number of collisions due to the agents’ additional skill of ‘steering’ to the right. In reality, counterclockwise spinning is more difficult as it involves advanced foot- work. Being able to dynamically switch between clockwise and counterclockwise is even more complicated as dance-partners need to communicate properly within fractions of MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 187 a second. Yet, well-trained dancers are able to perform these kind of movements and are sometimes seen to use this skill to avoid conflicts with other dancers, exemplary by making a swift counterclockwise spin when the path on the left is blocked. Anyway, our model would suggest, that it is definitely worth training. 5.2. Spin fast and reduce the distance CS3 also showed that it is beneficial for an agent if it is capable to spin large angles, which is achieved by creating a high angular velocity v at the start of the spin. The model is ang based on standard laws of kinetics and kinematics. Therefore, three factors influence the maximum angular velocity an agent is able to apply for spinning: mass (moment of inertia), kinetic energy and distance between the two rotation axes R and R . Therefore, A B pffiffiffi reducing the maximum kinetic energy by a factor of x, increasing the mass by a factor x, pffiffiffi or increasing the distance between axis by a factor of x yields the same kinetic effects in the model. This fact is easily translated to the real system: if dance partners double the space between themselves, they need to develop four times more kinetic energy – i.e. apply more force, be quicker, make more accurate steps, etc. – to spin with the desired angle. Thus, in order to dance in an energy efficient way (which can be important to stand through a long glittering ball night) a couple needs to shift the two rotational axes as closely together as possible. 5.3. Do not distribute uniformly CS2 showed that a uniform initial distribution on the dance-floor is not the most beneficial with respect to collision avoidance, but the one depicted on the right side of Figure 7 is. Clearly, in reality, this distribution is not easy to achieve. In case all dancers would collaborate, they could follow the strategy: distribute evenly and then make one step towards the border of the room. 5.4. Stay away from the centre In most cases, the crowd is huge and does not collaborate. Hereby, CS3 gives some answers to where dancers should position themselves in order to prepare for a peaceful waltz. Experienced dancers should always try to find a spot closer to the borders, while novice dancers have two options: In case they want to be good team-players, they should place themselves closer to the centre in order to give experience dancers more freedom. By doing so, they pay the price of more ‘traffic’ around themselves. More selfish novices would pick a spot closer to the border, risking experienced dancers complaining about them. 5.5. The more quadratic, the better The model also provides insights for the organizers of balls. The results of CA4 indicated that organizers of balls should, in case of direct conflict of two equivalently sized 188 M. BICHER ET AL. rectangular rooms, prefer the more quadratic dance-floor. Also, splitting one large and long rectangular room into two more quadratic ones might be helpful. 5.6. Avoid overcrowded ballrooms Last but not least, it is generally a good idea to avoid overcrowded ballrooms. Although there is no real need for a simulation model to give this advise, our agent-based approach helps to understand the dynamics of overcrowding a little bit better. As seen in Figure 6 the collision count increases quadratically with the dancer-density, an effect which becomes even more prominent when not only perfectly trained, but also less skilled dancers are involved. 6. Conclusion In this work, we presented a highly complex individual-based model for ballroom dancing which stands out from other ABMs by several interesting features. Firstly, due to the interplay of music cadence and physics, an innovative hybrid update strategy was chosen. The dynamic update of a dancing couple agent combines aspects of time-/state- events at which typical decision-making processes are executed, and the classic differ- ential equations of motion which allows the agent to move. Secondly, planning the couple’s future track on the dance-floor requires solving a difficult unique control problem: The path of the whole couple needs to be steered by adjusting the angular velocity. The presented simulation model is the perfect environment for analysis of waltz dancing in a virtual laboratory and a lot more case studies with potentially interesting outcomes could be made. With respect to limitations, the model cannot and should not be used to observe the dancing behaviour of a single dancing couple without the context of the dancing crowd. For this, the individual behaviour of the dancers (i.e. steps sequence, body movement, . . .) is modelled in a too simplified way and would require adifferent modelling approach. Ripped costumes and even injuries are possible outcomes of colliding waltz couples and in most cases, the less skilled couple is blamed. Even though the results of the case studies show that, though the skill level of the dancers matters, this parameter is on the macro- scopic scope less influential than the density of dancers on the parquet. Hence, in order to avoid collisions with other dancers, the easiest solution is to wait until after midnight when the dance-floor is empty enough to dance peacefully – no matter the skill set. Acknowledgments We would like to give our special thanks to Prof. Dr Felix Breitenecker for many decades of not only academic, but also cultural education. His annual invitations to the traditional TU Ball inspired us to implement this model. The authors acknowledge TU Wien University Library for financial support through its Open Access Funding Programme. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 189 Disclosure statement No potential conflict of interest was reported by the authors. References [1] S.J. Wolf, Beweis daß das Walzen eine Hauptquelle der Schwäche des Körpers und des Geistes unserer Generation sey, Hendel, Johann Christian, Halle, Germany, 1799. [2] S. Grobler, Thelifeand deathofthe pianowaltz, Ph.D. diss., University of Cincinnati, [3] S. Wolfram, Complex systems theory, Emerging syntheses in science: Proceedings of the founding workshops of the Santa Fe Institute, 1988, p. 182. [4] T. Shioya, Analysis of swing movement in ballroom dancing, Multidiscip. Digital Publ. Inst . Proc. 2 (2018), pp. 224. [5] K. Hachimura, K. Takashina, and M. Yoshimura, Analysis and evaluation of dancing movement based on LMA, in Robot and Human Interactive Communication, 2005. ROMAN 2005. IEEE International Workshop on, IEEE, Nashville, Tennessee, 2005, pp. 294–299. [6] P. Zaletel, G. Vuckovic, N. James, A. Rebula, and M. Zagorc, A time-motion analysis of ballroom dancers using an automatic tracking system, Kinesiologia Slovenica 16 (2010), pp. 46. [7] D. Gaertner, K. Clark, and M. Sergot, Ballroom etiquette: A case study for norm- governed multi-agent systems, in Coordination, Organizations, Institutions, and Norms in Agent Systems II, Springer, Hakodate, Japan, 2006, pp. 212–226. [8] T. Musha and H. Higuchi, Traffic current fluctuation and the burgers equation, Jpn. J. Appl. Phys. 17 (1978), pp. 811. doi:10.1143/JJAP.17.811. [9] F. Cucker and S. Smale, Emergent behavior in flocks, IEEE Trans. Automat. Contr. 52 (2007), pp. 852–862. doi:10.1109/TAC.2007.895842. [10] P. Deuflhard, J. Hermans, B. Leimkuhler, A.E. Mark, S. Reich, and R.D. Skeel, Computational Molecular Dynamics: Chellenges, Methods, Ideas, Springer, Berlin Heidelberg, Germany, 1999. [11] A. Sen and T.E. Smith, Gravity Models of Spatial Interaction Behavior, Springer, Heidelberg, Germany, 1995. [12] L.A. Pipes, An operational analysis of traffic dynamics, J. Appl. Phys. 24 (1953), pp. 274–281. doi:10.1063/1.1721265. [13] T. Schelhorn, D. O’Sullivan, M. Haklay, and M. Thurstain-Goodwin, STREETS: An agent-based pedestrian model, in CUPUM’99 Computer in urban planning and urban management on the edge of the millennium: Proceedings of the 6th International Conference, 01, Franco Angeli, Venice, Italy, 1999, pp. 1–12. [14] C.M. Macal and M.J. North, Tutorial on Agent-Based Modeling and Simulation Part 2: How to Model with Agents, in Proceedings of the 2006 Winter Simulation Conference, IEEE, Monterey, California, 2006, pp. 73–83. Available at http://www.informs-sim.org/wsc06pa pers/008.pdf. [15] E. Bonabeau, From Classical Models of Morphogenesis to Agent-Based Models of Pattern Formation, Artif. Life. 3 (1997), pp. 191–211. doi:10.1162/artl.1997.3.3.191 [16] R.M. Axelrod, The Evolution of Cooperation, Basic Books, New York, 1984. [17] V. Grimm, U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T. Grand, S.K. Heinz, G. Huse, A. Huth, J.U. Jepsen, C. Jrgensen, W.M. Mooij, B. Mller, G. Peer, C. Piou, S.F. Railsback, A.M. Robbins, M.M. Robbins, E. Rossmanith, N. Rger, E. Strand, S. Souissi, R.A. Stillman, R. Vab, U. Visser, and D.L. DeAngelis, A standard protocol for describing individual-based and agent-based models, Ecol. Modell. 198 (2006), pp. 115–126. Available at https://linkinghub.elsevier.com/retrieve/pii/S0304380006002043. 190 M. BICHER ET AL. [18] V. Grimm, U. Berger, D.L. DeAngelis, J.G. Polhill, J. Giske, and S.F. Railsback, The ODD protocol: A review and first update, Ecol. Modell. 221 (2010), pp. 2760–2768. Available at https://linkinghub.elsevier.com/retrieve/pii/S030438001000414X. [19] T.A. Henzinger, The Theory of Hybrid Automata,in Proceedings - Symposium on Logic in Computer Science, IEEE,Piscataway, NJ, United States, 1996, pp. 278–292. DOI:10.1017/ s0963180100007039 [20] M. Bicher, Classification of Microscopic Models with respect to Aggregated System Behaviour, Diss., TU Wien, Vienna, Austria, 2017. [21] D.T. Greenwood, Principles of Dynamics, Prentice-Hall, Englewood Cliffs, NJ, 1988. [22] E. Coumans, Bullet physics engine, Open Source Software, Available at http://bulletphysics. org1 (2010), p. 3. [23] M. Goslin and M.R. Mine, The panda3d graphics engine, Computer 37 (2004), pp. 112–114. doi:10.1109/MC.2004.180. [24] C. Lang, Panda3D 1.7 Game Developer’s Cookbook, Packt Publishing Ltd, Lincoln Road Olton Birmingham, UK, 2011. [25] E. Coumans, Bullet 2.80 physics sdk manual, Available at http://www.cs.uu.nl/docs/vakken/ mgp/assignment/Bullet (2012) [26] D.J. Tracy, S.R. Buss, and B.M. Woods, Efficient large-scale sweep and prune methods with AABB insertion and removal, in 2009 IEEE Virtual Reality Conference, IEEE, Lafayette, Louisiana, 2009, pp. 191–198. [27] S. Ross, Simulation, Knovel Library, Academic Press, Amsterdam, 2012. Appendix A. Detailed Description of Submodel T A.1. Finding the new target (3). As mentioned, every agent first needs to perceive its position to find out which of the super-ellipses (1) currently describes its track. We consider the i-th beat of the music piece and the new rotation axis of the agent is R with R ðt Þ¼ðxðt Þ; yðt ÞÞ (anything stated henceforth extends analogously A A i i i 4 4 for R ). Moreover, the agent’s current track calculates to ðxðt Þ=aÞ þðyðt Þ=bÞ ¼ c . B i i i Based on these observations, the agent solves the root-finding problem 4 4 ðxðt Þþ d cosðαÞÞ ðyðt Þþ d sinðαÞÞ i i c ¼ þ 4 4 a b with respect to the unknown α. It results in the two intersection points of the super-ellipse and a circle with radius d centred in R ðt Þ. Hereby, it gives the agent two possible target points to A i aim for to stay on its current track. As it is convention to move counterclockwise on Viennese ballrooms, the one with the larger polar angle, with respect to the centre of the ballroom, is the correct solution. This solution, henceforth denoted as R ðt Þ, forms the agent’s new target B iþ1 point. In order to calculate how much energy needs to be applied, the agent calculates the new target angle α :¼ \ R ðt ÞR ðt ÞR ðt Þ (measured clockwise). p B i A i B iþ1 A.2. Calculate the new energy and direction (4). As mentioned in 2.2.6, agents memorize their previously performed spins to learn from them. With E ; j 2f1; ... ; ig denoting the previously used kinetic energies for spinning, α ; j 2 j j f1; .. . ; ig denoting the successfully danced angles, and t ; j 2f0; .. . ; ig standing for the lengths of all previously danced measures, the agent forecasts the necessary energy E to achieve angle α p p with the next spin. This is done assuming that the square root of the kinetic energy, the resulting angle and the length of the measure are in a linear relationship: MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 191 Figure A1.. Diagram for deciding the amount of energy E and direction d applied for the i þ 1-st iþ1 iþ1 spin starting at time t . We assume that this spin has rotation axis R with current position R ðt Þ. i A A i A spin around R is calculated analogously. () pffiffiffiffi E ðt  t Þ j j j1 c :¼ Median : j 2f1; .. . ; ig; α  0 : (A1) Using the median eliminates outliers resulting from spins interrupted by collisions. Moreover, we use the median of the last four measured lengths dt :¼ Median ðt  t Þ; j 2fi  4; .. . ; ig (A2) m j j1 to get a feasible guess for the length of the next beat. The values c; dt and α are then used to m p calculate the target energy E using E :¼ c : (A3) dt Furthermore, the target energy E is compared to the maximum energy E the agent is able to p max use for spinning. If E  E , the target energy is applied as E to the new rotational axis in the p max iþ1 form of torque. Otherwise, a couple of escape strategies are attempted which are summarized in a decision-tree structure in Figure A1. 192 M. BICHER ET AL. In case every strategy fails, the agent aims for a default angle α of 70 degree. This angle was def chosen based on observations and experience to avoid a spinning deadlock. Variations of this parameter in the final simulations showed that its precise value has no influence on the simulation outcomes. Finally, in case i ¼ 1 when no previous learning data is available, the agent simply applies the default angle E into the default spinning direction (clockwise). def A.3. Energy to angular velocity (5). Finally, the applied kinetic energy is transformed into angular velocity via sffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffi jv j 2E 1 4E 1 tang kin kin v ¼ ¼ ¼ : ang d m d m d spin Hereby, we assume that the mass of the spinning object m is equivalent to the mass of one of the spin dance partners and therefore approximately half of the mass of the couple. The sign of the expression is determined by the direction of the spin. http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Mathematical and Computer Modelling of Dynamical Systems Taylor & Francis

Modelling a Viennese ballroom: agent-based simulation to investigate complex behaviour

Download PDF
24 pages

Loading next page...
 
/lp/taylor-francis/modelling-a-viennese-ballroom-agent-based-simulation-to-investigate-OahHEyV1mA

References

References for this paper are not available at this time. We will be adding them shortly, thank you for your patience.

Publisher
Taylor & Francis
Copyright
© 2020 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group.
ISSN
1744-5051
eISSN
1387-3954
DOI
10.1080/13873954.2020.1727930
Publisher site
See Article on Publisher Site

Abstract

MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 2020, VOL. 26, NO. 2, 169–192 https://doi.org/10.1080/13873954.2020.1727930 ARTICLE Modelling a Viennese ballroom: agent-based simulation to investigate complex behaviour a,b c c M. Bicher , S. Winkler and A. Körner a b Institute of Information Systems Engineering, TU Wien, Vienna, Austria; dwh Simulation Services, dwh GmbH, Vienna, Austria; Institute of Analysis and Scientific Computing, TU Wien, Vienna, Austria ABSTRACT ARTICLE HISTORY Received 4 March 2019 Dancing Viennese Waltz in one of the great historic ballrooms is an Accepted 6 February 2020 important and indispensable part of Austrian culture. This dance, while being tradition, is quite difficult to perform, especially if the KEYWORDS dance-floor is crowded. There, it is additionally challenging to avoid Agent-based simulation; collisions with other dancers, as they pace through the ballroom at ballroom dancing; complex a high velocity. Dependent on the dancer’s skill level, spinning speed system; Viennese Waltz can be adjusted to succeed. This paper presents an agent-based waltz model which makes it possible to investigate the influence of heterogeneously skilled dancers on the movement smoothness of the dancing crowd. Herein, each agent represents one dancing cou- ple in reality and it moves on the dance-floor by a rotatory motion with periodically switching rotation axes. Interaction between agents occurs via inelastic collisions. By performing a couple of case studies, we analyse and quantify the widespread rumour that the presence of only a few unskilled dancers disturbs the flow of the dancing crowd. 1. Introduction The waltz was first recorded in writing by Salomo Jakob Wolf in 1797 [1]. He described a novel dance that stands out by its dynamic spinning motion (old German walzen translates to to spin). This dance has evolved and became part of Austrian culture in the form of the famous Viennese waltz. It is part of every ball in Vienna and is traditionally danced at midnight on New Years Eve. On first glance, rhythm and step sequence seem fairly simple, but quick pace and spin make it more difficult than it looks. In contrast to other ballroom dances like Boogie or Cha-Cha-Cha, which are part of the Latin American dances, couples do not remain in a small dancing area, but explore the entire dance-floor. Moreover, in comparison to dances like Quickstep, Rhumba or Tango, dancing couples are much faster and con- stantly spinning. Hence, beside correct execution of the steps themselves, dancers need to be careful to avoid collisions with other dancers and furnishing. This is often challenging as typical ballrooms are usually quite crowded. All potential obstacles have to be identified and collision avoiding/handling strategies have to be made in a fraction of seconds while constantly spinning at more than 180 per second. Dizziness CONTACT M. Bicher martin.bicher@tuwien.ac.at Institute of Information Systems Engineering, TU Wien, Vienna, Austria © 2020 The Author(s). Published by Informa UK Limited, trading as Taylor & Francis Group. This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial-NoDerivatives License (http://creativecommons.org/licenses/by-nc-nd/4.0/), which permits non-commercial re-use, distribution, and reproduction in any med- ium, provided the original work is properly cited, and is not altered, transformed, or built upon in any way. 170 M. BICHER ET AL. is a typical side effect that makes the dance even more challenging [1,2]. This is the reason why the basic steps and simple movements of Viennese waltz can lead to complete chaos on the dance-floor. Collisions between dancers, bottlenecks and clusters of several couples are common phenomena. Yet, many effects are only understood on the scope of individual dancers so far. Detailed analysis of the behaviour of the crowd would require different research techniques, for example modelling and simulation. 1.1. Research objective In our study, we aim to model and simulate a crowded ballroom with waltz dancers to gain deeper insights and identify influential parameters that lead to clusters and bottle- necks on the dance-floor. The outcomes of the study should lead to general advises, that help dancers to avoid getting stuck on the dance-floor. In order to do so, we need to research a typical example of a complex system [3]: We observe the macro-behaviour of a system, the dancing crowd, which is the direct result of the micro-behaviour of many interacting elements, the dancers. Hereby, we need to model research problems on different scales: Micro-Scale: How does a dancing couple perform a technically correct waltz-spin maintaining the right distance between each other? Meso-Scale: How does a dancing couple adapt its movements in response to other dancers, spatial limitation and its targeted dance direction? Macro-Scale: How does the crowd-density change on the dance-floor? Where and when can we detect clusters of dancers? We will present a simulation model that combines all three levels of abstraction and present our attempts of validation. Based on chosen case studies, we will draw conclu- sions on how clusters and collisions on the dance-floor can be avoided. 1.2. Novelty of the research Apart from medically oriented ones, there have been only a few studies on this topic. Most of them only deal with one of the three previously mentioned scales but do not investigate the entire system. Some studies investigate the physics of one single couple dancing (the Micro- Scale) in order to get insights into competitive ballroom dancing [4,5]. TheMeso-Scaleis analysed in [6], where trajectories of competitive ballroom dancers are tracked, yet the interactions between dancing couples are not analysed. In [7] the Macro-Scale was studied but only investigated from a psychological and not a spatial perspective. In summary, literature that considers all three scales and their effects on each other could not be found. 1.3. Modelling approach We decided to tackle the research question using a model, which is best interpreted as an agent-based modelling approach (ABM) with continuous time update. Creating a model combining all three levels of complexity, as explained in the previous section, excluded any other alternatives: MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 171 First of all, the system can hardly be modelled by one (or a few) aggregated ordinary differential equations, as the model needs to depict temporal and spatial changes. The corresponding extension, partial differential equations, are sometimes used for modelling of systems with interacting, moving particles such as traffic flow [8]or flocking [9]. Unfortunately, the macroscopic description form leads to a loss of information on the Micro- and the Meso-scale. This also holds true for all other modelling strategies in which the individuality of the dancers is lost, for example cellular automata. Consequently, a microscopic modelling concept is required. With respect to microscopic modelling techni- ques, classical discrete event simulation is not handy as well because the system involves continuous movement which cannot be neglected. Therefore, the strategy needs to be continuous as well. Continuous microscopic models are used exemplary for simulating particles, so-called molecular dynamics models. These models would dedicate the whole set of necessary (ordinary) differential equations for each particle (see [10]). Other approaches that model the individuals in a similar way are so-called gravity models [11] or the classical car-following models [12]. These approaches have the same foundation, and the behaviour and interaction are based on (differential) equations for each entity. Unfortunately, modelling every dancer or dancing couple solely via a continuous (system of) kinematic differential equations is not enough: the beat of the music poses for a discontinuity at which dancers need to perform a complicated decision-making process. Consequently, we chose to extend the ideas of the mentioned three modelling approaches by complicated time and state events, that involve classical decision-making processes. In [13] a pedestrian dynamics model is presented, that is based on the latter idea. As the author denotes the presented model by the term ‘agent-based’, we decidedtodosoaswell. The generic concept of an agent-based modelling [14] is a suitable framework to describe such highly complex models. It does neither involve any restrictions with respect to time dynamics nor with respect to the decision-making processes involved. Therefore, it is easily capable of combining continuous and discrete processes on a microscopic scale. Moreover, this flexible strategy is already famous for depicting group dynamic effects [15,16] and stands out with a very natural and picturesque description of the modelled system. One could call the model presented in this work a hybrid agent-based model: time-continuously updated with time- and state-events. 2. Model description To increase comprehensibility of the agent-based model (ABM) its description follows the ODD (Overview, Design Concepts, Details) protocol [17,18]. To keep the description of this quite complicated model as short and readable as possible, certain details are found in Appendix A for reproducibility purposes. 2.1. Overview 2.1.1. Purpose The ABM model simulates an arbitrary number of couples dancing in a ballroom to several beats of waltz music. The model allows to perform scenario tests like mixing dancers with different skill levels, comparing different distributions on the dance-floor as 172 M. BICHER ET AL. well as test different shapes of the ballroom. The results of these tests can furthermore be used to draw qualitative and partially quantitative conclusions about the dynamics of a waltz-dancing crowd. 2.1.2. Entities, state variables and scales In this model, each agent represents one dancing couple in reality. As we are interested in spatial interactions it is important to consider the physical properties of the dancing couples. Therefore, the physical appearance of a dancing couple is modelled by a solid cuboid with specific measurements and mass M, representing the two dancers common mass. The agent’s cuboid has two possible rotational axes R and R which represent the A B physical body centres of the two dancers and a distance d in between, as shown in Figure 1. The rotational axes are defined as two-dimensional coordinates relative to a global coordi- nate system. Note, that there is currently no specific purpose for the vertical dimension of the cuboid as all movement takes place in the plane. Yet, this feature allows 3D visualiza- tions and possible model extensions in the future. The movement of agents is, of course, inspired by waltz spins in reality: Agents may lock the position of one of the axes and apply a kinetic energy E to perform a spin around it. Periodically switching the rotation axis, the agent moves forward (see Figure 2). Hence, the coordinates of the rotational axes R ¼ R ðtÞ and R ¼ R ðtÞ change over time. A A B B The agents’ spinning behaviour depends on their individual parameters. These include, most importantly, a maximal kinetic energy E an agent can ‘produce’ for max spinning, as well as a default spinning angle α , which purpose will be explained later. def By default, agents apply their energy to spin in the clockwise direction, yet some of the agents are capable of spinning counterclockwise as well. We will describe this ability by the value of a boolean variable left. The parameters M; E ; α , left and d as well as the max def dimensions of the cuboid stay constant during simulation, while the kinetic energy applied for each spin will vary. The dance-floor represents the environment of the agents. As its appearance also influences the crowd’sbehaviour,wewilltestdifferent shapes of it. It is fully con- tinuous and every position on the dance-floor is reachable. In order to make agents ‘dance’ a music piece, represented as a sequence of beat-times (i.e. time instances of the first beat of every measure), is chosen, defining the length of the simulation run. As mentioned, the physical positions R ðtÞ and R ðtÞ of the agents are updated A B R R A B . . Figure 1. Abstraction of a dancing couple via a cuboid (dashed line) with two rotational axes R and R A B seen from above. The black and white circles should visualize the leading and following dance partner. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 173 R R A B . . Figure 2. An Illustration of the periodically switching rotation axis which makes the dancing couple both spin and move forward. The shape of the arc described by the dancing couple is solely defined by the achieved spinning-angle α. continuously (see Section 2.1.3) which makes the model both time- and spatially continuous. A full collection of model parameters and parameter values is found in Table 1. 2.1.3. Process overview and scheduling All agents in this model are treated simultaneously and update their state/position time- continuously according to classical frictionless equations of motion. Within this process, two types of events interfere with the movement of the agents: Table 1. Full list of all model parameters as well as used values or value-ranges, respectively. A value- range indicates, that this parameter is either varied in a case study or is defined by an agent-collective (see last column). Symbol Parameter Value/value-range Unit N total number of agents 2f1; .. . ; 120g 1 Case Study 1 N total number of agents per collective 2f1; .. . ; 120g 1 Case Study 3 t ; i ¼ 1; .. . ; time instants of the music beats see Section 4 s a outermost super-ellipse, first axle 2½17:7; 25:0 m Case Study 4 b outermost super-ellipse, second axle 2½12:5; 17:7 m Case Study 4 p initial distribution power-parameter 2½0:6; 1:4 1 Case Study 2 E maximum kinetic energy 2½228; 240 J see Table 2 max left counterclockwise ability ftrue; falseg boolean see Table 2 d distance between R and R 0:5m A B ð, ; , ; , Þ cube dimensions ð1:0; 2d; 1:7Þ m 1 2 3 m mass of a couple 150 kg Emax E default kinetic energy J α default angle 70 deg. st 1 measure nd 2 measure 174 M. BICHER ET AL. T T 1 2 a a 1 2 Figure 3. Hybrid automaton representation of the simulation model for the dancing couple, where the two events for the scheduling procedure are illustrated. The time event T represents the beat of the music and the state event T represents collisions between agents. Collisions (state event). Whenever the cuboids of two agents touch, a collision occurs. Inelastic collision handling is applied and both agents are deflected. This process is handled by classic collision equations and the event’s outcome can not be influenced by agents. Music-beat (time event). On the first beat of every three-four- or six-eight-time of the waltz music, the movement of all agents is entirely stopped and a new spin is initiated. Agents decide the amount of energy and the direction used in performing their next spin. During these events, agents are actively influencing the simulation via decision processes. Following the dynamic hybrid system (DHS) approach (see [19]), a graphical diagram of the described process can be depicted. Hereby, each agent, explained above, can be described by an (Alur-Henzinger) hybrid automaton (see Figure 3). Vertex (control mode) A describes the physical dynamics of the spinning movement and uses the state vector x ¼ðR ; R ; φ; ωÞ ; x y where ðR ; R Þ is the position of the current rotation axis of the agent, φ is the angle of the x y movement and ω the angular velocity. The mathematical description of the hybrid automaton requires event formulations, which are given by T ¼ Aa A and T ¼ Aa A: 1 1 2 2 Transition T corresponds to the time event Music-beat and T to the state event Collisions. 1 2 While T can be described as a typical time-event, using the classic DHS formulation, the state event cannot, as agent-agent-interactions are involved. Hereby, different to a classical DHS model, the guard and jump of T need to consider the states of all agents simulta- neously: The guard is represented on the interaction rules of the agents (collision detection) and the jump corresponds to the collision handling of the two involved agents. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 175 All three components, A, T and T are explained in Section 2.3 in more detail. 1 2 Note, that for implementation of such a system, time is usually discretized again. This is done by an ordinary differential equation solver (see Section 3.1). 2.2. Design concept 2.2.1. Basic principles From the modelling point of view, the most interesting feature of waltz dancing is the unusual way dancing couples move forward on the dance-floor: After spinning for one measure of waltz music, the orientation of the two partners is more-less switched and the couple, seen as a whole, has moved into the direction of their first step. After one more measure the orientation of the partners is almost back to the original setting, yet the couple itself has continued moving into the (more-less) same direction, compare Figure 2. Consequently, the movement of the whole couple emerges from the individual spin- ning movement of the partners and it is quite challenging to influence it: By adapting the speed of the rotation, the partners can not influence the velocity of the whole couple’s movement, but their orientation after the spin – they can ‘steer’ this way, but can neither break nor accelerate. Also, changing the spinning direction from the easier clockwise ‘Rechtswalzer’ to the much more sophisticated counterclockwise ‘Linkswalzer’ will not influence it. The actual velocity can only be changed by modifying the length of the individual steps, a parameter which is very difficult to change without getting problems with the synchronicity of both partners, e.g. physiological limits, stepping on partners feet, problems with keeping the angular momentum, etc. This observation founds the main concept of the simulation model: As agents are unable to change the distance of their rotation axes on run-time, they cannot influence how fast they actually proceed on the dance-floor. Yet, modifying their angular momen- tum and the direction of their spin, they can decide (up to some extents) where the whole couple is facing after the next measure. This type of steering is extremely unusual and becomes even more strange considering that aiming to the right or to the left with a clockwise spin makes a huge difference. 2.2.2. Emergence There are two elements in the present model which can be interpreted as sources of so- called emergent behaviour. On the lower level, the aforementioned trajectory of the agents can be considered as emerging from the individual spins of the agent. Even though the mathematical relation between individual angle and orientation of the agent after one spin is quite easy, the achieved individual angle depends on the kinetic energy applied and, in particular, the length of the music beat. For a music piece, especially the latter is usually not constant and the agents need to adjust the kinetic energy to compensate it (see Section 2.2.3). On the higher level, the interaction between the agents via collisions is a second source for emergent behaviour. 2.2.3. Objectives In the present model, agents only follow one objective: to stay on track. The tracks of the agents are defined as a set of steady curves that lead around the centre of the ballroom. In the real system, dancers are trying to follow more or less concentric paths around the 176 M. BICHER ET AL. centre of the ballroom in the counterclockwise direction. Based on observations, we found that neither ellipses nor rectangles sufficiently describe the shapes of these paths. The prior leaves the corners of the rectangular ballroom completely unused while the latter unrealistically over-exploits them. Hence, we decided to model the tracks as a family of fourth-order curves, so-called super-ellipses 4 4 x y þ ¼ c; c 2½0; 1: (1) 4 4 a b Hereby 2a and 2b correspond to the length and width of the ballroom. Moreover, any point ðx ; y Þ2 R on the dance-floor corresponds to one unique p p curve of this family. Hence, whenever an agent makes a decision about its spinning behaviour, it determines the super-ellipse that contains its current position, and aims to remain on this curve after the next spin. For detailed information on how this is done, the reader is referred to Section 2.3.4. Finally, the modelled agents do not aim to avoid collisions with other agents, though we would find such behaviour with dancers in reality. We found it feasible to neglect this feature because (1) modelling would be too complex, and (2) it can be substituted by a correct interpretation of the results: A collision of two agents in the model can be interpreted in the real system by either (a) at least one of the two couples is unable to finish its spin, (b) at least one of the two couples is performing some collision-avoidance move, or finally (c) an actual collision of two couples. 2.2.4. Sensing As mentioned in Section 2.1.3 an agent can only actively influence its motion at the music- beat events. Hereby the agent senses its velocity and its current position in the ballroom as well as thegraph of thecorrespondingsuper-ellipse(seeabove). The agent uses this information to calculate a target angle which it aims to spin in the course of the music- beat. Yet, in the course of this decision process, the agent is not aware of any other agents around. In case of a collision event, the agent is aware of the other agent and its motion, but cannot use this information for any future decisions. 2.2.5. Interaction In this model, the only form of interaction between agents is the collisions mentioned above (for more details on collisions, see 2.3.3). 2.2.6. Learning Due to occasional tempo changes in classical music, it is difficult for an agent to stay on track. In particular, the mapping α : energy applied ¼ E 7! αðE Þ¼ angle spinned kin kin depends on many additional parameters, namely (1) length of the current measure, (2) mass of the couple, (3) whether the spin is clockwise or counterclockwise, (4) distance between the rotation-axis and (5) whether a collision occurred or not. Although the physical relations between (1), (2), (3) and (4) are well known, (1) and (5) are difficult to predict. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 177 Table 2. The different skill level characterized through maximal kinetic energy and ability to switch spinning direction. The rightmost column contains the estimated maximal spinning angle per second. Collective name E (maximum energy) left (counterclockwise spin) est. angle range per second max perfect agents E ¼ 240 J true \ ½188; 188 max good agents E ¼ 240 J false \½0; 188 max slow agents E ¼ 228 J false \ ½0; 183 max So we decided to model this function as a black box combined with a learning process: At the beginning of the simulation, agents are unaware of the parameters of the black-box and perform their initial spin with a default energy E . Based on the outcome of each passed spin, they learn how much energy they have to apply to achieve a certain target angle. This process is described in 2.3.4 in more detail. 2.2.7. Stochasticity The model is a perfect example of an initial-value-stochastic microscopic model (com- pare with the classification proposed in [20]). While the initial positions and orientations of the agents are stochastic, the model dynamics, i.e. the decision rules of the agents as well as the kinematic movement, are deterministic. Hence, the same initial configuration of agents lead to the same simulation outcome. 2.2.8. Collectives In order to test the heterogeneity of the crowd more easily, we decided to summarize agents with equivalent traits to collectives. Hereby, we specifically refer to two para- meters, that define the agent’s ‘dancing skills’: The maximum kinetic energy E max applicable for spinning and the agent’s ability to perform counter-clockwise spins left (see Table 2). Note, that the estimated angle-range in this table is given for spins of precisely 1 s. Hence, a slow agent might get troubles performing 180 degree spins, if the measure of the waltz music is shorter than 1 s. To avoid misunderstanding, we want to clarify that these groupings are designed for analysis purposes only and do not have any impact on the model logic. 2.2.9. Observation First, the time-dependent positions R ðtÞ and R ðtÞ are recorded at each beat-event. A B After the simulation, this information is investigated for distribution of dancers with different skill levels. Second, the total number of agent-contacts, an agent has within the course of one simulation run, is tracked. We analyse the latter with respect to average collisions per agent and per skill-level. 2.3. Details 2.3.1. Initialization We initialize the model with a total number of N ; i 2f1; 2; 3g agents of each of the three collectives, explained in Section 2.2.8 and their skill parameters set as specified. We define N ¼ N þ N þ N as the total number of agents in the model. Furthermore, initial 1 2 3 positions for the agents are sampled. Given a uniformly drawn orientation in ½0; 2π of 178 M. BICHER ET AL. the agent, the centre point of the agent, i.e. ðR ð0Þþ R ð0ÞÞ, is set according to A B a specific distribution: Let X ¼ðX ; X Þ be a uniformly distributed point in 1 2 4 4 x x 1 2 ðx ; x Þ : þ  1 ; (2) 1 2 4 4 a b then the couple’s initial position is given by TðXÞ with p1 4 4 4 * * X X 1 1 TðXÞ¼ X  þ ; (3) 4 4 a a with a power-parameter p 2 R . Consequently, agents are distributed fully uniformly on the super-elliptical region for p ¼ 1. For p > 1 agents are more likely to be close to the centre of the room and for p< 1 vice versa. We will investigate the impact of this parameter in one of the case studies. 2.3.2. Submodel a – physical description of agent-movement The solely physical process is modelled frictionless and with full respect to energy con- servation. By applying angular velocities to the agents’ collision bodies, all angular and translational accelerations are and remain equivalent to zero at all times. In between any two events, the angular and translational velocity is constant and all positions change uniformly. Directly after the execution of a music-beat event (see 2.3.3) only the angular velocity around the rotation axis is non-zero, leading to a rotation with constant speed. In contrast, after the execution of a collision event also translational velocities can be non- zero (see 2.3.4) and the movement becomes a mixture of non-accelerated translation and rotation around a non-centred axis. 2.3.3. Submodel T – music beats Whenever a first beat of the three-four or six-eight time of the music occurs, i.e. at every time instant t , every agent starts a new spin around the switched spinning axis. Hereby, it decides on the energy it applies for the next spin and, in case the agent has an advanced skill level, whether it performs the spin clockwise or counterclockwise. This process can be split into the following steps: (1) The agent stops moving and all velocities are set to zero. (2) The rotation axis is switched. (3) Agents perceive their current position in the environment and find their target path. (4) Agents decide on how much energy they want/can apply for their next spin to reach the target and if they want/need to change spinning direction to succeed. (5) The kinetic energy is transformed to a new angular velocity and applied to the cuboid. A detailed description of steps 3–5 can be found in Appendix A. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 179 2.3.4. Submodel T – collisions In case the cuboids of two agents get in contact with each other, i.e. the surfaces of any two cuboids of agents intersect, inelastic collision-equations, i.e. collisions with restitu- tion coefficient equal to zero, are evaluated to solve it. This is feasible as collisions between humans usually lead to bending of body parts that absorb almost all the energy of the push. Equations for inelastic collisions in 2d and 3d are well known (e.g. see [21]) therefore they will not be described any further. 3. Simulation model analysis 3.1. Implementation In order to simulate the time-continuous elements of the model, the open-source 3D multi-body physics engine Bullet [22] was used. This powerful engine makes it possible to evaluate the collisions of hundreds of agents in reasonable time. As visual output is very important for face-validation purposes, we used Bullet’s interface provided by the graphics engine Panda3D [23]. Panda3D is primarily intended for the development of PC games (e.g. see Pirates of the Caribbean Online [24]) therefore the output of the simulation is nicely visualized in 3D and can be observed on run-time. In order to get a more realistic impression, the agents’ cube-shapes are replaced by more realistic 3D models of dancing couples, as seen in Figure 5. Both Panda3D and Bullet are freely available and are used as importable Python3 modules. A code for running the base model is available on GitHub https://github.com/MBicher/waltz_simulation. 3.2. Algorithms, complexity and numerics As the solution of the frictionless equations of motions does not involve any problems, Bullet’s internal semi–implicit Euler method is sufficient for integrating the occurring equations with a fixed time-step dt ¼ 0:025 ½s. Moreover, Bullet stands out with its fantastic features of collision detection and handling. Hereby it uses a Sweep and Prune algorithm to detect, and a variant of the projected Gauss–Seidel method to sequentially handle all occurring collisions [25]. Alike most simulations that involve the interaction of individuals, the computational and temporal complexity of the algorithm rises non-linearly when the number of agents is increased. Considering N as the total number of agents in the simulation, this procedure would lead to the computational effort of OðÞ N . Luckily, Bullet combines a Sweep and Prune algorithm and axis-aligned bounded boxes (AABB) to pre-process collisions and therefore reduces the computational effort to OðÞ N [26]. These observations do not take into account what happens in the course of the simulation on the dance-floor. The number of collisions that actually take place increases quadratic as long as the space of the dance-floor is kept constant (see next section). Consequently, the computational effort will increase by square as well. Apart from the long-term behaviour, our experiments revealed that up to an agent- count of about 50  70 the simulation time almost remains constant. Up to this number, the standard processes of Panda3D and Bullet, which occur during time-evolution, 180 M. BICHER ET AL. dominate the efforts for collision detection, evaluation and other iterative procedures. Thereafter, the computational effort increases as expected. 3.3. Value calibration Besides parameters actually measured from observing the real system (mass of a couple, size of the ballroom, . . .), Table 1 also shows one particular parameter, calibrated from simulation runs, namely the skill-dependent maximum kinetic energy E . This was max done by a simple parameter sweep comparing the applied kinetic energy with the maximum achievable angle. In other words, the second column of Table 2 was calibrated using the values in the fifth. 3.4. Verification To ensure that all implemented procedures work as modelled, a couple of verification scenarios has been set up and tested. These include: ● scenarios with only one dancing couple following a given trajectory (circle, ellipse, super-ellipse, sine-wave . . .), to make sure that the goal-seeking behaviour works, ● scenarios with two colliding dancing couples in different settings, to make sure that collisions are treated correctly, and ● scenarios with several dancing couples perfectly aligned in a circle, to investigate the propagation of numerical errors. 3.5. Validation Like most ABMs that include complex human behaviour, our dancing model is very difficult to validate. So far, we restricted our attempts to qualitative face validation together with dancing experts comparing open-source videos and pictures of famous balls in Vienna: Opernball, Zuckerbäckerball, Kaffeesiederball and, of course, the TU Ball, hosted annually by our university. Feedback of the face-validation was fully positive, as many group-dynamic effects, observed in reality, could also be reproduced in the virtual ballroom (see Section 5). In the future, we aim to perform additional validation attempts via trajectory-tracking devices. We plan to use motion capturing software to track a few number of dancers in a laboratory setting filmed from above, as well as GPS-based devices to track single trajectories of volunteers in a real ballroom situation. 4. Simulation scenarios and results In order to draw conclusions for the real system, we specified four case studies in which we varied specific model parameters and track the defined observables: the average number of collisions per agent and their time-dependent positions. Although the model only contains randomness within the initial conditions, the model is not only a stochastic one, but we expect the model dynamics to depend heavily MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 181 on the initial setting as well. Consequently, we need to perform a Monte Carlo simulation to make valid statements about the outcome: Each defined scenario is simulated 80 times with different random-number generation seeds. Around the gathered sample mean of the results of the simulation runs we will always state estimates for a confidence interval that contains the unknown real mean of the stochastic simulation. This estimate is based on the number of samples, the conver- gence speed of the sample mean, the sample variance of the results, and a Gaussian distribution assumption: With Φ being the cumulative probability function of the standard normal distribution and X and s being the sample mean and variance of N independent Monte Carlo samples, it is shown [27] that ! rffiffiffiffiffi 1  p s 1 N P X  μ Φ  p: (4) 2 N By adding these confidence interval estimates, we want to convince the reader that observed trends seen in the scenarios are not coincidental, exemplary due to fluctuations. We perform all scenarios for manually extracted measures of Waltz No. 2 by Dmitri Shostakovitch, which is 3:55 minutes long and has a total of 222 measures. The time- spans between consecutive measures vary between 0:83 s and 1:25 s and average in 0:97 s. Consequently, slow agents are not able to dance 180 degree spins all the time. 4.1. Case study 1: varying the number of dancing couples Our first subject of interest is the agent-density. Though it is clear that increasing the total number of dancers has a negative influence on the overall performance of the waltz, the rate of influence can only be quantified by simulation. Our simulation model reveals that the total number of agents N on the (same) dance- floor influences the average collision count per agent with approximately quadratic order. This can be justified analytically by thinking about the average number of collisions as the product of number of agents and the fraction of occupied space on the dance-floor. As the latter depends linearly on the number of agents N, the product grows approximately with N . Figure 4 shows a graph of the collision count for N from 30 to 130 with 80 simulation runs each. Hereby, only perfect agents have been considered (see Section 2.2.8). That means that all modelled agents are sufficiently skilled to perform spins with at least 180 degrees and are able to spin counterclockwise as well. The initial distribution is uni- form (p ¼ 1). The shaded area around the sample mean, calculated from 80 Monte Carlo runs, shows the mentioned estimated confidence interval that contains the real, unknown mean of the actual distribution of the simulation observable with 95% probability. Moreover, Figure 5 gives a snapshot of the simulation for N ¼ 100. The upper part of the figure shows the picturesque 3D-output of the simulation provided by Panda3D. In the lower part of the figure, the positions of the rotation axis R ðt Þ of all agents are B end depicted as red dots, and the corresponding trajectories R ðt Þ; i 2f0; .. . ; t g are given B i end as lines fading from black (very recent) to white (long time ago). For the regarded simulation setting with N ¼ 100, it shows quite smooth super-elliptic lines which 182 M. BICHER ET AL. empiric mean value of 80 simulation runs 95% confidence interval for the actual mean 30 40 50 60 70 80 90 100 110 120 130 140 150 number of agents N Figure 4. Average number of collisions per agent within 222 measures of waltz-music under variation of the total numbers of agents. The shaded area shows a confidence interval around the sample mean of the Monte Carlo analysis that contains the real, unknown mean with 95% probability (see (4)). Figure 5. Snapshot of the simulation with 100 perfect agents after spinning 222 measures. The upper image shows a snapshot of the visualization from Panda3D for case study 1, the bottom image shows the corresponding traces of the agents on the dance-floor, to be precise, the positions of R ðt Þ; i 2f1; .. . ; 222g. B i number of collisions per agent MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 183 indicate that agents hardly collided. It seems, the crowd circulated quite fluidly around the centre in this simulation run. 4.2. Case study 2: varying the initial distribution In the second case study, we investigate the impact of agents’ initial positions. Hereby the power-parameter p is varied (see Equation (3)) to find out, which initial distribution results in the minimal number of collisions. Parameter p is varied in steps of 0:1 between 0:3 and 1:3 for N ¼ 100 perfect agents. The mean value of collisions is drawn out of a Monte Carlo analysis with 80 simulation runs each. Results are shown in Figure 6. Although we would assume that a fully uniform distribution (p ¼ 1) of dancers would result in the most homogeneous crowd, the least number of collisions occurs with a parameter around p ¼ 0:6. Consequently, being placed slightly more densely to the border is advantageous with respect to avoiding conflicts. Figure 7 shows the difference between the optimal distribution suggested by the simulation compared with a fully uniform one. 4.3. Case study 3: varying the skill level Before performing tests with the heterogeneity among the crowd, the impact of reducing the skill level of all dancers, is investigated. It seems feasible that the skill to switch spin direction and to achieve larger maximum angles, should increase the couple’s ability to react to disturbances. Hence, we would expect the total number of collisions to increase when these skills are reduced. empiric mean value of 80 simulation runs 95% confidence interval for the actual mean 0.30.40.50.60.70.80.91.01.11.21.3 initial distribution coefficient c Figure 6. Impact of the initial distribution on the number of collisions. The average number of collisions per agent (gathered from 80 simulation runs each) is plotted against the power-parameter p. The shaded area shows a confidence interval around the sample mean of the Monte Carlo analysis that contains the real, unknown mean with 95% probability (see (4)). number of collisions per agent 184 M. BICHER ET AL. 10 10 5 5 0 0 −5 −5 −10 −10 −20 −10 010 20 −20 −10 010 20 ball-room length [m] ball-room length [m] Figure 7. Comparison of 3000 points drawn according to a fully uniform distribution (left, p ¼ 1) and the experimentally determined ‘optimal’ distribution for dancing (right, p ¼ 0:6). Table 3. Average number of collisions per agent for different skill levels and total agent counts. The expressions in the square brackets describe an estimated 95% confidence interval around the depicted sample mean of 80 simulation runs (see Equation (4)). Skill-level N ¼ 60 N ¼ 80 N ¼ 100 N ¼ 120 N perfect agents 48.96 [±3.40] 79.49 [±4.02] 121.42 [±6.48] 170.52 [±7.24] N good agents 58.20 [±4.76] 101.37 [±7.14] 178.31 [±9.88] 292.67 [±14.47] N slow agents 152.89 [±11.47] 237.50 [±11.74] 384.65 [±12.13] 562.70 [±16.59] mix of all three 119.81 [±6.23] 184.10 [±7.46] 304.68 [±12.18] 429.70 [±14.89] The simulation setting of Section 4.1 is slightly changed, (1) removing the agents’ ability to switch the spin direction (good agents) and moreover (2) reducing their maximum achievable kinetic energy (slow agents). Afterwards, heterogeneity is tested as we perform simulations with a homogeneous mix of agents with each of the three skill levels: A total of bN=3c perfect, bN=3c good, and bN=3c slow agents was homogeneously distributed (p ¼ 1) on the dance-floor for different total numbers N. The most important model outcomes are summarized in Table 3. The results clearly show that worsening the skill level leads to vastly more collisions occurring in the model. This effect is particularly significant as soon as the total number of agents on the dance-floor increases. Interestingly, the collision counts for the homo- geneously mixed scenarios are much more similar to the slow agent scenarios than the perfect or good agent ones. Apparently, this is a situation in which, proverbially, good apples, put into a bad situation, partially become bad apples. Anyway, in the course of this study we recognized that good agents and in particular slow agents tend to cluster in the centre of the ballroom due to their missing ability to steer to the right, i.e. towards the border of the ballroom. This situation is depicted via a Gaussian kernel-density estimation (KDE)of8000samples (80simula- tion runs with 100 agents each) in Figure 8. The upper heat-maps show the KDE for the agents’ positions at the beginning of each simulation, the lower show the KDE at the end. As a direct application of the previous finding, the influence of the initial positions of 50 slow agents together with 50 perfect agents is analysed. We assumed that, while on their way to the middle, slow agents would drag perfect agents with them leading to unneces- sarily many collisions. ball-room width [m] MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 185 good agents slow agents t =00 : 00 t =03 : 55 Figure 8. Kernel density estimation (KDE) of the positions of 100 agents in 80 simulation runs. The upper plots show the KDE at the start of the simulation while the lower plots depict the situation at the end of the music piece. The two left plots show the KDE for good agents, the two right plots for slow agents. Both collectives tend to cluster in the middle of the ballroom. Table 4. Average number of collisions compared for two different initial settings for a heterogeneous crowd. The expressions in the square brackets describe an estimated 95% confidence interval around the depicted sample mean of 80 simulation runs (see (4)). Scenario Avg. collisions per perfect agent Avg. collisions per slow agent Avg. collisions per agent ‘reference’ 325.05 [±13.68] 333.18 [±11.52] 329.12 [±10.55] ‘slow start at centre’ 74.26 [±4.51] 543.83 [±25.18] 309.04 [±12.82] In our ‘reference’ scenario the perfect and slow agents are mixed homogeneously and their position is initialized uniformly (p ¼ 1) on the dance-floor. In the alternative scenario (‘slow start at centre’), the initial positions are determined the same way, yet, afterwards, agents are resorted so that the initial positions of all slow agents are closer to the centre than those of all perfect agents. The results are depicted in Table 4 and partially confirm the assumptions: Leaving the outside tracks to perfect agents significantly reduces the number of collisions for them, but also vastly increases them for the slow agents. Interestingly, though, with respect to the overall average number of collisions per agent, both scenarios are almost equivalent. As the latter was initially defined as a ‘measure for fluidity’ in the ballroom, the alternative setting is not convincingly better than the reference. 4.4. Case study 4: varying the ballroom Run-time observations of simulation indicate that agents show immense problems staying on their track while spinning along the longer axis of the super-ellipses, as they need to achieve almost 180 degree per spin. Therefore, in this final case study, we analyse the impact of the dimensions of the ballroom on the collision count. 186 M. BICHER ET AL. Table 5. Average number of collisions compared for different ratios of the ballroom. The expressions in the square brackets describe an estimated 95% confidence interval around the depicted sample mean of 80 simulation runs (see (4)). Agents a : b ¼ 2 : 1 a : b ¼ 3 : 2 a : b ¼ 1 : 1 N ¼ 100 perfect agents 121.42 [±6.48] 104.29[±5.52] 101.75 [±5.84] N ¼ 100 slow agents 384.65 [±12.13] 258.70[±10.96] 202.63 [±8.14] a : b =1 :1 a : b =2 :1 Figure 9. Comparison of a quadratic and a rectangular ballroom. The colours of the displayed super- ellipses show how much angle an agent needs to spin per measure in order to stay on track. Keeping the total area of the ballroom constant, the ratio a : b of the two axis is varied and the total number of collisions is observed. Results are displayed in Table 5. The results show that the quadratic (a : b ¼ 1 : 1) ballroom helps especially slow agents as the regions in which they are forced to spin angles close to 180 degree are much shorter, as visualized in Figure 9. 5. Discussion Though especially quantitative validation is difficult, analysis of the model and the performed case studies (henceforth, short CS) leave enough scope for interesting con- clusions for the real system. Some of them confirm already known properties about the real system, while some others are, up to some extent, surprising. 5.1. Learn to spin counterclockwise In CS3 the feature of dynamically switching spinning direction, from clockwise to counterclockwise and vice versa, is an essential advantage compared to other dancers. It significantly decreases the total number of collisions due to the agents’ additional skill of ‘steering’ to the right. In reality, counterclockwise spinning is more difficult as it involves advanced foot- work. Being able to dynamically switch between clockwise and counterclockwise is even more complicated as dance-partners need to communicate properly within fractions of MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 187 a second. Yet, well-trained dancers are able to perform these kind of movements and are sometimes seen to use this skill to avoid conflicts with other dancers, exemplary by making a swift counterclockwise spin when the path on the left is blocked. Anyway, our model would suggest, that it is definitely worth training. 5.2. Spin fast and reduce the distance CS3 also showed that it is beneficial for an agent if it is capable to spin large angles, which is achieved by creating a high angular velocity v at the start of the spin. The model is ang based on standard laws of kinetics and kinematics. Therefore, three factors influence the maximum angular velocity an agent is able to apply for spinning: mass (moment of inertia), kinetic energy and distance between the two rotation axes R and R . Therefore, A B pffiffiffi reducing the maximum kinetic energy by a factor of x, increasing the mass by a factor x, pffiffiffi or increasing the distance between axis by a factor of x yields the same kinetic effects in the model. This fact is easily translated to the real system: if dance partners double the space between themselves, they need to develop four times more kinetic energy – i.e. apply more force, be quicker, make more accurate steps, etc. – to spin with the desired angle. Thus, in order to dance in an energy efficient way (which can be important to stand through a long glittering ball night) a couple needs to shift the two rotational axes as closely together as possible. 5.3. Do not distribute uniformly CS2 showed that a uniform initial distribution on the dance-floor is not the most beneficial with respect to collision avoidance, but the one depicted on the right side of Figure 7 is. Clearly, in reality, this distribution is not easy to achieve. In case all dancers would collaborate, they could follow the strategy: distribute evenly and then make one step towards the border of the room. 5.4. Stay away from the centre In most cases, the crowd is huge and does not collaborate. Hereby, CS3 gives some answers to where dancers should position themselves in order to prepare for a peaceful waltz. Experienced dancers should always try to find a spot closer to the borders, while novice dancers have two options: In case they want to be good team-players, they should place themselves closer to the centre in order to give experience dancers more freedom. By doing so, they pay the price of more ‘traffic’ around themselves. More selfish novices would pick a spot closer to the border, risking experienced dancers complaining about them. 5.5. The more quadratic, the better The model also provides insights for the organizers of balls. The results of CA4 indicated that organizers of balls should, in case of direct conflict of two equivalently sized 188 M. BICHER ET AL. rectangular rooms, prefer the more quadratic dance-floor. Also, splitting one large and long rectangular room into two more quadratic ones might be helpful. 5.6. Avoid overcrowded ballrooms Last but not least, it is generally a good idea to avoid overcrowded ballrooms. Although there is no real need for a simulation model to give this advise, our agent-based approach helps to understand the dynamics of overcrowding a little bit better. As seen in Figure 6 the collision count increases quadratically with the dancer-density, an effect which becomes even more prominent when not only perfectly trained, but also less skilled dancers are involved. 6. Conclusion In this work, we presented a highly complex individual-based model for ballroom dancing which stands out from other ABMs by several interesting features. Firstly, due to the interplay of music cadence and physics, an innovative hybrid update strategy was chosen. The dynamic update of a dancing couple agent combines aspects of time-/state- events at which typical decision-making processes are executed, and the classic differ- ential equations of motion which allows the agent to move. Secondly, planning the couple’s future track on the dance-floor requires solving a difficult unique control problem: The path of the whole couple needs to be steered by adjusting the angular velocity. The presented simulation model is the perfect environment for analysis of waltz dancing in a virtual laboratory and a lot more case studies with potentially interesting outcomes could be made. With respect to limitations, the model cannot and should not be used to observe the dancing behaviour of a single dancing couple without the context of the dancing crowd. For this, the individual behaviour of the dancers (i.e. steps sequence, body movement, . . .) is modelled in a too simplified way and would require adifferent modelling approach. Ripped costumes and even injuries are possible outcomes of colliding waltz couples and in most cases, the less skilled couple is blamed. Even though the results of the case studies show that, though the skill level of the dancers matters, this parameter is on the macro- scopic scope less influential than the density of dancers on the parquet. Hence, in order to avoid collisions with other dancers, the easiest solution is to wait until after midnight when the dance-floor is empty enough to dance peacefully – no matter the skill set. Acknowledgments We would like to give our special thanks to Prof. Dr Felix Breitenecker for many decades of not only academic, but also cultural education. His annual invitations to the traditional TU Ball inspired us to implement this model. The authors acknowledge TU Wien University Library for financial support through its Open Access Funding Programme. MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 189 Disclosure statement No potential conflict of interest was reported by the authors. References [1] S.J. Wolf, Beweis daß das Walzen eine Hauptquelle der Schwäche des Körpers und des Geistes unserer Generation sey, Hendel, Johann Christian, Halle, Germany, 1799. [2] S. Grobler, Thelifeand deathofthe pianowaltz, Ph.D. diss., University of Cincinnati, [3] S. Wolfram, Complex systems theory, Emerging syntheses in science: Proceedings of the founding workshops of the Santa Fe Institute, 1988, p. 182. [4] T. Shioya, Analysis of swing movement in ballroom dancing, Multidiscip. Digital Publ. Inst . Proc. 2 (2018), pp. 224. [5] K. Hachimura, K. Takashina, and M. Yoshimura, Analysis and evaluation of dancing movement based on LMA, in Robot and Human Interactive Communication, 2005. ROMAN 2005. IEEE International Workshop on, IEEE, Nashville, Tennessee, 2005, pp. 294–299. [6] P. Zaletel, G. Vuckovic, N. James, A. Rebula, and M. Zagorc, A time-motion analysis of ballroom dancers using an automatic tracking system, Kinesiologia Slovenica 16 (2010), pp. 46. [7] D. Gaertner, K. Clark, and M. Sergot, Ballroom etiquette: A case study for norm- governed multi-agent systems, in Coordination, Organizations, Institutions, and Norms in Agent Systems II, Springer, Hakodate, Japan, 2006, pp. 212–226. [8] T. Musha and H. Higuchi, Traffic current fluctuation and the burgers equation, Jpn. J. Appl. Phys. 17 (1978), pp. 811. doi:10.1143/JJAP.17.811. [9] F. Cucker and S. Smale, Emergent behavior in flocks, IEEE Trans. Automat. Contr. 52 (2007), pp. 852–862. doi:10.1109/TAC.2007.895842. [10] P. Deuflhard, J. Hermans, B. Leimkuhler, A.E. Mark, S. Reich, and R.D. Skeel, Computational Molecular Dynamics: Chellenges, Methods, Ideas, Springer, Berlin Heidelberg, Germany, 1999. [11] A. Sen and T.E. Smith, Gravity Models of Spatial Interaction Behavior, Springer, Heidelberg, Germany, 1995. [12] L.A. Pipes, An operational analysis of traffic dynamics, J. Appl. Phys. 24 (1953), pp. 274–281. doi:10.1063/1.1721265. [13] T. Schelhorn, D. O’Sullivan, M. Haklay, and M. Thurstain-Goodwin, STREETS: An agent-based pedestrian model, in CUPUM’99 Computer in urban planning and urban management on the edge of the millennium: Proceedings of the 6th International Conference, 01, Franco Angeli, Venice, Italy, 1999, pp. 1–12. [14] C.M. Macal and M.J. North, Tutorial on Agent-Based Modeling and Simulation Part 2: How to Model with Agents, in Proceedings of the 2006 Winter Simulation Conference, IEEE, Monterey, California, 2006, pp. 73–83. Available at http://www.informs-sim.org/wsc06pa pers/008.pdf. [15] E. Bonabeau, From Classical Models of Morphogenesis to Agent-Based Models of Pattern Formation, Artif. Life. 3 (1997), pp. 191–211. doi:10.1162/artl.1997.3.3.191 [16] R.M. Axelrod, The Evolution of Cooperation, Basic Books, New York, 1984. [17] V. Grimm, U. Berger, F. Bastiansen, S. Eliassen, V. Ginot, J. Giske, J. Goss-Custard, T. Grand, S.K. Heinz, G. Huse, A. Huth, J.U. Jepsen, C. Jrgensen, W.M. Mooij, B. Mller, G. Peer, C. Piou, S.F. Railsback, A.M. Robbins, M.M. Robbins, E. Rossmanith, N. Rger, E. Strand, S. Souissi, R.A. Stillman, R. Vab, U. Visser, and D.L. DeAngelis, A standard protocol for describing individual-based and agent-based models, Ecol. Modell. 198 (2006), pp. 115–126. Available at https://linkinghub.elsevier.com/retrieve/pii/S0304380006002043. 190 M. BICHER ET AL. [18] V. Grimm, U. Berger, D.L. DeAngelis, J.G. Polhill, J. Giske, and S.F. Railsback, The ODD protocol: A review and first update, Ecol. Modell. 221 (2010), pp. 2760–2768. Available at https://linkinghub.elsevier.com/retrieve/pii/S030438001000414X. [19] T.A. Henzinger, The Theory of Hybrid Automata,in Proceedings - Symposium on Logic in Computer Science, IEEE,Piscataway, NJ, United States, 1996, pp. 278–292. DOI:10.1017/ s0963180100007039 [20] M. Bicher, Classification of Microscopic Models with respect to Aggregated System Behaviour, Diss., TU Wien, Vienna, Austria, 2017. [21] D.T. Greenwood, Principles of Dynamics, Prentice-Hall, Englewood Cliffs, NJ, 1988. [22] E. Coumans, Bullet physics engine, Open Source Software, Available at http://bulletphysics. org1 (2010), p. 3. [23] M. Goslin and M.R. Mine, The panda3d graphics engine, Computer 37 (2004), pp. 112–114. doi:10.1109/MC.2004.180. [24] C. Lang, Panda3D 1.7 Game Developer’s Cookbook, Packt Publishing Ltd, Lincoln Road Olton Birmingham, UK, 2011. [25] E. Coumans, Bullet 2.80 physics sdk manual, Available at http://www.cs.uu.nl/docs/vakken/ mgp/assignment/Bullet (2012) [26] D.J. Tracy, S.R. Buss, and B.M. Woods, Efficient large-scale sweep and prune methods with AABB insertion and removal, in 2009 IEEE Virtual Reality Conference, IEEE, Lafayette, Louisiana, 2009, pp. 191–198. [27] S. Ross, Simulation, Knovel Library, Academic Press, Amsterdam, 2012. Appendix A. Detailed Description of Submodel T A.1. Finding the new target (3). As mentioned, every agent first needs to perceive its position to find out which of the super-ellipses (1) currently describes its track. We consider the i-th beat of the music piece and the new rotation axis of the agent is R with R ðt Þ¼ðxðt Þ; yðt ÞÞ (anything stated henceforth extends analogously A A i i i 4 4 for R ). Moreover, the agent’s current track calculates to ðxðt Þ=aÞ þðyðt Þ=bÞ ¼ c . B i i i Based on these observations, the agent solves the root-finding problem 4 4 ðxðt Þþ d cosðαÞÞ ðyðt Þþ d sinðαÞÞ i i c ¼ þ 4 4 a b with respect to the unknown α. It results in the two intersection points of the super-ellipse and a circle with radius d centred in R ðt Þ. Hereby, it gives the agent two possible target points to A i aim for to stay on its current track. As it is convention to move counterclockwise on Viennese ballrooms, the one with the larger polar angle, with respect to the centre of the ballroom, is the correct solution. This solution, henceforth denoted as R ðt Þ, forms the agent’s new target B iþ1 point. In order to calculate how much energy needs to be applied, the agent calculates the new target angle α :¼ \ R ðt ÞR ðt ÞR ðt Þ (measured clockwise). p B i A i B iþ1 A.2. Calculate the new energy and direction (4). As mentioned in 2.2.6, agents memorize their previously performed spins to learn from them. With E ; j 2f1; ... ; ig denoting the previously used kinetic energies for spinning, α ; j 2 j j f1; .. . ; ig denoting the successfully danced angles, and t ; j 2f0; .. . ; ig standing for the lengths of all previously danced measures, the agent forecasts the necessary energy E to achieve angle α p p with the next spin. This is done assuming that the square root of the kinetic energy, the resulting angle and the length of the measure are in a linear relationship: MATHEMATICAL AND COMPUTER MODELLING OF DYNAMICAL SYSTEMS 191 Figure A1.. Diagram for deciding the amount of energy E and direction d applied for the i þ 1-st iþ1 iþ1 spin starting at time t . We assume that this spin has rotation axis R with current position R ðt Þ. i A A i A spin around R is calculated analogously. () pffiffiffiffi E ðt  t Þ j j j1 c :¼ Median : j 2f1; .. . ; ig; α  0 : (A1) Using the median eliminates outliers resulting from spins interrupted by collisions. Moreover, we use the median of the last four measured lengths dt :¼ Median ðt  t Þ; j 2fi  4; .. . ; ig (A2) m j j1 to get a feasible guess for the length of the next beat. The values c; dt and α are then used to m p calculate the target energy E using E :¼ c : (A3) dt Furthermore, the target energy E is compared to the maximum energy E the agent is able to p max use for spinning. If E  E , the target energy is applied as E to the new rotational axis in the p max iþ1 form of torque. Otherwise, a couple of escape strategies are attempted which are summarized in a decision-tree structure in Figure A1. 192 M. BICHER ET AL. In case every strategy fails, the agent aims for a default angle α of 70 degree. This angle was def chosen based on observations and experience to avoid a spinning deadlock. Variations of this parameter in the final simulations showed that its precise value has no influence on the simulation outcomes. Finally, in case i ¼ 1 when no previous learning data is available, the agent simply applies the default angle E into the default spinning direction (clockwise). def A.3. Energy to angular velocity (5). Finally, the applied kinetic energy is transformed into angular velocity via sffiffiffiffiffiffiffiffiffiffi rffiffiffiffiffiffiffiffiffiffi jv j 2E 1 4E 1 tang kin kin v ¼ ¼ ¼ : ang d m d m d spin Hereby, we assume that the mass of the spinning object m is equivalent to the mass of one of the spin dance partners and therefore approximately half of the mass of the couple. The sign of the expression is determined by the direction of the spin.

Journal

Mathematical and Computer Modelling of Dynamical Systems – Taylor & Francis

Published: Mar 3, 2020

Keywords: Agent-based simulation; ballroom dancing; complex system; Viennese Waltz

References

  • The life and death of the piano waltz
  • Complex systems theory
  • Analysis of swing movement in ballroom dancing
  • Analysis and evaluation of dancing movement based on LMA
  • A time-motion analysis of ballroom dancers using an automatic tracking system
  • Ballroom etiquette: A case study for norm- governed multi-agent systems
  • Traffic current fluctuation and the burgers equation
  • Emergent behavior in flocks
  • An operational analysis of traffic dynamics
  • STREETS: An agent-based pedestrian model
  • Tutorial on Agent-Based Modeling and Simulation Part 2: How to Model with Agents
  • From Classical Models of Morphogenesis to Agent-Based Models of Pattern Formation
  • A standard protocol for describing individual-based and agent-based models
  • The ODD protocol: A review and first update
  • The Theory of Hybrid Automata
  • Classification of Microscopic Models with respect to Aggregated System Behaviour
  • Bullet physics engine
  • The panda3d graphics engine
  • Bullet 2.80 physics sdk manual
  • Efficient large-scale sweep and prune methods with AABB insertion and removal