Design and Validation of Cyber-Physical Systems Through Co-Simulation: The Voronoi Tessellation Use Case

This paper reports on the use of co-simulation techniques to build prototypes of co-operative autonomous robotic cyber-physical systems. Designing such systems involves a mission-specific planner algorithm, a control algorithm to drive an agent performing its task; and the plant model to simulate the agent dynamics. An application aimed at positioning a swarm of unmanned aerial vehicles (drones) in a bounded area, exploiting a Voronoi tessellation algorithm developed in this work, is taken as a case study. The paper shows how co-simulation allows testing the complex system at the design phase using models created with different languages and tools. The paper then reports on how the adopted co-simulation platform enables control parameters calibration, by exploiting design space exploration technology. The INTO-CPS co-simulation platform, compliant with the Functional Mock-up Interface standard to exchange dynamic simulation models using various languages, was used in this work. The different software modules were written in Modelica, C, and Python. In particular, the latter was used to implement an original variant of the Voronoi algorithm to tesselate a convex polygonal region, by means of dummy points added at appropriate positions outside the bounding polygon. A key contribution of this case study is that it demonstrates how an accurate simulation of a cooperative drone swarm requires modeling the physical plant together with the high-level coordination algorithm. The coupling of co-simulation and design space exploration has been demonstrated to support control parameter calibration to optimize energy consumption and convergence time to the target positions of the drone swarm. From a practical point of view, this makes it possible to test the ability of the swarm to self-deploy in space in order to achieve optimal detection coverage and allow unmanned aerial vehicles in a swarm to coordinate with each other.


I. INTRODUCTION
Cyber-Physical Systems (CPS) are a large class of systems characterized by a complex interplay of hardware and software components, often operating in a largely unpredictable environment.CPSs are usually designed with model-based development (MBD) techniques [1]: design begins with a The associate editor coordinating the review of this manuscript and approving it for publication was Nasim Ullah .
high-level model of the system, which is validated, corrected, and refined in several cycles, until a model is obtained that developers judge fit to be taken as a production blueprint.
System validation is used to check the compliance of each system element and of the overall system with its purpose and functions.Validation activities are planned and carried out throughout the development process.
Validating CPS models is complex, as they must reflect the interactions of several subsystems and different aspects of each subsystem: For example, the physical aspects of an actuator include mechanics and electromagnetism, mechanics in turn require modeling such concepts as torque, inertia, viscosity, and so on.Such complexity makes it desirable, and often necessary, to rely on several modeling languages and tools.Actually, most tools offer a large set of specialized libraries that span a wide range of application fields, nevertheless it is often the case that more specific tools are preferred.Also, CPSs may be developed in parallel by various teams, possibly belonging to different organizations, which may have expertise with different modelling tools.
In cyber-physical systems, a co-simulation approach can be adopted, which allows different system elements to be simulated each by dedicated simulators, with a co-simulation master orchestrating their execution.The Functional Mockup Interface [2] is a public-domain standard that defines a container and an interface to exchange dynamic simulation models using various languages.FMI is supported by more than one hundred modelling tools [3].
In this work, co-simulation is used to study the behaviour of a swarm of unmanned aerial vehicles (UAV) (or drones) that must cover uniformly a given area, which is often useful in many applications, such as search and rescue or environmental monitoring.Due to their maneuvrability, UAVs are convenient tools for monitoring areas difficult to reach for land vehicles [4], [5].In particular, among different types of UAV's, quadcopters have become the most popular for their flexibility, due to very low moments of inertia, greater stability, hovering capability, as well as lower takeoff requirements [6], [7].
Designing a swarm of drones involves three levels of control: (i) a high-level, mission specific planner algorithm defining the final position of each drone; (ii) an intermediatelevel control algorithm to drive each drone along its prescribed trajectory; and (iii) a low-level plant model to simulate the drone dynamics.
In this schema, the crucial part is the design of the intermediate-level algorithm.In fact, in a space-coverage task the final position of the drones is entirely determined by the geometry of the area to be covered and by the number of drones, so that designing the high-level control reduces to choosing an appropriate geometric algorithm.In turn, simulating the individual drones relies on standard dynamics models whose parameters depend on the physical characteristics of the drone.
The intermediate-level control, instead, has a more complex task, involving the optimization of various contrasting figures of merit (e.g,, speed and fuel consumption), whose relative weight depends on the particular application with its operational constraints.
This paper discusses how co-simulation and design space exploration support the development of this kind of systems.In particular, the approach based on centroidal Voronoi tessellation has been used for the high-level control.
Its Python implementation was validated with a prototype simulation based on a simplified physical model.The validated implementation was then integrated in a full-scale co-simulation, along with the two other control algorithms.Co-simulation was performed on the INTO-CPS framework [8].
The Voronoi algorithm was first verified using typical values for controller parameters.Then the INTO-CPS design space exploration tool was used to calibrate the nonlinear controller parameters to optimize power consumption and responsiveness.
The main contribution of the paper consists in (i) demonstrating the application of co-simulation to the problem of quadcopter control for space coverage; (ii) showing how design-space exploration coupled to co-simulation supports control parameter calibration to optimize energy consumption and convergence time.A technique to extend the tessellation algorithm to a bounded area is also introduced.
This paper is organized as follows: Section II reports related work on cyber-physical systems design, including modeling/simulation tools and co-simulation.Section III introduces background notions (quadcopters, Voronoi tessellation, and co-simulation), for the selected case study.Section IV presents the technique developed for centroidal tessellation of a bounded area.Section V reports the cosimulation architecture, the simulation scenario, and the preliminary results from co-simulation.Section VI shows and discusses the calibration of controller parameters for energy consumption and convergence time.Finally, Section VII concludes the paper and reports on future work.

II. RELATED WORK
System models are typically built with block-based languages, which describe a system as an assembly of functional blocks, each representing a possibly complex mathematical operation, interconnected by data flows.For example, MATLAB/Simulink [9], [10] is a commercial environment for the model-based development of cyber-physical systems.It provides a graphical model editor and functions to generate demonstrative prototypes and to analyse relevant model properties.
Some works define ad-hoc solutions for the integration of different simulators, dedicated to different system components, e.g., the hardware and the software parts.For example, in [11] Simulink and a simulator of logic specifications (the pvsio interpreter of the Prototype Verification System [12]) have been integrated to simulate and formally verify properties of a pacemaker.
Farhat et al. [13], instead, use Matlab and Simpack to compare the performance of mechatronically-driven railway vehicles' guidance and steering to that of a conventional vehicle.A non-linear Simpack vehicle model with the specific mechatronic actuation and sensing, and a simplified linear control model in Matlab/Simulink are co-simulated.
Reinhart and Weissenberger [14] address multibody simulation in the context of numerical control machines.
The multibody model considers flexible machine structural components, guideways, feed drive dynamics and axis controllers together with the motion trajectory generation of the control.
In [7], Simulink/Gazebo tools are used for the validation of an approach to quadcopter control where wind disturbance is modeled by unknown exogenous inputs, exploiting the Robot Operating System (ROS) middleware in the simulations.In [15], the architecture ROS/Gazebo has been extended with the possibility of simulation of co-operative UAVs.
Co-simulation is a technique in which simulators of heterogeneous models are executed simultaneously while exchanging data during run-time.The FMI (functional mockup interface) [2] is a standardized interface for model exchange and co-simulation of dynamic models across different modeling and simulation tools.
Co-simulation has been applied extensively in different application fields.In [16], a co-simulation approach is used for a brushless motor, with the electrical and mechanichal parts modeled in Simulink and the feedback linearization control modeled with a function written in C.
In [17], a co-simulation open source software for operation and planning studies of distributed energy resources has been presented.The software allows users to perform largescale high-fidelity simulations for bulk power system (BPS) planning and operation.
Papers [18] and [19] show how human performance models can be incorporated into models of CPSs.In [18], a train driver model is coupled to models of the rolling stock and of the movement authority; in [19] human-machine interfaces of an integrated clinical environment are considered.
In [20], co-simulation was applied to analyse Model Predictive Control systems for autonomous driving.This requires an accurate analysis of the interplay among three main components: the plant, the model predictive control algorithm, and the processor where the algorithm is executed.Co-simulation of the three components was used to determine if the controller running on the chosen hardware meets the time requirements determined by the response time of the plant.Satisfactory tradeoffs between algorithm complexity and processor performance could be studied.
In [21], the co-simulation-based framework Vico for marine crane onboard operations is presented.Simulation and analysis of sub-systems is performed with different software tools and the framework enables the digitalization of marine operations.
In [22], co-simulation is applied to co-operative mobile robots.Simulations for multi-robot systems are executed by independent tools, such as Gazebo for physics and mobility, ROS2 for software development, and ns-3 for communications and the networking infrastructure [23].The framework integrates such simulators and allows running experiments that combine all the involved robotic systems keeping the synchronization time between the simulators consistent.
In [24], an approach to the formal verification of a variant of a well-known consensus protocol of UAVs is presented, using a co-simulation framework for system validation.The present paper builds on [24], discussing control parameter calibration by design space exploration, using a more complex case study.
In [25], Cho et al. present an advanced co-simulation platform that concurrently performs UAV simulation and wireless network simulation.This platform is taylored to the simulation of centrally controlled UAV swarms, whereas the approach of the present paper can be applied to distributed control applications, although only a central control case study is discussed.Also, the implementation reported by Cho et al. does not rely on the FMI standard.
An extensive survey on co-simulation has been published by Gomes et al. [26], providing definitions of the fundamental concepts and a taxonomy of the literature based on the discrete events and continuous time computational models.
A paper by Fitzgerald et al. [27] introduces foundational and process management issues for the design of CPSs and cites several methods, languages, and tools, using a twowheel self-balanced personal transporter as an example.
An approach to CPS design extending the software-defined networking (SDN) concept is the software-defined cyberphysical systems (SD-CPS) framework [28], wherein a workflow of microservices (web services), coordinated by an SDN controller, realises the support for CPS operation.An SD-CPS federated controller deployment is proposed, which comprises several controllers from different domains or organizations at the edge, in a case study with embedded mobile systems in connected cars, edge cloud nodes and servers.Another approach inspired by SDN and service-oriented architecture is Cyber Physical System as a Software-defined Service (CPSS) [29], wherein the CPS control platform is divided into the Infrastructure-as-a-Service, Network-asa-Service, and Aplication-as-a-Service layers.CPS design methods as exemplified by the above references can be smoothly integrated with co-simulation approaches such as the one discussed in the present paper, due to their reliance on a distributed software architecture.

III. BACKGROUND
This section introduces notions used in the rest of the paper: quadcopter unmanned vehicles, Voronoi tessellation and the INTO-CPS framework for co-simulation.

A. QUADCOPTERS
Quadcopters have become very popular UAVs for their flexibility.However, combined control of rotational and translational motions is required to obtain accurate path tracking and autonomous flight capacity, which results in highly nonlinear modeling [30].
The dynamics of a quadcopter can be described with reference to Fig. 1, where the craft is modeled as a crossshaped rigid frame supporting four rotors, each turning at speed ω i (i = 1, . . ., 4), whose thrust and torque determine 1066 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
the velocity and attitude of the craft.A moving reference frame F ′ has its origin at the center of mass, with the arms lying on axes x ′ and y ′ , whereas F, with axes x, y, and z, is a fixed reference frame.The attitude of the craft is defined by its Euler angles ψ, φ, and θ (yaw, roll, and pitch), which describe the rotation that brings the axes of F ′ to coincide with those of F after an appropriate translation.
Equations ( 1) and ( 2) below express the dynamics for position and attitude, respectively, as functions of the rotational speeds of each rotor: where g is the gravity acceleration, ψ d is the desired yaw angle, C φ , C θ , C ψ , and C z are physical parameters.The quadcopter controller must compute the rotor speeds that will drive it to the desired position (x d , y d , z d ) with the desired yaw ψ d .This can be achieved with the design shown in Fig. 2, where the position controller commands the attitude controller, which computes the rotor speeds.The latter are fed to each rotor's low-level controller, included in the plant block, that translates the desired rotor speeds into current levels for the motors.
An attitude controller can be defined by equations of the following form: In the above equations, λ P and λ A are the controller gain parameters, which affect the responsiveness of the controller to tracking errors.Specifically, λ A is related to the desired responsiveness of attitude stabilization, while λ P is related to that of position control.Since the attitude control loop is internal to the position control one, the frequencies of the eigenvalues associated with the attitude dynamics should be higher than those of the position dynamics.As a rule of thumb, the two gain parameters should ideally satisfy the condition λ A /λ P > 10.In this work, the physical parameters of Erle-copters (Fig. 3) are used in the simulations [7].

B. VORONOI DIAGRAMS
Given a set S of k points in a plane , a Voronoi diagram [31] is a partition of into k Voronoi regions (or cells) V i , each associated with a distinct point (seed or generator) s i ∈ S, and consisting of all the points in the plane closer to s i than to any other point s j , with j ̸ = i, in S, i.e. [32]: In general, the Voronoi region V i of a point s i ∈ S = {s 1 , . . ., s k } can be constructed as follows: (i) the perpendicular bisectors l ij of the segments between s i and the other points of S are drawn; (ii) if ij is the half plane bounded by l ij and containing s i , V i is the intersection of the ij 's. Figure 4 shows the construction of the cell for seed s 1 in a set of six seeds.The Voronoi diagram for set S is then the union of the Voronoi regions, as shown in Figure 5.Each region is bounded by a finite number of segments or half lines, called ridges.In the figure, dashed lines represent infinite ridges delimiting unbounded regions.The points where two or more ridges meet are called vertices.

1) AN IMPLEMENTATION OF THE VORONOI TESSELLATION ALGORITHM
The implementation of the Voronoi algorithm used in this work relies on the Python class SciPy.Spatial.Voronoi [33].In order to show how a Voronoi diagram is represented in the Python code, we use the simple example of Figures 4 and 5.The Voronoi class is initialized with an array containing the tessellation seeds, as shown in Listing 1.
The class constructor computes the tessellation vertices and ridges, stored in class attributes, as shown in Listing 2. Each vertex is identified by its index in the array vor.vertices, ranging from 0 to 5 in the example.The vertex indices are used in the list vor.ridge_vertices to define the ridges.Each ridge is identified by its extremes: if an extreme is a vertex it is referred to by its index in the vor.vertices array, while an extreme at infinite distance is represented by −1.For example, the pair [−1, 0] represents a half line originating from vertex 0, i.e., point [8.98840764, 4.78471338], whereas [0, 1] represents the ridge between vertices 0 and 1, i.e., points [8.98840764, 4.78471338] and [3.36882475, 5.11037441].
The above data can be compared with Figure 5, produced by the program.The blue dots are the seeds, the solid lines are the ridges of the only finite region, and dashed lines represent (semi)infinite ridges.

2) CENTROID COMPUTATION
In many applications of drone swarms, it is convenient to place each drone in the geometric center, or centroid of the respective cell.A Voronoi tessellation whose generating points are the centroids of each cell is called a centroidal Voronoi tessellation [31].A centroidal Voronoi tessellation can be obtained by Lloyd's method [34]: Given an initial set of k seeds {s 0,i } k i=1 ; 1) construct the initial Voronoi tessellation {V 0,i } k i=1 for the seeds {s 0,i } k i=1 ; 2) at each step n+1, compute the centroids of the Voronoi regions {V n,i } k i=1 found at step n; these centroids become the new set of seeds {s n+1,i } k i=1 ; 3) If this new set of seeds meets some convergence criterion, terminate; otherwise, return to step 2.
Step 2 above is performed by computing the centroid coordinates (C x , C y ) with the formulae for a polygon with n 1068 VOLUME 12, 2024 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
vertices [35], [36]: This method produces a set of drone placements that optimize the average distance of each drone from the points of its assigned region (centroid property) and the spacing between drones (Voronoi tessellation property).

C. CO-SIMULATION
Co-simulation techniques allow a complex and heterogeneous simulated system to be split into simple sub-systems.The simulation of each sub-system can be executed by a dedicated simulator, possibly supporting a different modeling language.
The Functional Mockup Interface (FMI) [37] is a standard for co-simulation.Each sub-system is exported as a Functional Mockup Unit (FMU) and a co-simulation is the execution of a set of FMUs that exchange data and signals under the co-ordination of a master program.Figure 8 in Section V shows an example of co-simulation architecture, where the Co-simulation Orchestration Engine (COE) is the master program.Various frameworks for FMI-compliant cosimulation have been developed, including INTO-CPS [8], [38] and CoSim20 [39].In particular, the INTO-CPS framework is used in this work.
The INTO-CPS toolchain includes the Design Space Exploration (DSE) tool [40], used to find optimal combinations of parameter values within specified ranges or discrete sets searched by a DSE engine that executes the cosimulations for each combination.In particular, the Pareto method is applied to ranking the results of the DSE.The set of combinations with the highest rank (rank 1) represents the best compromises of a pair of objective functions.The set contains all combinations of parameter values where it is not possible to find other combinations which improve both objective functions.

IV. CENTROIDAL TESSELLATION OF A BOUNDED REGION
The Python implementation of the Voronoi algorithm described in Section III-B1 was taken as a starting point for the simulation of a drone swarm.That implementation computes a tessellation over an unbounded plane, while the envisioned application involves covering a bounded region.A workaround was found to tessellate a convex polygonal region: dummy points at appropriate positions outside the bounding polygon are added to the original set of seeds, then the Python implementation is applied to the extended set to tessellate the infinite plane containing the polygon.More precisely, for each seed one dummy point is placed symmetrically to the seed with respect to each side of the polygon.so that, for each seed, as many symmetric points as the number of sides of the polygon are created.
As an example, Figure 6 shows the generation of the tessellation of a square with four seeds, with sides parallel to the cartesian axes.Figure 6  This modified algorithm was used to perform Step 1 of Lloyd's method (Section III-B2).In this way, the behaviour of a swarm of drones that must optimally cover an area bounded by a polygon can be described as follows: 1) The initial positions of the drones are the seeds of the initial tessellation; 2) the (modified) Voronoi tessellation for the current seeds is computed; 3) for each cell, the centroid is computed; 4) each drone is moved to the centroid of the respective cell; 5) if the new position of all drones is within a given (small) distance from the previous one, stop; otherwise, repeat Step 2 with the new positions as current seeds.
This behavioral model was implemented in Python, obtaining a high-level simulator of a space coverage strategy.Figure 7 shows the initial, an intermediate, and the final configurations of a ten-drone swarm deployed over a hexagonal area.In the pictures, produced by the Python simulation, the blue dots are the initial drone position at each cycle of the algorithm, and the red ones are the corresponding centroids.
The model abstracts from the physical model of the drones, thus ignoring finite speeds, control-loop delays, and communication delays.A more accurate model, implemented by co-simulation, is discussed in the next section.

V. CO-SIMULATION OF A UAV SWARM
Co-simulation was used to study the behavior of a drone swarm with the task of reaching an optimal coverage of a plane surface bounded by a polygon.To this purpose, the planner algorithm, implemented in Python, uses the centroidal Voronoi tessellation described above to assign, at each simulation step, a destination point to each drone, starting from arbitrary initial positions.The related FMU has been generated using UniFMU [41].The drone's destination is the command input to the intermediate-level drone controller, implemented in C, that computes the rotor speeds according to (3), (4), and (5) in Section III-A.A low-level plant model, implemented in Modelica [42], simulates the drone according to (1) and (2).More details on the control FMU and plant FMUs are discussed in [24].
Figure 8 shows the architecture of the co-simulation of a swarm of drones, where all the FMUs are connected to the master algorithm of the co-simulation, the Co-simulation Orchestration Engine (COE), while Figure 9 shows the logical connections between the planner and one drone.The logical connections can be explicitly provided to the COE through a graphical user interface made available by the INTO-CPS Association.
The co-simulation has two synchronization granularities: the fixed co-simulation step T at which the plant and the intermediate-level controller interact, and a step T ′ at which the planner FMU is executed.Step T ′ is an integral multiple of T , i.e., T ′ = τ T , with τ a positive integer.

A. SIMULATION SCENARIO
In the simulated scenario, ten drones must be placed over a hexagonal area of 1600 m 2 .The co-simulation step T (0.05 s) was chosen to be of the same order as the period between successive updates of the outputs computed by actual controllers in drone systems.
Figure 10 shows the initial displacement of the drones as shown by the graphical user interface created with the technology discussed in [43] that allows users to track the movement of the drones through the evolution of the co-simulation and Figure 10b shows the final placement achieved by the Voronoi tessellation.
In the co-simulation, the system evolves as follows: 1) the Voronoi tessellation is computed every τ cosimulation steps;   2) each drone executes the control algorithm every co-simulation step T, moving towards the centroid of its region for τ co-simulation steps, then a new tessellation is computed; 3) the simulation ends at time T s , This value is chosen large enough to let the drones reach the desired final ''honeycomb'' formation.Informally, when the drones reach such a configuration, they keep executing the algorithm, which does not further change their positions (a formal definition of convergence is given in subsection VI).

B. PRELIMINARY TESTS FOR λ A AND λ P
A preliminary series of co-simulations was performed to find starting points and ranges for the parameters to be optimised by exploring the design space.The following results were obtained using a trial and error approach.For example, with T = 0.05 s, the co-simulation reveals divergent behaviors, even for values of λ A and λ P that are close to satisfying the ideal condition λ A /λ P > 10 (Section III-A).Specifically, Figures 11 and 12 show that the coordinates x and y begin to oscillate around the desired values and finally diverge.In addition, it was found that the acceptable ranges of values for controller parameters depend also on the controller update rate T .Indeed, in the instances of Figures 13 and 14, the x and y coordinates oscillate and diverge, even though the ideal ratio between the two control gain parameters is satisfied.In general, it was found that, for T = 0.05 s, the controller parameters must satisfy the two constraints λ A ≤ 9.0 and λ P ≤ 0.9. Figure 15 shows the result of a co-simulation where these last constraints are met, and as expected the x and y coordinates move towards the values of the assigned centroid, without oscillations; the final position is reached in approximately two minutes of simulated time.

VI. DESIGN SPACE EXPLORATION
The DSE tool was then used to find optimal values for controller parameters λ A and λ P , according to the criteria of minimizing (i) time t c to attain convergence, and (ii) energy consumption.
Convergence was defined to be attained at the first instant since when each drone remains within a circle of radius   R from its target position for a minimum interval of δ seconds.The chosen values for R and δ were 0.01 m and 5 s.Convergence times were computed off-line from simulation traces.
Energy consumption for a motor is computed with the following formula: where E is the total energy until t c , T is the simulation step, and P i is the motor power at step i: where K f is a constant of the motor and ω d is the difference between motor speed and the hovering speed ω h , a characteristic of the drone.We observe that energy consumption is affected by the update frequency of the desired co-ordinates (τ T ).
Table 1 shows the parameter values defining the design space considered in the co-simulations, with the constraint λ A ≥ 5λ P , resulting in twenty-six combinations.The simulation time T s is 300 s for all simulation runs.  2  and 3, by energy consumption in Tables 4 and 5, and by Pareto rank in Tables 6 and 7. Tables 2 and 3 confirm that larger values of τ cause longer convergence times.As an example, consider the case of λ A = 7 and λ P = 0.9.For τ = 6 (Table 2) the update frequency of desired co-ordinates is 0, 3 s and the convergence time is 155.25 s.For τ = 20, (Table 3) the update frequency of desired co-ordinates is 1 s and the convergence time is 244, 85 s.In particular, for τ = 20 (Table 3), the values of t c show that the maximun simulated time of 300 s was reached   Tables 4 and 5 show that energy consumption increases with higher values of λ A and λ P (with few exceptions).It may be observed that the impact of λ P on energy consumption is greater than the impact of λ A and that τ has a relatively low impact on energy consumption.The ordered data series is split into two tables for readability.

Results shown ordered by convergence time in Tables
Finally, Tables 6 and 7 compare the optimality of the various combinations of values for τ , λ A and λ P by respective Pareto rank.Here, too, the ordered data series is split into two tables.
In Figure 17, it is possible to observe the combinations with the best rank (rank 1), shown with green points joined by a green line.The combinations that provide the best tradeoff   The controller defined by (3) is nonlinear, which makes the characterization of optimal combinations of parameter values more difficult than in the case of linear systems.For this reason, simulation-based design space exploration is an effective tool to deal with this problem.

VII. CONCLUSION
This paper proposes an approach based on co-simulation to design and validate co-operative robotic systems since the early phases of the system design.In particular, this work shows how design space exploration can be successfully applied to the expensive task of parameter calibration for the control algorithms of a co-operative UAV-based cyberphysical system.
One of the advantages of this approach is the easy interaction between multidisciplinary teams (e.g., information engineers, automation engineers, etc), due to the support for heterogeneous languages and simulation tools afforded by co-simulation.For example, in the use case discussed in this work, Python was used for the co-ordination algorithm, C for the drone controller and Modelica for the drone plant.Another advantage is found in model reuse, as models already available can be imported as black boxes in the cosimulation framework, provided that the original modeling tools support the FMI standard.Moreover, this modularity eases system design refinement.The model presented in this paper, for example, could be refined by including network communications and delays or extended by considering alternative system architectures.
In general, this approach makes it possible to develop a consistent set of models for a given system at different levels of abstraction, enabling a wide spectrum of analyses to be performed while minimizing development effort.It may be added that the same approach can be used in more complex applications.For example, instead of a swarm of drones of the same type, the deployment of vehicles of different kinds, possibly both air-and land-faring, could be considered.
As further work, it is planned to consider (i) different system architectures for the co-ordination of vehicles; and (ii) study UAVs equipped with further sensors, for example cameras, for activities such as disaster management and recovery.

FIGURE 4 .
FIGURE 4. Construction of the Voronoi region for seed s 1 .

FIGURE 5 .
FIGURE 5. Voronoi diagram for the six points of Fig. 4.

LISTING 1 .LISTING 2 .
Set of seeds.Data structures for example of Fig.5.
(a) shows the seeds and the dummy points corresponding to the leftmost seed and Figure 6(b) shows the tessellation computed by the Voronoi algorithm.The tessellation of the bounding square is obviously obtained by deleting all cells outside the square, as shown in Figure 6(c).

FIGURE 6 .
FIGURE 6.How to build a Voronoi diagram inside a polygon.

FIGURE 7 .
FIGURE 7. How to build a Voronoi diagram inside a polygon.

FIGURE 8 .
FIGURE 8. Co-simulation architecture of a swarm.

FIGURE 10 .
FIGURE 10.Drone positions shown by the GUI.
are the ones closest to the bottom left corner of the graph, i.e, τ = 6, λ P = 0.8, and λ A = 7 or λ A = 8.

TABLE 3 .
Ordered by convergence time (τ =20).before reaching convergence.Further, for a given value of τ , parameter λ P determines the speed of convergence.