Distributed area coverage control with imprecise robot localization: Simulation and experimental studies

This article examines the area coverage problem for a network of mobile robots with imprecise agents' localization. Each robot has uniform radial sensing ability, governed by first order kinodynamics. The convex-space is partitioned based on the Guaranteed Voronoi (GV) principle and each robot's area of responsibility corresponds to its GV-cell, bounded by hyperbolic arcs. The proposed control law is distributed, demanding the positioning information about its GV-Delaunay neighbors. Simulation and experimental studies are offered to highlight the efficiency of the proposed control law.


I. INTRODUCTION
A REA coverage of a planar region by a swarm of robotic agents is an active field of research with great interest and potential. The task involves the dispersal inside a region of interest of autonomous, sensor-equipped robotic agents in order to achieve sensor coverage of that region. Distributed control algorithms are a popular choice for the task because of their computational efficiency compared to centralized controllers, their ability to adapt in real-time to changes in the robot swarm and their reliance only on local information. Thus they can be implemented on robots with low computational power and low power radio transceivers, since each agent is required to communicate only with its neighbors. The downside of distributed control schemes is that they lead to local optima of their objective function because of their lack of information on the state of all agents.
Area coverage can be either static [1]- [3] in which case the agents converge to some static final configuration or sweep [4]- [6] in which the agents constantly move in order to satisfy a constantly changing coverage objective.
Other factors that need to be taken into account are the type of the region surveyed [7], [8], the sensing pattern and performance of the sensors [9]- [12], or the particular dynamics of the agents [13].
There have been varied approaches to area coverage such as distributed optimization [14], [15], model predictive control [16], [17], game theory [18] and optimal control [19]. This  All localization systems have an inherent uncertainty in their measurements, thus there has been significant work on taking this uncertainty into account either by data fusion [20], safe trajectory planning [21] or spatial probability [22]. Our approach creates a suitable space partioning scheme in lieu of the agents' positioning uncertainty and a distributed gradientbased control law.
This article examines the problem of persistent area coverage of a convex planar region by a network of homogeneous agents equipped with isotropic sensors. In order to take into account the localization uncertainty, each agent is assumed to lie somewhere inside a circular positioning uncertainty region and Guaranteed Voronoi partitioning [23] of the region of interest is constructed from these uncertain regions. Subsequently, a distributed control law is derived based on the aforementioned partitioning so that a coverage objective increases monotonously. Because of the complexity of that law, a suboptimal one is proposed and both laws are compared through simulations. Additionally, the suboptimal control law is evaluated by two experiments using differential drive robots. This article is an extension of [24] which did not contain simulation results from the optimal control law and lacked the experimental evaluation of the suboptimal control law.
The layout of the article is as follows. The required mathematical preliminaries are presented in Section II. The Voronoi and the Guaranteed Voronoi space partitions are examined in Section III. Section IV contains the definition of the coverage objective and the derivation of the gradient ascent based control law that maximizes it. Additionally it contains the proposed suboptimal control law. Simulation studies are presented in Section V in order to compare and evaluate both control laws. Section VI contains experimental results from the implementation of the suboptimal control law and is followed by concluding remarks.

II. MATHEMATICAL PRELIMINARIES
Let us assume a compact convex region Ω ∈ R 2 under surveillance and a network of n identical dimensionless mobile agents (nodes). Each agent's dynamics are governed bẏ where u i is the control input for each agent and I n = {1, 2, . . . , n}. The agents' positions are not known precisely and thus in the general case each agent has a positioning uncertainty of arXiv:1612.04704v1 [cs.SY] 14 Dec 2016 different measure. We assume that an upper bound r u i for the positioning uncertainty of each agent is known and thus its center lies within a disk where · is the Euclidean metric and q i the agent's position as reported by its localization equipment. All agents are equipped with identical omnidirectional sensors with circular sensing footprints and as such the sensed region of each agent can be defined as where r s is the common range of the sensors. We also define for each agent a 'Guaranteed Sensed Region' (GSR) as the region that is sensed by the agent for all of its possible positions inside C u i . The GSRs are defined as and since C u i and C s i are disks, the previous expression can be further simplified into provided that r s ≥ r u i . When an agent's position is known precisely, i.e. r u i = 0, its GSR is equivalent to its sensed region, whereas when an agent's positioning uncertainty is too large compared to its sensing capabilities, i.e. r u i > r s , then C gs i = / 0.

III. SPACE PARTITIONING
Most distributed area coverage schemes assign an area of responsibility to each agent based on information from its neighboring agents. Then each agent is responsible for maximizing the coverage solely on its own responsibility region. The most common partitioning scheme is the Voronoi partitioning, described in Section III-A. However in this article, because of the localization uncertainty of the agents, the guaranteed Voronoi diagram described in Section III-B is used for the assignment of areas of responsibility.

A. Voronoi Diagram
The Voronoi diagram for a set of points Q = {q 1 , q 2 . . . q n }, q i ∈ Ω is defined as follows Each V i ⊂ Ω is called the Voronoi cell of node i. A Voronoi diagram of 6 points constrained in a subset Ω of R 2 is shown in Figure 1 (left), where the cell boundaries are shown in blue.
It is important to note that the Voronoi diagram is a complete tessellation of Ω, that is i∈I n V i = Ω.
Another useful property of Voronoi diagrams is their duality to Delaunay triangulations. Let us define the Delaunay neighbors of a point q i as When constructing the Voronoi cell of a point q i , only its Delaunay neighbors need to be considered, thus simplifying the definition of its Voronoi cell into Using this property, a Voronoi diagram can be constructed by first finding the Delaunay triangulation of the input points and then for each cell V i by finding the intersection of the halfplanes defined by point q i and all points in N i .

B. Guaranteed Voronoi Diagram
The Guaranteed Voronoi (GV) diagram is defined for a set of uncertain regions D = {D 1 , D 2 , . . . , D n }, D i ⊂ Ω. Each uncertain region D i contains all the possible positions of a point q i ∈ R 2 whose localization cannot be measured precisely. Thus the GV diagram is defined as An equivalent definition of the GV diagram is as follows This definition is analogous to the construction of the classic Voronoi cell by the halfplane intersection method described previously, with the difference that the regions being intersected are no longer halfplanes. As such, each cell Vi g contains the points on the plane that are closer to the corresponding point q i of region D i for all possible configurations of the uncertain points. A Guaranteed Voronoi diagram of 6 regions constrained in a subset Ω of R 2 is shown in Figure 1 (right), where similarly the cell boundaries are shown in blue.

Remark 1.
In contrast to the classic Voronoi diagram, the Guaranteed Voronoi diagram does not constitute a tessellation of Ω, since i∈I n V g i ⊂ Ω. We define a neutral region O so that O ∪ i∈I n V g i = Ω. Similarly to the Voronoi diagram, we can define the Guaranteed Delaunay neighbors N g i of a region D i so that in order to construct V g i , only the regions in N g i need be considered. Thus the Guaranteed Delaunay neighbors are defined as Having defined the Guaranteed Delaunay neighbors, the definition of the GV diagram can be simplified into It is important to note that when two or more uncertain regions overlap, none of them are assigned GV cells. In addition, if the uncertain regions degenerate into points, the GV diagram converges to the classic Voronoi diagram.

C. Guaranteed Voronoi Diagram of Disks
The GV diagram is examined in the case when the uncertain regions D i are disks C u i as defined in Equation (2). In that case it has been shown that H i j are regions on the plane bounded by hyperbolic branches. For any two nodes i and j, ∂ H i j and ∂ H ji are the two branches of a hyperbola with foci at q i and q j and semi-major axis equal to the sum of the radii of C u i and C u j . The parametric equation of these hyperbola branches can be found in the Appendix.
The dependence of the GV cells of two disks on the distance of their centers d i j = d(q i , q j ) is as follows. When the disks C u i ,C u j , overlap both of the cells V g i ,V g j are empty as shown in Fig. 2 (a). If the disks represent the possible positions of the agents, this is an undesirable configuration as it can lead to collisions. When the disks C u i ,C u j are outside tangent (i.e., d(q i , q j ) = r u i + r u j ), the resulting cells are rays starting from the centers of the disks q i , q j and extending along the direction of q i q j as seen in Fig 2 (b). When the disks are disjoint, the GV cells are bounded by the two branches of a hyperbola. As d i j increases further, the eccentricity of the hyperbola increases and the distance of the disk centers from the hyperbola vertices (the points of the hyperbola closest to its center) increases as seen in Fig 2 (c), (d). The distance between the hyperbola's vertices remains constant at r i + r j as d i j increases.
The GV cells of two disks also depend on the sum of their radii r i + r j as seen in Figure 3. As r i + r j decreases, the eccentricity of the hyperbola increases and the result on the cells is the same as if the distance of the disks' centers increased, as explained previously. When r i + r j = 0, the GV cells are the classic Voronoi cells.
It has been also shown that the Guaranteed Delaunay neighbors N g i are a subset of the classic Delaunay neighbors of the uncertain disks' centers. The two sets N i and N g i are equal when the GV diagram is constructed for the whole plane, however when the GV diagram is constrained within a subset Ω of the plane, N g i can be a subset of N i . Such a case is shown in Figure 1 where nodes 1 and 5 are Delaunay neighbors but not Guaranteed Delaunay neighbors. This way, the GV diagram of disks can be constructed in a manner similar to the classic Voronoi diagram. First the Delaunay triangulation of the disks' centers is constructed and then each cell V g i is constructed by intersecting the regions bound by hyperbolic arcs H i j ∀ j ∈ N g i . This process is shown in Figure 4. It is useful to define Guaranteed-Sensing Voronoi (GSV) cells as which are the parts of the GV cells that are guaranteed to be sensed by their respective nodes. It can be shown that since GV cells are disjoint, GSV cells are also disjoint regions.

IV. PROBLEM STATEMENT -DISTRIBUTED CONTROL LAW
By taking into account the GV partitioning of the space, the following area coverage criterion is constructed where the function φ : R 2 → R + is related to the a priori knowledge of importance of a point q ∈ Ω indicating the probability of an event to take place at q in a coverage scenario. This criterion expresses the area that is guaranteed to be closest to and covered by the nodes. A distributed control law is designed in order to maximize the coverage criterion (9) while taking into account the nodes' guaranteed sensing pattern (3) and the GV partitioning (7). It is shown that this control law leads to a monotonic increase of the coverage criterion.
Theorem 1. For a network of mobile nodes with uncertain localization as in (2), uniform range-limited radial sensing performance as in (3), governed by the individual agent's kinodynamics described in (1), and the GV-partitioning of Ω defined in (7), the coordination scheme where n i is the outward unit normal on ∂V gs i and α i a positive constant, maximizes the performance criterion (9) along the nodes' trajectories in a monotonic manner, leading in a locally area-optimal configuration of the network.
Proof. In order to guarantee monotonous increase of the coverage criterion (9), we evaluate its time derivative By using the gradient-based control law monotonous increase of H is achieved. We now evaluate the partial derivative of H with respect to q i as Considering the second summation term, infinitesimal motion of q i may only affect ∂V gs j at ∂V gs where O is the neutral region defined in Remark 1, since both hyperbola branches are affected by alteration of one of the foci. In addition, only the Guaranteed Delaunay neighbors (6) of i are considered in the summation, as a major property of GVpartitioning. Therefore, the former expression can be written via the generalized Leibniz integral rule [25] (by converting surface integrals to line ones) as where υ i i , υ i j stand for the transpose Jacobian matrices with respect to q i of the points q ∈ ∂V gs i , q ∈ ∂V gs j , respectively, i.e.
The boundary ∂Ṽ gs i can be decomposed in disjoint sets as These sets represent the parts of ∂V gs i that lie on the boundary of Ω, the boundary of the node's guaranteed sensing region, and the boundary of the unassigned neutral region of Ω, respectively. This decomposition is shown in Figure 5 with ∂V gs i ∩ ∂ Ω in green, ∂V gs i ∩ ∂C gs i in solid red and ∂V gs i ∩ ∂ O in solid blue. Hence ∂ H ∂ q i can be written as It is apparent that υ i i = 0 at q ∈ ∂V gs i ∩ Ω since all q ∈ ∂ Ω remain unaltered by infinitesimal motions of q i , assuming no alteration of the environment. Considering the second integral, for any point q ∈ ∂Ṽ gs i ∩ ∂C gs i it holds that υ i i (q) = I 2 , where I stands for the identity matrix, since they translate along the direction of motion of q i at the same rate. As far as the sets ∂V gs i ∩ ∂ O and ∂V gs j ∩ ∂ O, j ∈ N g i are concerned, they can be merged in pairs via utilization of the left and right hyperbolic branches, as introduced in GV-partitioning (7), as follows However, for any two Delaunay neighbors i, j it can be observed that • The hyperbolic branches ∂ H i j and ∂ H ji are symmetric with respect to the perpendicular bisector of q i , q j , and are governed by the same set of parametric equations (left and right branch). • The vectors n j are mirrored images of the corresponding n i with respect to the perpendicular bisector of q i , q j .
Taking into account the above, ∂ H ∂ q i can be written as The unit normal vectors n i , n j and Jacobian matrices υ i i , υ i j in the second part of (17) can be evaluated via utilization of the parametric representations of the sets over which the integration takes place. These sets are parts of the left (∂ H i j ) and right (∂ H ji ) branch of the hyperbola that assigns the space among two arbitrary nodes i, j. The second term of the sum in (17) requires knowledge of all the nodes in j∈N g i N g j ⊆ I n and so the control law is distributed.
The decomposition of the set ∂V gs i ∩ ∂ O of node i into mutually disjoint hyperbolic arcs ∂V gs i ∩ ∂ H i j and ∂V gs i ∩ ∂ H ik is shown in Figure 6. Figure 6 also shows the other domains of integration used in (17)   The proposed law (10) leads to a gradient flow of H along the nodes trajectories, while H increases monotonically, since Remark 2. It is not immediately clear from the control law (10) the information node i requires in order to implement it. The integrals over ∂V gs i ∩ ∂C gs i and ∂V gs i ∩ ∂ H i j require knowledge of the Guaranteed Delaunay neighbors N g i in order to compute V gs i and H i j . However the integral over ∂V gs j ∩ ∂ H ji requires knowledge of the whole cell V gs j of node i and as such requires informations form the Guaranteed Delaunay neighbors N g j of node j. Thus node i requires information from all nodes inÑ which are the Guaranteed Delaunay neighbors of all the Guaranteed Delaunay neighbors of node i.
A. Constraining Nodes Inside Ω When using this control law, it is possible that a portion of some node's positioning uncertainty region lies outside the region Ω, i.e. that C u i ∩ Ω = C u i for some node i. As a result, it would be possible for that node to be outside the region Ω given that it may be anywhere within C u i . In order to prevent this situation, instead of Ω, a subset Ω s i ⊆ Ω will be used, in the general case different for each node since their uncertainty radii are allowed to differ. This subset is calculated as the Minkowski difference of Ω with the disk C (r u i ) = {q ∈ Ω : q ≤ r u i }.
Thus by constraining the center of each uncertainty disk q i inside Ω s i , it is guaranteed that no part of the uncertainty disk will ever be outside Ω. This is achieved by projecting the velocity control input u i on ∂ Ω s i . This is only needed if q i ∈ ∂ Ω s i and the contorl input u i points towards the outside of Ω s i . As a result, the control law becomes where ε is an infinitesimally small positive constant.

B. Suboptimal Control Law
Because of the high complexity of the integrals over H i j and H ji of the control law (10), the large amount of information required in order to implement it, as described in Remark 2, as well as the problems described in Section IV-A, a simplified version of the control law is proposed. leads the nodes to suboptimal trajectories.
The control law enhancement described in Section IV-A is used in conjunction with this simplified control law.

V. SIMULATION STUDIES
Simulations have been conducted in order to evaluate the efficiency of the control law as well as compare the complete (10) with the simplified (19) law. All nodes have common positioning uncertainty regions C i , ∀q i ∈ Ω as well as sensing performance C gs i . The a priori importance of points inside the region was φ (q) = 1, ∀q ∈ Ω. In this case, the maximum possible value for H is where A (·) is the area function, given that n guaranteed sensing disks can be packed inside Ω. In all simulations, the control law extension as described in Section IV-A is used.

A. Case Study I
In this case study the convergence speed between the optimal and suboptimal laws is examined. A network of three nodes in a large region is simulated. Figure 7 shows the initial and final states of the network, with the uncertainty disks shown in black, the guaranteed sensing disks in red and the GV cells in blue. In Figure 8 the trajectories of the nodes when using the optimal control law are shown in red and when using the suboptimal control law in blue, while the initial positions are represented by circles and the final by squares. It is observed that the final configurations in both cases are very similar since the average distance between nodes' final states is 0.57% of the region's diameter. This is due to the low number of nodes in the network and is a result that can not be generalized, as seen in the following case study. Figure 9 shows the evolution of the coverage objective H with time, again shown in red for the optimal control law and in blue for the suboptimal. As it is expected, the optimal law converges faster than the suboptimal. The final value of the objective function is the same for both control laws in this case since the network reaches a globally optimal state in both cases.

B. Case Study II
In this case study network of ten nodes is examined. Figure  10 shows the initial and final states of the network, with the uncertainty disks shown in black, the guaranteed sensing disks in red and the GV cells in blue. In Figure 11 the trajectories of the nodes when using the optimal control law are shown in red and when using the suboptimal control law in blue, while the initial positions are represented by circles and the final by squares. It is observed that the final configurations in both cases seem very similar at first glance, however the average distance between nodes' final states is 10.12% of the region's  diameter. Careful examination of the trajectories shows that not all nodes starting from the same initial position converge to the same final position, explaining the large average distance. . Figure 12 shows the evolution of the coverage objective H with time, again shown in red for the optimal control law and in blue for the suboptimal. As it is expected, the optimal law converges faster than the suboptimal. The final value of the objective function is the same for both control laws in this case because of the similarity of their final configurations. However, since both control laws converge to local maxima, there are no guarantees that one or the other will always achieve better performance.
This case study highlights the need for the control enhancement in Section IV-A, especially in the case of the optimal control law. It can be seen in Figure 11 that a node moving under the optimal law (red trajectory) was stopped from moving outside the region on the lower left corner.

VI. EXPERIMENTAL STUDIES
In order to evaluate the simplified control law (19), two experiments were planned and conducted. The experiments   11. The trajectories of the nodes using the optimal (red) and suboptimal (blue) control laws.
consist of 3 robots, an external pose tracking system and a computer, communicating with the robots and the pose tracking system via a WiFi router. The robots used are the differential-drive AmigoBots by ActiveMedia Robotics, which provide a high level command set, such as setting their translational and rotational velocities directly. These robots are also equipped with encoders, used for self-estimating their pose, but due to the drifting error, an external tracking system was implemented. Although the external pose estimation system is used for the control law implementation, its performance is compared with that of the encoders. Since the robots are not equipped with actual sensors, circular sensing patters with a radius r s = 0.3 m were assumed in the control law implementation. The pose estimation system consists of a camera and an ODROID-XU4 octa core computer, running an ArUco based pose estimation server. ArUco [26] is a minimal library for Augmented Reality applications and provides pose estimation of some predefined sets of fiducial markers. Hende, by attaching a fiducial marker on top of each robot, the ArUco library provides an estimation of each robot's pose.
In order to estimate the positioning uncertainty of the experimental setup, a fiducial marker was placed on the centroid and on each of the vertices of the region Ω. A 3D triangular mesh was created from these measurements as shown in Figure 13. The maximum uncertainty value from that mesh r u max = 0.032 m was used for the experiments. The control law was implemented algorithmically as a loop with a period T s = 0.1 sec. After the end of each iteration the computer sends velocity commands to the robots. At first, the robots' current positions are used to create the Guaranteed Voronoi diagram of the uncertain disks and calculate the control inputs according to the simplified control law (19). Subsequently, a target point q t i ∈ R 2 is created for each robot based on its current position q i and its control inputũ i . Once the target point for each robot has been found, and given each robot's position q i and orientation θ i , the translational v i and rotational ω i velocity control inputs are sent to each robot.
where dθ i = ∠(q t i − q i ) − θ i and v max , ω max are constraints on the maximum translational and rotational velocity respectively of the Amigobot robots. Once all robots are within a predefined distance d t = 0.02 m of their respective target points, the Guaranteed Voronoi diagram and the control law are calculated again and new target points are produced.
The region shown in all figures regarding the experiments is Ω s ⊆ Ω i.e. the region inside which the centers of the uncertain disks are constrained, as explained in Section IV-A.
Both experiments are evaluated against simulations with the same initial robot positions, positioning uncertainty and sensing performance.

A. Experiment I
The initial and final positions of the robots for both the experiment and the simulation are shown in Figures 14 and  15 respectively with the uncertainty disks shown in black, the guaranteed sensing disks in red and the GV cells in blue. Additionally, photos of the experiment initial and final configurations can be seen on Figure 16. It can be seen that in the final configuration all guaranteed sensing disks are contained within their respective GV cells in both the experiment and the simulation, thus leading the network to a globally optimal configuration. By comparing the final configurations, as well as the robot trajectories for the experiment and simulation shown in Figure 17, it can be seen that they differ significantly. The mean distance between the nodes final positions in the experiment and simulation is 8.68% of the diameter of Ω s . Because in this current case there are multiple globally optimal configurations the network can possibly reach and since the implemented control law differs from the theoretical one, it is natural that the final configurations between the experiment and the simulation differ. Nevertheless, a globally optimal configuration was reached in both cases. The coverage objective H did not increase monotonously due to the implemented control law. It did increase however from 83.7% to 99.9% of its maximum possible value. The positioning data from the ArUco library and the robot encoders are compared in Figure  19 [Left]. As it is expected, the further a robot moves and the more it rotates, the larger the positioning error of the encoders grows. The mean error between the ArUco and encoder final robot positions is 4.04% of the diameter of Ω s . Figure 19 [Right] shows the trajectories of the target points used in the control law implementation.

B. Experiment II
In this experiment, half the diameter of the robots was added to their positioning uncertainty, so that each robot was completely contained within its positioning uncertainty disk      simulation. By comparing the final configurations, as well as the robot trajectories for the experiment and simulation shown in Figure 23, it can be seen that they do not differ significantly. The mean distance between the nodes final positions in the experiment and simulation is 1.54% of the diameter of Ω s . Because of the increased positioning uncertainty in this case, which results in smaller GV cells, the number of local maxima of the objective function is significantly smaller than in the previous experiment. Thus it is quite possible that, despite their differences, the theoretical and implemented control laws converge to the same final configuration. The coverage objective H did not increase monotonously due to the implemented control law. It did increase however from 38.8% to 98.3% of its maximum possible value. The positioning data from the ArUco library and the robot encoders are compared in Figure  25 [Left]. As it is expected, the further a robot moves and the more it rotates, the larger the positioning error of the encoders grows. The mean error between the ArUco and encoder final robot positions is 1.54% of the diameter of Ω s . Figure 25 [Right] shows the trajectories of the target points used in the control law implementation.
A video of the simulations and the experiments can be found on http://anemos.ece.upatras.gr/images/stories/videos/ PTGS IEEETAC.mp4

VII. CONCLUSION
This article examines the area coverage problem by a homogeneous team of mobile agents with imprecise localization.  A gradient ascent based control law was designed based on a Guaranteed Voronoi partition of the region. Because of the complexity of the resulting control law, a simpler suboptimal control law is proposed and the performance of both is compared in simulation studies. Additionally, two experiments were conducted to highlight the efficiency of the suboptimal control law.

APPENDIX
We assume n disks D i , i ∈ I n with centers q i = [x i , y i ] T ∈ R 2 and radii r u i , i ∈ I n .

A. Hyperbola Definition and Properties
We will use the GV diagram definition given in (5) for two disks D i and D j and shown that the boundaries of their GV cells are hyperbola branches.
Let us find the boundary of the cell of node i, V g i = H i j , since we examine only two nodes. We have that ∂ H i j = q ∈ Ω : max q − b i = min q − b j , ∀b i ∈ D i , ∀b j ∈ D j , however since D i and D j are disks, the maximum and minimum distances can be calculated and thus ∂ H i j = q ∈ Ω : q − q i + r u i = q − q j − r u i , where q i and q j are the centers of D i and D j respectively. The equation ∂ H i j = q ∈ Ω : q − q j − q − q i = r i + r j , defines one branch of a hyperbola since it is the locus of points whose difference of distances from two foci q i and q j is equal to a positive constant r i +r j . The other branch of the hyperbola is ∂ H ji = q ∈ Ω : q − q i − q − q j = r i + r j , and corresponds to the boundary of V g j . In a hyperbola, the positive constant is called the major axis and is usually denoted as 2a whereas the distance between its foci is denoted as 2c. A third parameter, the minor axis denoted as 2b is then defined by the equation a 2 + b 2 = c 2 . As such in this case we have that 2a = r i + r j and 2c = q i − q j . Additionally, the closest points on the two branches are called the vertices of the hyperbola and their distance is 2a.

B. Hyperbola Parametric Equation
Assuming the foci of the hyperbola are q i = [−c 0] T and q j = [c 0] T , its parametric equation is where a, b and c are those defined in Appendix A and the +(−) sign corresponds to the West (East) branch of the hyperbola which is H ji (H i j ). By rotating the hyperbola foci around the origin by an angle θ , the parametric equation becomes γ(t) = cos(θ ) − sin(θ ) sin(θ ) cos(θ ) ±a cosh(t) b sinh(t) , t ∈ R Finally, given two foci q i = [x i y i ] T and q j = [x j y j ] T anywhere on the plane, the general hyperbola parametric equation is γ(t) = cos(θ ) − sin(θ ) sin(θ ) cos(θ ) ±a cosh(t) b sinh(t) + x i +x j 2 y i +y j 2 , t ∈ R (20) where θ = arctan( y j −y i x j −x i ).

C. Outward Unit Normal Vectors
For any curve γ(t) the normal vector is defined aŝ The magnitude of this vector is the curvature at that particular point on the curve and it points towards the center of curvature.
Since the hyperbolic branches define convex regions on the plane, the outward unit normal vector is given by However, since the resulting expressions are too large to show here, an indicative plot of the outward unit normal vectors of one branch of a hyperbola is shown in Figure 27.

D. Jacobian Matrix
The Jacobian matrix υ i j T shows the change in the curve γ j (t) = γ jx (t) γ jy (t) caused by the movement of node i, thus its transpose is The expressions for the Jacobian matrix elements are too large to show here and as such some indicative plots of them with respect to the parameter t are shown in Figure 28.