Mathematical Modelling of Turning Delays in Swarm Robotics

1.9 MB·

by  Jake P. Taylor-King

SubmittedtoIMAJournalofAppliedMathematics.
Mathematical Modelling ofTurning Delays inSwarm Robotics
4
1 JAKEP. TAYLOR-KING, BENJAMIN FRANZ, CHRISTIAN A.YATESAND RADEK ERBAN∗
0
2 MathematicalInstitute,UniversityofOxford,

t AndrewWilesBuilding,RadcliffeObservatoryQuarter
c
WoodstockRoad,Oxford,OX26GG,UnitedKingdom
O

1 Weinvestigatetheeffectofturningdelaysonthebehaviourofgroupsofdifferentialwheeledrobotsand
showthatthegroup-levelbehaviourcanbedescribedbyatransportequationwithasuitablyincorporated

] delay.Theresultsofourmathematicalanalysisaresupportedbynumericalsimulationsandexperiments
O
withE-Puckrobots. Theexperimentalquantitywecomparetoourrevisedmodelisthemeantimefor
R robots to find the target area in an unknown environment. The transport equation with delay better
. predictsthemeantimetofindthetargetthanthestandardtransportequationwithoutdelay.
s
c
Keywords:Velocityjumpprocess,Swarmrobotics,Transportequationwithdelay
[

2
1 Introduction
v
2
Much theory has been developed for the coordination and control of distributed autonomous agents,
1
where collections of robots are acting in environments in which only short-range communication is
6
7 possible (ReifandWang, 1999). By performingactions based on the presence or absence of signals,
. algorithmshavebeencreatedto achievesome greatergroupleveltask; forinstance, to reconnoitrean
1
0 areaofinterestwhilstcollectingdataormaintainingformations(Desaietal., 2001). Inthispaper,we
4 willinvestigateanimplementationofsearchingalgorithms,similartothoseusedbyflagellatedbacteria,
1 inaroboticsystem.
:
v ManyflagellatedbacteriasuchasEscherichiacoli(E.coli)usearun-and-tumblesearchingstrategy
i in which movement consists of more-or-less straight runs interrupted by brief tumbles (Berg, 1983).
X
Whentheirmotorsrotatecounter-clockwisetheflagellaformabundlethatpropelsthecellforwardwith
r
a roughly constant speed; when one or more flagellar motors rotate clockwise the bundle flies apart
a
andthecell‘tumbles’(Kimetal.,2003). Tumblesreorientthecellinamore-or-lessuniformly-random
direction (with a slight bias in the direction of the previous run) for the next run (BergandBrown,
1972). In the absenceof signalgradientsthe randomwalk is unbiased, with a meanrun time ∼1sec
andatumbletime∼0.1sec. However,whenexposedtoanexternalsignalgradient,thecellresponds
byincreasing(decreasing)therunlengthwhenmovingtowards(awayfrom)afavourabledirection,and
thereforetherandomwalkisbiasedwithadriftinthatdirection(Berg,1975;Koshland,1980). Similar
behaviourcanbeobservedinswarmsofanimalsavoidingpredatorsandcoordinatingthemselveswithin
agroup(Couzinetal.,2002).
ThebehaviourofE.coliisoftenmodelledasavelocityjumpprocesswherethetimespenttumbling
isneglectedasitismuchsmallerthanthetimespentrunning(Othmeretal.,1988;ErbanandOthmer,
2004). In such a velocity jump process, particles follow a given velocityu from a set of allowed ve-
locitiesV ⊂Rd, d =2,3, for a finite time. The particle changes velocity probabilistically according
to a Poisson processwith intensity l , i.e. the mean run-durationis 1/l . A new velocity v is chosen
∗e-mail:[email protected]
(cid:13)c InstituteofMathematicsanditsApplications2013;allrightsreserved.
2of20 TAYLOR-KING,FRANZ,YATESANDERBAN
accordingtotheturningkernelT(v,u):V×V →R. Formallytheturningkernelrepresentstheproba-
bilityofchoosingvasthenewvelocitygiventhattheoldvelocitywasu. Therefore,itisnecessarythat
T(v,u)dv=1andT >0.
V
R Denotingby p(t,x,v)thedensityofbacteriawhichare,attimet,atpositionxwithvelocityv,the
velocityjumpprocesscanbedescribedbythetransportequation(Othmeretal.,1988)
¶ p
(t,x,v)+v·(cid:209) p(t,x,v)=−l p(t,x,v)+l T(v,u)p(t,x,u)du. (1.1)
¶ t x ZV
Assumingthatl andT areconstant,onecanshowthatthelong-timebehaviourofthedensity̺(t,x)=
p(t,x,v)dv is given by the diffusion equation (HillenandOthmer, 2000). If l depends on an ex-
V
Rternal signal (e.g. nutrient concentration), then the resulting velocity jump process is biased and its
long time behaviour can be described by a drift-diffusion equation for ̺ (OthmerandHillen, 2002;
ErbanandOthmer,2005).
Inthispaper,wewillstudyanexperimentalsystembasedonE-Puckrobots(BonaniandMondada,
2004). Weprogrammethesedifferentialwheeledrobotstofollowarun-and-tumblesearchingstrategy
inordertofindagiventargetset. Inthefirstsetofexperiments,weconcentrateonthesimplestpossible
scenario:anunbiasedvelocityjumpprocessintwospatialdimensionswiththefixedspeeds∈R+,the
constantmeanruntimel −1∈R+,andtheturningkernelwhichisindependentofu
d (||v||−s)
T(v,u)= . (1.2)
2p s
A special feature of the E-Puck robots is that they can perform turns on the spot as in the classical
velocity jump process described by (1.1). In this paper, we will investigate in how far (1.1) presents
a good description of the behaviour of the robotic system and we will develop an extension of (1.1)
thatresultsin a bettermatchbetweenexperimentaldata andmathematicalmodel. We thenapplythis
extended velocity jump theory to a biased random walk through the incorporationof signals into the
experimentalsetup.
The paper is organizedas follows: in Section 2, we introduce the experimentalsystem as well as
theobtaineddata. Thisdataiscomparedtotheclassicalvelocityjumptheory. InSection3,weextend
the velocityjump theory to includefinite turningtimes for unbiasedrandomwalks and compareit to
ourexperimentaldata, showinga muchimprovedmatch. ThisnewtheoryisinSection4appliedtoa
situationwithanexternalsignalandthereforeabiasedrandomwalk.Weconcludeourpaper,inSection
5,bydiscussingtheimplicationsofourresults.
2 Velocityjumpprocessesinexperimentswithrobots
Equation(1.1)introducedthedensitybehaviourofthegeneralvelocityjumpprocessthatweareaiming
to investigate using the experimental set-up described in Section 2.1. In particular, we will initially
concentrate on a simple unbiased velocity jump process with the fixed speed s∈R+, the mean run
durationl −1∈R+andtheturningkernel(1.2). InSection4wewillpresentsituations,wheretheturn-
ingfrequencychangesaccordingtoanexternalsignal,asisindeedcommoninbiologicalapplications
(ErbanandOthmer,2005). Thisfixed-speedvelocityjumpprocesscanbeviewedasastartingpointfor
consideringmorecomplexsearchingalgorithms.Wewilldemonstratethatbyincludingasmallmodifi-
cation(theintroductionofadelaytotheturningkernel),wecanalterthissimplevelocityjumpprocess
sothatitmodelsthebehaviouroftheE-Puckrobots.
MathematicalModellingofTurningDelaysinSwarmRobotics 3of20
Weareinterestedincomparingtheidealisedvelocityjumpprocess,givenin(1.1)–(1.2),torobotic
experiments. Due to a restriction in numbersof robots, one cannotfeasibly talk abouta “density” of
robots that could be compared to p(t,x,v) as given in (1.1). Therefore, our experiments concentrate
ontheescapeofrobotsfromagivendomain. Wemayinterpretthisasthetargetfindingabilityofthe
E-Puck robots. Usingthese experiments,we caninferdatabothonthe fluxatthebarrierandtheexit
timesandcancomparethosetonumericalresultsofvelocityjumpprocessesinSections2.3and2.4.
2.1 Experimentalset-upandprocedure
Toobtaintheempiricaldata,anexperimentalsystemconsistingof16E-Puckrobotswasused. E-Puck
robots are small differential wheeled robots with a programmable microchip (BonaniandMondada,
2004).Thediameterofeachrobotise =75mmwithaheightof50mmandweightof200g. Throughout
the experiments, the speed was chosen to be s=5.8×10−2m/sec. The robots turn with an angular
velocityw =4.65/sec. FullspecificationsalongwithapicturearegiveninAppendixA.
Intheexperiments,weusearectangulararenaW withwallsonthreeofthe4edgesandanopening
to the targetarea T alongthe fourthedge1. A diagramof thearenaalongwiththe notationusedcan
be seen in Figure 1 and a photo is shown in Figure 5(b) in AppendixA. When considering such an
arena, onehas to distinguishbetween the size ofthe physicalarenaand the effectivearena (shownin
blueinFigure1)thattherobotcentrescanoccupy. Theeffectivearenausedintheexperimentshasthe
dimensionsL =1.183mandL =1.145m=L −e /2. Thereflective(wall)boundaryandthetarget
x y x
boundarywillbedenotedas¶W R and¶W T,respectively,andcanbedefinedas
¶W T =W ∩T , ¶W R =¶W \¶W T . (2.1)
Throughouttheremainderofthepaper,wewillalsousenR (resp.nT)todenotetheoutwardspointing
normalonthereflective(resp.target)boundary.
Duringtheexperiments,robotswereinitialisedinsidearemovablesquarepenW ofeffectiveedge
0
lengthL =0.305m,showninFigure1andFigure5(b)inAppendixA.Ashortperiodoffreemovement
0
withinthepenbeforeitsremovalallowedustoreliablyreleaseallrobotsintothefulldomainW atthe
same time as well as randomising their initial positions within the pen. We recorded the exit time
for each of the robots, when its geometric centre entered the target area T. Each repetition of the
experimentwascontinuedfor300secoruntilall16robotshadleftthearena.
The robots were programmedusing C and a cross-compilingtool, with the firmware being trans-
ferred onto the robots via bluetooth. A pseudo-code of the algorithm implemented on the robots is
showninTable1. ThisalgorithmrepresentsavelocityjumpprocessinthelimitasD t→0(Erbanetal.,
2006),andgivesagoodapproximationaslongaslD t≪1.Intheexperimentsweusedl =0.25sec−1
implyingameanrundurationof4secandD t=0.1sec,resultinginlD t=2.5×10−2≪1. Notethat,
s, w andl canbechangedonasoftwarelevelonanE-Puck. Forw wechosethemaximumpossible
value, whilst for s we chose a value below the physicalmaximum. Choosing a lower velocitymeans
thatwemitigatetheeffectsofaccelerationanddecelerationtotherunningspeedsincetherobotscan-
notdothisinstantaneouslyasthebasicvelocityjumpmodelassumes. Inapracticalsetting,onecould
interpretsandw asgivencharacteristicsofthesystem,whilstl canbechoseninawaythataccelerates
the targetfindingprocessforthe givenapplicationwith thechoiceof l likely torepresenta trade-off
betweensamplinganareaandtimespentreorienting.
1WehaveW ∩T =0/,butW ∩T 6=0/,i.e.W andT touchbutdonotoverlap.
4of20 TAYLOR-KING,FRANZ,YATESANDERBAN
L =1.183m
x
¶W
R
e =7.5cm
T
W ¶W
Ly=1.145m 0 0.380m T 1.220m
L0=0.305m W
FIG.1. Schematicshowingoftheexperimentalset-upalongwiththenotationusedthroughoutthispaper. Dottedborderlines
correspondtotheeffectivearenaandboldlinestotheactualarena.Forfurtherdetailsseethetext.
InadditiontothealgorithminTable1,robotswerealsomadetoimplementanobstacle-avoidance
strategyusingthe fourproximitysensorsplacedat angles±17.5◦ and ±47.5◦ fromthecentre axisin
the front part of the E-Puck. Reflective turns were carried out based on the signals received at these
sensors. Astherobotsareincapableofdistinguishingbetweenwallsandotherrobots,thosereflections
occurwhetherarobotcollideswiththewall¶W R oranotherrobot. Asa consequencewe discussthe
importanceofrobot-robotcollisionsontheexperimentalresultsinthenextsection.
2.2 Relevanceofcollisionsforlownumbersofrobots
Fornon-interactingparticleswhichcanchangedirectioninstantaneously,equation(1.1)accuratelyde-
scribesthemesoscopicdensitythroughtime. However,inourexperimentstherobotsundergoreflective
collisionswhentheycomeintoclosecontact,ratherthanpassingthroughorovereachother. Foralow
numberofparticles,weusedMonteCarlosimulationstodemonstratethatcollisionsarenotthedomi-
nantbehaviourandhavelittleeffectonthedistributionofparticles.Inpanels(a)and(b)ofFigure2,we
comparetwoMonteCarlosimulations: (a)inwhichparticlesareallowedtopassthroughoneanother
and(b)inwhichcollisionsaremodelledexplicitly. InFigure2(c)wepresentthesolutionofequation
(1.1). Thiscomparisondemonstratesthatthemeandensityoftheunderlyingprocessconvergestothe
solutionoftransportequation(1.1). Theparametersemployedinthismodelcomparisonaretakendi-
rectlyfromtheequivalentrobotexperiment;(s,l ,e )=(5.8×10−2m/sec, 0.250sec−1, 7.5×10−2m).
In Figure 2(c), for the differentialequation, we use a first-order numericalscheme with Dq =p /20,
MathematicalModellingofTurningDelaysinSwarmRobotics 5of20
[1]Robotisstartedatpositionx(0)∈W . Generater ∈[0,1]uniformlyatrandom,sett=0and
0 1
cos2p r
v(0)=s 1 .
(cid:18) sin2p r1 (cid:19)
[2]Positionisupdatedaccordingtox(t+D t)=x(t)+D tv(t).
[3] Generate r ∈[0,1] uniformly at random. If r <lD t, then generate r ∈[0,1] uniformly at
2 2 3
randomandset
cos2p r
v(t+D t)=s 3 .
(cid:18) sin2p r3 (cid:19)
[4]Sett=t+D t andcontinuewithstep[2].
Table1.Analgorithmicimplementationofthevelocityjumpprocess.
(a) (b) (c)
Ly 1.5 Ly 1.5 Ly 1.5
1 1 1
Ly/2 Ly/2 Ly/2
0.5 0.5 0.5
0 0 0
0 0 0
0 Lx/2 Lx 0 Lx/2 Lx 0 Lx/2 Lx
FIG.2. Comparisonofindividual-based simulations with(1.1). Eachplotshowstheresultingdensityatthefinaltimeofthe
simulation,20sec. (a)Individual-basedsimulationusing16×4×104pointparticles. (b)Individual-basedsimulation,average
over4×104runsusing16particleswithhard-sphereinteractions. (c)Numericalsolutionto(1.1)usingafinitevolumemethod
withparametersgiveninthetext.
D x=L /200andD t=10−2sec.
y
In the Monte Carlo simulations we initialise particles in the effective pen for 20sec where they
undergohard-spherecollisions. They are then released into the larger arena where in one simulation
they are point-particlesand in the other they undergoreflective collisions as hard-spheres. Instead of
removing particles at the target boundary as shown in Figure 1 (as we do in the robot experiments),
this edge of the domain is closed so that all edges correspondto reflective boundaryconditions. For
transportequation(1.1),wemodeltheinitialconditionasastepfunctionoverthepen. Thesedensities
arevisualisedinFigure2. Formally,thisinitialconditioncanbewrittenas
c W d (||v||−s)
p(0,x,v)= 0 , (2.2)
L22p s
0
wherec W 0denotestheindicatorfunctionoftheinitialregionW 0.Thecorrespondingboundarycondition
is p(t,x,v)=p(t,x,v′)forx∈¶W R wherethereflectedvelocityv′isdefinedas
v′=v−2(v·nW )nW , (2.3)
wherenW istheoutwardpointingnormalatthepositionx∈¶W .
6of20 TAYLOR-KING,FRANZ,YATESANDERBAN
After 20sec, we record the density in each of the scenarios and present the results in Figure 2.
There is minimal visible discrepancybetween the Monte Carlo simulations presented in Figure 2 for
our choice of parameter values. In order to compare the three simulations given in Figure 2 we also
employedapairwiseKolmogorov-Smirnovtest(Peacock,1983). Avalue(oftheKolmogorov-Smirnov
metric)closetozerodenotesagoodfitbetweenthetwosimulations. Itcorrespondstotheprobability
thatonecanrejectthehypothesisthatthedistributionsareidentical. Whencomparingthetwo Monte
Carlosimulations,avalueof2.37×10−2wasobtained;whencomparingequation(1.1)withthehard-
sphereMonteCarlosimulation,avalueof5.65×10−2wasobtained;finallywhencomparingequation
(1.1)withthepoint-particleMonteCarlosimulation,avalueof3.40×10−2wasobtained.Thissupports
thevisualobservationthatallthreedensitydistributionsareallhighlysimilar.
InthelimitwhereN→¥ ,forN beingthenumberofrobots,transportequation(1.1)canbealtered
bytheadditionofa Boltzmann-likecollisionterm(Harrisi,1971;Cercignani, 1988). Itcanbeshown
thattheeffectsofcollisionsbetweenrobotsarenegligibleforthepresentedstudy(Franzetal.,2014).
2.3 Comparisonbetweentheoryandexperiments: lossofmassovertime
In thisand subsequentsectionswe comparethe resultsof 50repetitionsof the experimentsdescribed
in Section 2.1 with numerical results obtained by solving the correspondingmathematical equations.
Onewayofinterpretingtheexperimentalexit-timedataisbyconsideringtheexpectedmassremaining
insidethearenaW atagiventime. Fortheexperimentaldatathisquantityisplottedasasolid(black)
line in Figure 3(a). We compare this result to the variation of the remaining mass with time from a
numericalsolutionof(1.1)combinedwiththefollowingboundaryconditions:
p(t,x,v)=0, x∈¶W T,v·nT <0,
(2.4)
p(t,x,v)=p(t,x,v′), x∈¶W R,
wherethereflectedvelocityv′isdefinedby(2.3).AsdemonstratedinSection2.2,suchacomparisonis
reasonablesincecollisionsdonothaveamajorimpactintheparameterregimechosenhere. Theinitial
condition for transport equation (1.1) is identical to the condition given in equation (2.2). The mass
remaininginthedomainisthendefinedas
m(t)= p(t,x,v)dxdv,
ZW ZV
and is plotted as a dotted (red) line in Figure 3(a). The initial mass is normalized to 1. An obvious
observationfromFigure3(a)isthatthetransportequationdescriptiondoesnotmatchtheexperimental
data well, with the robotsexiting the arena significantly slower than predicted. In this figure, we use
a first-orderfinitevolumemethodwithDq =p /20,D x=1.183m/200andD t =10−3sec inorderto
solvetransportequation(1.1).
2.4 Comparisonbetweentheoryandexperiments: meanexittimeproblem
An alternative way to interpretthe experimentaldata is to consider mean exit times. Throughoutthe
experimentsonly708of the 800(=50×16)robotsleft the arena beforet =300sec. The average
end
exit time of those 708 robotswas 121.92sec. In order to be able to compareexperimentalexit times
withthemeanexittimeproblems,itisnecessarytoestimatethemeanexittimeofall800robots.Using
thebestexponentialfitonthemassovertimerelation(cf. Figure3(a)),wecanestimatethemeanexit
MathematicalModellingofTurningDelaysinSwarmRobotics 7of20
(a) (b)
150
1
0.8
100
c
s0.6 e
s s
Ma n
i
0.4 t
50
0.2
0 0
0 50 100 150 200 250 0 0.2 0.4 0.6 0.8 1
t in sec x in m
FIG.3. (Colouravailable online.) Comparingtheexperimental results(solidblackline)tosolutions of(1.1)and(3.2)–(3.3).
Panel(a)showstherelativemassofthesystemovertime. Thedottedline(red)showsthenumericalsolutiontoequation(1.1)
withboundaryconditions(2.4),thedashedline(blue)showsthenumericalsolutionofthesystemofequations(3.2)–(3.3)with
boundaryconditions(3.4).Panel(b)showsthemeanexittimeaveragedoverallvelocitydirectionsvsthex-coordinatealongthe
arenaedge. Theadsorbingboundaryisatx=1.183m. Thedottedline(red)showsthemeanexittimecomputedusingequation
(2.5) with boundary conditions (2.6), the dashed line (blue) shows the mean exit time computed using equation (3.21) with
boundaryconditions(3.22). Inordertoallowdirectcomparisonwiththeexperimentaldata,theshorterboldlinesrepresentthe
averageoftheoreticallyderivedexittimesovertheregionW 0,fromwhichtherobotswerereleasedintheexperimentalscenario.
Forbothplotsparametersandnumericalmethodsaredescribedinthetext.
time of the remaining92 robotsto be 424.69sec. The approximatemean exittime established in the
experimentsistherefore156.74sec;thisvalueisplottedasthesolid(black)lineinFigure3(b). Inorder
tobeabletocomparethisvaluetoanalyticresults,onehastoreformulatethetransportequation(1.1)
intoameanexittimeproblem. Letusthereforedefinethemeanexittimet =t (x ,v )ofarobotthat
0 0
startsatpositionx ∈W withvelocityv ∈V. Thismeanexittimesatisfiesthefollowingequation
0 0
v ·(cid:209) t (x ,v )−lt (x ,v )+l T(u ,v )t (x ,u )du =−1. (2.5)
0 x0 0 0 0 0 ZV 0 0 0 0 0
InSection3,inwhichdelaysaremodelled,aderivationisgivenforthemeanexittimeproblem;setting
the delay term to zero allows one to see how equation (2.5) is derived. This so-called “backwards
problem”satisfiesthefollowingboundaryconditions
t (x0,v0)=0, x0∈¶W T ,v0·nT >0,
(2.6)
t (x0,v0)=t (x0,v′0), x0∈¶W R,
wherev′ isagainthereflectedvelocitywithrespecttov asdefinedin(2.3). Duetothearenashape,by
0 0
takingthespatialaverageinthey-direction
1 Ly/2
t (x ,v )= t (x ,y ,v )dy , (2.7)
x 0 0 0 0 0 0
LyZ−Ly/2
onecanfurthersimplifythemeanexittime problem. Inthecasewheretheturningkernelisgivenby
equation(1.2),onecanobtainaproblemwithtwoparametersx andq ,whereq ∈(−p ,p ]istheangle
0
8of20 TAYLOR-KING,FRANZ,YATESANDERBAN
definingthevelocityv byv =s(cos(q ),sin(q )).Fort =t (x ,q )
0 0 x x 0
¶t l p
scos(q ) x −lt + t (x ,f )df =−1,
¶ x0 x 2p Z−p x 0
t (0,q )=t (0,p −q ), (2.8)
x x
p p
t (L ,q )=0, q ∈ − , .
x x
2 2
h i
Wheninitialdirectioncannotbespecified,themean-exittimefromagivenx-positionisgivenby
1 p
t (x ,q )dq ,
2p Z−p x 0
wheret isthesolutionof(2.8). Thisisplottedasthedotted(red)lineinFigure3(b). Thenumerical
x
solution was performed using an upwind-scheme in the x-direction with D x=1.1825m/200 and an
angulardiscretisation of Dq =p /20. Additionally, we take the spatial averageof the mean-exittime
fromtheinitialregionW andplotthisasthebolddashedlineinFigure3(b). Thislinedoesnotmatch
0
well with the corresponding average mean-exit time found in the robot experiments. The numerical
solution of equation (2.8) predicted a mean exit time of 137.49sec, meaning an underestimation of
19.25secor12.3%comparedto theexperimentalexittimeof156.74sec. Inthefollowingsectionwe
willextendtheclassicalvelocityjumptheorytoimprovethismatchwiththeexperimentaldata.
3 Modellingturningdelays
InSection2.2,weobservedthatcollisionsbetweenrobotsdoesnotplayamajorroleinexplainingthe
discrepancybetweenthetransportequation(1.1)andtheexperimentaldatapresentedinSections2.3and
2.4. Aswellasassumingindependentlymovingparticles,thetransportequation(1.1)isalsopredicated
ontheassumptionthatthereorientationphasetakesanegligibleamountoftimecomparedtotherunning
phase. Sincethisassumptionisnotsatisfiedinourrobotexperiments,thissectionextendstheoriginal
modelthroughtheinclusionoffiniteturningtimes.
3.1 Introductionofarestingstate
Letusinitiallystatetwoassumptionsthatapplytotherobotexperiment,butmightnotextendtovelocity
jump processesin biologicalsystems, like the run-and-tumblemotion of E. Coli (Berg, 1983), which
hasmotivatedthesearchingstrategiesimplementedonrobots:
(a) a new direction v′ ∈V is chosen as soon as the particle enters the reorientation (“tumble”)
phase;
(b)thetimeittakesforaparticletoreorient(“tumble”)fromvelocityv∈V tov′∈V isspecified
bythefunctionK(v′,v):V×V 7→R+.
Assumption (b) implies that the turning time is constant in time and equal for each particle and, in
particular, does not depend on the particle’s history. For the robots studied in this paper, we can ad-
ditionally assume that reorientation phase is equivalent to a directed rotation with a constant angular
velocityw ∈R+. Therefore,theturningtime dependsonlyonthe anglebetweenthe currentvelocity
v∈V andthenewvelocityv′∈V andKtakestheform
1 v·v′
K(v′,v)= arccos . (3.1)
w (cid:18)||v|| ||v′||(cid:19)
MathematicalModellingofTurningDelaysinSwarmRobotics 9of20
Wenowextendtheclassicalmodel(1.1)throughtheintroductionofarestingstater(t,x,v,h )thatfor-
mallydefinesthenumberofparticlescurrently“tumbling”(turning)towardstheirnewchosenvelocity
vandremainingturningtimeh . Thedensity p(t,x,v)willnowonlydenotetheparticleswhichareat
timet intherunphase. Theupdateoftheextendedsystemisgiventhrough
¶ p
(t,x,v)+v·(cid:209) p(t,x,v)=−l p(t,x,v)+r(t,x,v,0+), (3.2)
¶ t x
¶ r ¶ r
(t,x,v,h )− (t,x,v,h )=l p(t,x,u)T(v,u)d (h −K(v,u))du. (3.3)
¶ t ¶h ZV
In(3.2)wecanseethatrunningparticleswillenteratumblephasewithratel andparticlesthathave
finished the tumblesignified throughh =0 will re-enterthe run-phase. Equation(3.3) representsthe
linear relation between t and h and shows that particles enter the tumble phase depending on their
newlychosenvelocitydirection. Inordertoguaranteeconservationofmassthroughoutthesystem,we
introducethenon-negativityconditionforh through
r(t,x,v,h )=0, for t>0, x∈W , v∈V and h <0.
Additionally,theboundaryconditionsforthesystem(3.2)–(3.3)aregiventhrough
p(t,x,v)=0, x∈¶W T , v·nT <0,
p(t,x,v)=−r(t,x,v,0+)/(v·nR), x∈¶W R, v·nR<0,
¶ r ¶ r (3.4)
¶ t(t,x,v,h )−¶h (t,x,v,h )=d (h −K(v,v′))(v′·nR)p(t,x,v′), x∈¶W R, v·nR<0,
r(t,x,v,h )=0, x∈¶W R, v·nR>0,
where v′ is the reflected velocity of v given by (2.3). In order to show that the system (3.2)–(3.3) is
actuallyconsistent,weprovethatmassinthesystemisconservedifnotargetispresent.
LEMMA3.1 Thetotalmassinsystem(3.2)–(3.3)withtheboundaryconditionsgivenin(3.4)inthecase
ofreflectiveboundarieseverywhere(¶W R =¶W ,¶W T =0/)giventhrough
¥
M(t)= p(t,x,v)dvdx+ r(t,x,v,h )dh dvdx,
ZW ZV ZW ZVZ0
isconserved.
Proof. Wedefineforeverypointx∈¶W R thetwosubsetsV+andV− ofV asfollows
V+(x)={v∈V : v·nR>0}, V−(x)={v∈V : v·nR<0}. (3.5)
Additionally,letusdefine
¥
R(t,x,v)= r(t,x,v,h )dh .
Z0
Integrating(3.3)withrespecttoh ∈[0,¥ ),weobtainafterreorderingforx∈/¶W :
¶ R
(t,x,v)=−r(t,x,v,0+)+l p(t,x,u)T(v,u)du.
¶ t ZV
10of20 TAYLOR-KING,FRANZ,YATESANDERBAN
Hence,foreverypointx∈/¶W weobtain
¶
[p(t,x,v)+R(t,x,v)]=−l p(t,x,v)+l p(t,x,u)T(v,u)du−v·(cid:209) p(t,x,v).
¶ t ZV x
Integratingthiswithrespecttox∈W andv∈V gives
¶
[p(t,x,v)+R(t,x,v)]dvdx=− v·(cid:209) p(t,x,v)dvdx. (3.6)
ZW ZV ¶ t ZW ZV x
Usingthedivergencetheorem,wecanevaluatetheintegralontherighthandsidetobe
v·(cid:209) xpdvdx = (v·nR(x))p(t,x,v)dvdx
ZW ZV Z¶W ZV
= (v·nR)p(t,x,v)dvdx+ (v·nR)p(t,x,v)dvdx
Z¶W ZV+(x) Z¶W ZV−(x)
= (v·nR)p(t,x,v)dvdx− r(t,x,v,0+)dvdx,
Z¶W ZV+(x) Z¶W ZV−(x)
wherewe haveusedthesecondboundaryconditionin (3.4)inthelaststep. Additionally,forx∈¶W
andv∈V−(x),weobtainbyintegratingthethirdboundaryconditionin(3.4)withrespecttoh ∈[0,¥ )
¶ R
¶ t (t,x,v)=−r(t,x,v,0+)+(v′·nR)p(t,x,v′).
Integrating this with respect to x∈¶W and v∈V and using the last boundary condition in (3.4) we
obtain
¶ R ¶ R
dvdx = dvdx
Z¶W ZV ¶ t Z¶W ZV−(x) ¶ t
= − r(t,x,v,0+)dvdx+ (v′·nR)p(t,x,v′)dvdx
Z¶W ZV−(x) Z¶W ZV−(x)
= − r(t,x,v,0+)dvdx+ (v·nR)p(t,x,v)dvdx. (3.7)
Z¶W ZV−(x) Z¶W ZV+(x)
Summinguptheresultsfrom(3.6)and(3.7),weobtaindM/dt=0andhencethetotalmassM(t)inthe
systemisconserved. (cid:3)
3.2 Transportequationwithturningdelays
We eliminate the resting state from system (3.2)–(3.3) and derive the generalization of the transport
equation(1.1)toa transportequationwithasuitablyincorporateddelay. Thiscanbedonebysolving
(3.3)forrusingthemethodofcharacteristics,whichresultsin
r(t,x,v,0)=r(0,x,v,t)+l T(v,u)p(t−K(v,u),x,u)H(t−K(v,u))du, (3.8)
ZV
whereH istheHeavisidestep function. LetusassumethatK(v,u)isgivenby(3.1). ThenK(v,u)6
p /w . Consideringtimest >p /w , we have r(0,x,v,t)=0. We cannow substitute (3.8) into (3.2) to
obtain
¶ p
(t,x,v)+v·(cid:209) p(t,x,v)=−l p(t,x,v)+l T(v,u)p(t−K(v,u),x,u)du, (3.9)
¶ t x ZV