Chapter 8 The DEVS Integrator: Motion in Space

1. Chapter 8 The DEVS Integrator: Motion in Space • From pulse-based to quantum-based integrators • Mapping Ordinary Differential Equations into DEVS • one-D DEVS Integrators • two-D DEVS integrators • Second Order Systems • Vectors • Motion in Space • Examples – Pursuer Evader, Bird Flocking

2. out out out in out out in rate rate rate varGen varGen Ratefn1 Ratefn2 DEVS Integrator – alternative representation of continuous systems with discrete events varGen sum pulses Instead of summing up pulses, the integrator Receives the fate and adds quanta on its own – with internal transition function When the rate changes, the Ratefn tells the integrator the new rate and it updates its sum based on the old rate and elapsed time. Ratefn Integrator Ratefn rate Converting from sum to integrate Ratefn = Ratefn1 +Ratefn2 Combining rateFns for the same sum into one equivalent

3. input d s(t) / dt = x(t) x y(t) = s(t). x ( t ) ò contents = input d t y y ( t ) x(t ) t s d s /dt s ò 1 1 f x x 1 s d s /dt s ... ò 2 2 f x 2 s x d s /dt s ò n n f n Continuous Integration concepts and Differential Equations d s1(t)/dt = f1(s1(t), s2(t), ..., sn(t), x1(t), x2(t),..., xm(t)) d s2(t)/dt = f2(s1(t), s2(t), ..., sn(t), x1(t), x2(t),..., xm(t)) ... d sn(t)/dt = fn(s1(t), s2(t), ..., sn(t), x1(t), x2(t),..., xm(t))

4. in x newInp x out =s out =s+qleft ta=|q/inp| qleft =q 0 s =s+qleft ta=|q/inp| qleft =q in x newInp x s s+e*inp qleft =qleft - |e*inp| ta = |qleft/newInp| Mapping Differential Equation Models into DEVS Integrator Models DEVS Integrator DEVS instantaneous function

5. output= s+qleft; input>0 q Generate output ta(s) = |quantum/input| Time advance = |qLeft/lastInput)| input<0 s’ s Make a transition state s = s + qLeft qleft = quantum; One Dimensional Integrator:Internal Transition /Output Generation -- integrator reaches predicted threshold and generates this threshold as output; the integrator predicts its next boundary crossing. output output threshold level q time Onde Dim Integrator uses input (derivative) to predict next boundary crossing

6. input Make a transition s = s + e*lastInput; qleft =qleft - |e*lastiinp| lastInput = input; ta = |qleft/input| elapsed time, e Time advance, ta One Dimensional Integrator:Response to External Input – an integrator receives a new input, updates its state, predicts a next state and reschedules itself accordingly to reach the next threshold Residual: QuantumLeft Method:-- limits time to next output by reducing quantumLeft with each external input. quantization. oneDFrom2DInt inp>0 s s+e*lastInp qleft =qleft - |e*lastInp| ta = |qleft/nput| q s inp<0 e First Order Integrator uses last input (derivative) to update to current state Residual: QuantumLeft Method: sigma = |qleft/inp| qleft =qleft - |e*inp| qleft =qleft – e*qleft/sigma =qleft (1– e*sigma) time of current external event time of last event

7. Second order system instantaneous function (identity) in in out out ò ò dv/dt = a ds/dt = v input (a = acceleration) output (v = velocity) output (s =position) restoring force (k = spring constant) friction (air resistance, damping) f(a,v,s) a v s ò ò dv/dt = a – b*v – k*s = f(a,v,s) ds/dt = v -b*v -a*s

8. Bunjee jumper modeled by a second order system 0 f(a,v,s) a v s ò ò -b*v The spring restoring force Only takes effect when The jumper reaches the stretch length L s dv/dt = a – b*v - F(s) = f(a,v,s) ds/dt = v L

9. Composite 1D integrator: encapsulates both position and velocity in one atomic model; updates both components at once. oneDCompInt in out input (acceleration) sate.s state = (s,v) after elapsed time, e, update both s and v s = s + v*e, v = v+ a*e ta = min {q/a, 10q/v} dv/dt = a – b*v – k*s ds/dt = v

10. Modeling a one dimensional car accelerate: dv/dt = f(a,v) (v,a > 0) ds/dt = v out speed bounded: f(a,v) = a if v < vmax = 0 else v goes to 0 brake: dv/dt = - b*v (b>0) ds/dt = v (dv/dt = 0 implies v => 0) brake accelerate

11. Class Vector public double dotProd(vect2DEnt v){ return x*v.x + y*v.y; } public vect2DEnt perpendicular(){ return new vect2DEnt(-y,x); } public vect2DEnt scalarMult(double d){ return new vect2DEnt(x*d,y*d); } public vect2DEnt maxLimit(double max){ if (norm() <= max) return this; else{ vect2DEnt na = normalize(); return na.scalarMult(max); } } public double norm(){ return Math.sqrt(x*x + y*y); } public vect2DEnt normalize(){ double norm = norm(); if (norm == 0) return ZERO; else return new vect2DEnt(x/norm,y/norm); } public class vect2DEnt extends entity { public double x,y; public static vect2DEnt ZERO = new vect2DEnt(0,0); public vect2DEnt(double x, double y){ super("vect2DEnt"); this.x = x; this.y = y; } public vect2DEnt add(double x, double y){ return new vect2DEnt(this.x+x,this.y+y); } public vect2DEnt add(vect2DEnt v){ return add(v.x,v.y); } public void addSelf(vect2DEnt v){ x += v.x; y += v.y; } public vect2DEnt subtract(vect2DEnt v){ return add(-v.x,-v.y); } quantization vector2DEnt

12. vect2DEnt.ZERO v = new vect2DEnt(1,1) v.add( new vect2DEnt(1,1-)) v’.y v.y Check: v.perpendicular().dotProd(v) = 0 new vect2DEnt(1-,1) v.perpendicular) v.x v’.x v.dotProd(v’) Check: v.normalize() .norm() = 1 v.perpenducar().norm() =v.norm(); v.scalarMult(d).norm() = d*norm() v.scalarMult(d)) v.norm()

13. output= nextstatel Generate output Make Transition; sigma = timeAdvance(quantum); remainingQuant = quantum; nextState = state.add(input.scalarMult(sigma)); s’ s 2 Dimentional Integrator -Internal Transition /Output Generation -- integrator reaches predicted state and generates this state as output; the integrator predicts its next state. timeAdvance = quantum/input.norm() input vector quantum is a radius quantization twoDint

14. input Make a transition state =state.add(lastInput.scalarMult(e)) lastInput = input; quantumLeft -= lastInput.scalarMult(e).norm() ta = quantumLeft/input.norm(); Response to External Input – an integrator receives a new input, updates its state, predicts the next state and reschedules itself to reach the next state. elapsed time, e Time advance, ta Input arrives after e update state by vector adding lastInput*e quantumLeft -= lastInput.scalarMult(e).norm() ta = quantumLeft /input.norm()

15. Pursuer Evader – vehicles controlled by perceiving agents • evader moves • perpendicular to • line joining vehicles • agents send position • updates at their • quantum threshold • crossings • pursuer moves • in direction of • line joining vehicles larger quantum size for an agent gives it the advantage over its opponent

16. Pursuer Evader (cont’d) • pursuer’s acceleration vector is the difference • between its position and the desiredPosition coming • in from evader • evader’s acceleration vector is perpendicular to • its difference vector public class driver extends ViewableAtomic{ protected vect2DEnt position, desiredPosition, difference; protected double tolerance = 2,maxPos = 100; public void deltext(double e ,message x){ Continue(e); for (int i=0; i< x.getLength();i++){ if (messageOnPort(x,"position",i)) { entity en = x.getValOnPort("position",i); position = vect2DEnt.toObject(en); position.addSelf(new vect2DEnt(r.uniform(-10,10),r.uniform(-10,10))); } else if (messageOnPort(x,"desiredPosition",i)) { entity en = x.getValOnPort("desiredPosition",i); desiredPosition = vect2DEnt.toObject(en); } } difference = desiredPosition.subtract(position); if (name.startsWith("evade")) difference = difference.perpendicular(); if (difference.norm()<= tolerance) holdIn("doBrake",0); else holdIn("doAccel",0); } }

17. ? v2 p1 p2 Class VectAim – computes time of intersection, and the velocity required to intersect in a given time, with a moving target p to intersect: v must be in same direction from origin as p let vect = vect2DEnt public boolean vect ::sameDir (vect v): has my direction public static double timeToIntersect (vect p,vect v): if v is in same direction as p then return p.norm()/v.norm(), else infinity public static double timeToIntersect (vect p1,vect p2, vect v1,vect v2){ return timeToIntersect(p2.subtract(p1),v1.subtract(v2)); //note sign diff } public static vect predictLoc (vect p2, vect v2,double ta): //predict the location vehicle will be in time ta {return p2.add(v2.scalarMult(ta));} public static vect dir (vect p1,vect p2p) //the direction from p1 to p2p {return p2p.subtract(p1).normalize();} public static double speed (vect p1,vect p2p,double ta) //the speed required to reach p2p from p1 in time ta {return p2p.subtract(p1).norm()/ta;} public static vect aim (vect p1,vect p2, vect v2,double ta) //the velocity required to intersect with target (p2,v2) in time ta { vect p2p = predictLoc (p2, v2,ta); vect dir = dir(p1,p2p); double speed = speed(p1,p2p,ta); return dir.scalarMult(speed); } v will (v1,p1) collide with (v2,p2)? v1 v2 p1 p2 p2’ = p2-p1 v2’=v1-v2 • place the origin at p1 • p2’ = p2-p1 is location of p2 in this frame • v2’ = v1-v2 is the relative velocity aiming at p2 • for the relative velocity v1’ • to be in the direction of p2’, the perpendicular components must • cancel

18. Missiles– pursuer uses sensing and time-critical prediction to continually re-aim for target • dumb pursuer – flies fixed trajectory • toward target • smart pursuer -- adjusts velocity to intersect with target’s predicted trajectory in a given time. The time left is reduced as time passes without reaching target • target has collision detector which • continually checks for intersection • of trajectories Pursuer uses this logic: if (myPos.subtract (otherPos).norm() < takeEffectDist){ vect vHit = vectAim.aim (myPos, otherPos, otherVel, Math.max (1,timeLeft)); difference = vHit.subtract (myVel). scalarMult (power); holdIn("doAccel",0); } dumb pursuer smart pursuer

19. 3 vehicleSpace vehicles cellGridPlot vehicle driver vehicleDyn • Start with vehicleSpace which is the coupled model for • Pursuer-Evader. • Modify it as described below. • flockSpace contains any number of birds with all-to-all • coupling of “out” to “in”. • birdDriver moves toward a desiredPosition that is a • combination of the positions of the other birds as shown below. • a lead bird provides direction and speed Bird flocking behavior same as vehicleSpace except for multiple birds flockSpace birds cellGridPlot bird same as vehicle except for driver birdDriver vehicleDyn NOTE: all operations use Vect2DEnt define birdDriver by adapting driver as follows: when receive a bird’s position: if (it is farther from me than a threshold){ desiredPosition = 0.5*(desiredPosition + bird’s position) difference = desiredPosition - position If (name.startsWith("bird_0")) difference = difference + vector in direction NE with magnitude = 100 (ALL OTHER ASPECTS OF driver REMAIN THE SAME)

20. Bird flocking behavior Problem: For 6 birds find value of threshold that keeps flock closest together while following lead bird Bird location: start with lead bird named “bird_0”, heading diagonally NE start in SW corner – birds are horizontally arrayed with spacing of 10 Note that all birds have the same information processing algorithm – they cannot distinguish among themselves (in particular, the lead bird is not explicitly known to the others. Motion Plot 15 draws the current location of all birds, the leader’s color is red. The number, 15, indicates the threshold class flockSpace(..)allows simulations with different thresholds and a different number of birds. The values can be set calling the constructor.

21. birdPos birdPos birdPos birdPos birdPos birdPos Define an experimental frame to measure flock coherence and lead following public class centerTrans extends realDevs{ protected centerFn cs; public void deltext(double e,message x){ Continue(e); if (phaseIs("passive")) for (int i = 0; i < x.getLength(); i++) if (messageOnPort(x, "inOldNewPos", i)) { Pair p = (Pair)x.getValOnPort("inOldNewPos", i); cs.addVoldVnew(p); holdIn("output",0); } else if (messageOnPort(x, "inLead", i)) { entity en = x.getValOnPort("inLead", i); vect2DEnt leadPos = vect2DEnt.toObject(en); leadDist = leadPos.subtract(cs.center()).norm(); passivate(); }} Transducer to measure coherence and lead following lead birdPos center of mass distance to lead bird class centerFn extends Queue{ public vect2DEnt sum,avg; public centerFn(int numElements){ super(); …. } public void addVoldVnew(vect2DEnt vold ,vect2DEnt vnew){ remove(vold); add(vnew); if (size()>numElements)remove(); computeSum(); computeMax(); } enclosing radius center of mass (w/o lead bird) = sum of all position vectors /number of birds enclosingRadius = max distance of all birds from center distanceToLead bird = |lead position – center| measure of coherence+close following = w1*enclosingRadius +w2*distanceToLead Bird where w1,w2 can be chosen to weight coherence vs close following

22. Results of some representative threshold settings threshold = 0 threshold = 60 threshold = 22 enclosing radius for large threshold, birds can pick up the lead bird’s signal only after large distrance. distance to lead bird for small threshold, birds stay close together but lose the signal of the lead bird’s direction

23. Parallel Simulation Coupled Model for Optimization Study equally weighted sum of radius and distance to lead bird class flockSpaceOpt(..)allows simulations of flockSpaces with different thresholds at the same time and displays the values of the means (the current and the total). Breadth of Optimum Run alternative experiments in parallel Optimum Threshold = 25