Relative kinematics of an anchorless network

Re vised manuscript 1 Relati v e Kinematic s of an Anchorless Network Raj Thilak Rajan, Geert Leus and Alle-Jan van der V een Abstract Estimating the location of N coordinates in a P dimensional Euclidean space from pairwise distances (or proximity mea- surements), is a principal challenge i n a wide variety of fields. Con ven tionally , when localizing a static netw ork of immobile nodes, n on-linear dimension al reduction techniques are applied on the measured Euclidea n distan ce matrix (ED M) to obtain the relativ e coord inates upto a rotation and translation. In this article, we focus on an anchorless network of mobile nodes, where the distance measurements between the mobile nodes are time-varying in nature. Furthermore, in an anchorless network the absolute kno wledge of any node positions, motion or reference frame is absent. W e deri ve a nov el data mode l which relates t he ti me- v arying EDMs to the time-vary ing relative positions of an anchorless network. Using this data model, we estimate the relative position, relativ e v elocity and h igher o rder deriv ati ves, which are collectiv ely termed a s the relati ve kinematics of t he anchorless network. T he derived data model i s inherently ill - posed, howe v er can b e solved using certain r elat ive immobility constraints. W e propose elegan t closed form solution s to recursi vely estimate the relati ve kinematics of the network. For the sak e of completeness, estimators are also proposed to fi nd the absolute kinematics of the nodes, gi ven kno wn reference anchors. Cram ´ er-Rao bou nds are deriv ed for t he new data m odel an d simu lations are performed to analyze the pe rformance of the proposed solutions. Ke yword s: L y apuno v-like eq uation, relati ve velocity , relativ e acceleration, mu ltidimensional sca ling, time-v arying distance. I . I N T RO D U C T I O N Estimating the relative coord in ates of N po ints (or nodes) in a P d imensional E uclidean spa c e using prox imity measure ments is a f undam e ntal problem spanning a br oad range of applications. These ap plications include, but are not limited to, psy cho- metric analy sis [2], perceptu a l mapping [3], rang e-based anchorless localization [4], co mbinato r ial-chemistry [5], polar-based navigation [6], sensor arr ay calibration [7] and in general explorator y data analysis [8]. In anchorless localization scenarios for instance, n odes heavily rely on c o -oper a ti ve estimation of relative co ordina te s. Such ancho rless networks naturally arise when nodes a r e inaccessible o r only intermittently monitored , as is the case in space-based satellite arrays [9], un derwater networks [10] or indoor wireless sen sor networks [11]. In such reference-f ree scenarios, the prox imity inf ormation , often measured as pairwise d istances between the no des, form a key inpu t in e stima tin g the relativ e coo rdinates of n odes. These relativ e coor dinates are typically estimated using Non- linear dimensaion ality reduction algorithms (such as Multidimension al scaling (MDS)), which ha ve been s tudied rigoro usly over the past decades [ 8], [12]. Howe ver , con siderably less attention has been directed tow ards anchor less mobile s cenario s. Our prim ary focus in this article is on an anch orless n etwork o f mobile nodes, where we u se the term an chorless to in dicate no absolute knowledge of the nod e positions, motion or r eference frame. Furtherm ore, since the nod es are mo bile, both the node positions and the pairwise distance measuremen ts b etween the nodes are time-varying in n ature. Our motiv e is to relate the time-varying pairwise distance measurements to time-deriv ati ves of the node co ordina te s. For an an chorless network, these include the relative position , relative velocity , relative acceleration and higher-order derivati ves which we cumu lati vely refer to as relative kinematics in this article. It is worth n oting that the u niv ersal definitio n of relati ve kinematics inherently relies on the infor m ation in th e absolu te referen ce frame. F or example, the no n -relativistic relati ve velocity between two objects is rightly defined as the dif ference between their respective ab solute velocity vectors [13]. I n an anchorless frame work h owe ver , a natural question arises on whether th e relati ve kin ematics can be estimated , gi ven o nly time-varying distan ce mea surements. Ergo, we wish to understand the relatio nship between th e time-varying distance measurements and th e relati ve kinematics of mobile node s, which is the prime foc u s of this article. A. Pr evious work A key challenge in our pursuit is that both the time-varying distance and the time-varying relative po sitions are no n-linear in nature. In particular, the Euclidean d istance b e tween a pair of mobile nod es is almost always a non- linear fu n ction of time, ev en if the nodes are in linear independ ent motion [14]. Theref ore, it is p erhaps not surprising that traditional methods to solving such a pro blem have been state-spac e based appro aches with th e assistance of known anch ors [ 15]. T h e in itial p osition of the n odes are estimated using MDS like algorithm s, which u se the E uclidean d istance m atrix (EDM) at a single time-instant to estimate th e relati ve node positions. Given this initial estimate, the relative positions are tracked over a per io d of time with Doppler measur ements and kn own ancho rs [1 6], or via subspa ce tracking me th ods [17]. Unfortunately , Dop pler measurem ents and ancho r inf ormation are no t always a vailable. Secondly , subspace tracking is applicable only for small perturb ations in motion and therefor e of fer little insight on the kinematics of the motion itself. Submitted : 27th July 2016 , In re vision A part of this work is published in [1] R.T . Rajan, G.J.T . Leus and A.-J. van der V een are with TU D elft, Delft, The Nethe rlands (email: r . t.rajan@t udelft.nl; g.j.t.leus@tudelft.nl; a.j.v anderv een@tude lft.nl) Re vised manuscript 2 T ABLE I: Notations Notation Description P Number of dimension s N Number of nodes ( N > P ) D ( t ) ∈ R N × N Euclidean distance matrix at t ime t S ( t ) ∈ R P × N Absolute positions at time t S ( t ) ∈ R P × N Relativ e positions at time t X ∈ R P × N Absolute instantaneous positions at time t 0 X ∈ R P × N Relativ e instantaneous positions at t ime t 0 Y m ∈ R P × N m th order absolu te kinema tics at t 0 Y m ∈ R P × N m th order relati ve kinematics at t 0 H m ∈ R P × P Rotation matrix of the m t h order kin ematics h m ∈ R P × 1 T ranslational vector o f t he m th order kine matics In ou r p revious study , we p r oposed a two-step solution to estimate relative velocities of the no d es from time-varying distan c e measuremen ts [18]. Firstly , the d e r iv atives of the time-varying distances were estimated by solving a V ande r monde -like system of linear equations. The estimated regression coef ficients (called range parameter s) jointly yield the relati ve velocities an d the relativ e positions, u sing MDS-like algo rithms. Howe ver, the proposed solution is v alid o nly for lin ear motion , which is not always practical. Fur thermor e , the pr eviously p roposed MDS-ba sed relative velocity estimator hea vily relies on th e seco nd- order time-derivati ve of distance, and und er Gaussian n o ise assump tions, it performs worse th a n th e relative position estimator . Thus, desig n ing mo re optima l estimators f or the relati ve veloc ity is on e of the key motiv ations for th e pu rsuit of a generalized framework p resented in this article. Moreover, understand ing the higher order relative kinematics o f m otion in Euclidean space via time-varying d istance measuremen ts is crucial for next-generation localization tec h nolog ies. B. Contributions and overvie w W e presen t a n ovel data model in Section II , wh ich r elates the time-varying distanc e s to the kinematics of the mo bile nod es. More concretely , this relationship is established via the d eriv ativ es of th e time- varying distan c e (called range parameters), wh ich is estimated in Section III using d y namic ranging. In Section IV we show that the relationship between the range parameter s and the relati ve kinema tics takes the form of a L yapu nov-like set o f equations, wh ich is in herently ill-p osed. In p u rsuit of unique solutions, we propose elegant least squares algorithm s, which c an be solved un der certain assumptions. For the sake of comp letion, in Section V, we a lso propo se similar alg orithms for estimating the ab so lute kinematics of the nodes, g iv en known refe rence parameters in the cluster . T o compar e the per forman ce of our estimators, we deriv e constrained C ram ´ er-Rao bound s (CRBs), under Gaussian n oise assump tion on the data. A optimal choice of the weigh ting matrix en sures the proposed estimator is th e best lin ear unbiased estimato r (BLUE) f o r th e gi ven d ata model. In addition , unconstraine d o racle bounds are also deri ved in Section VI, a s a benchm ark for next g eneration estimator s. In Section VII , we conduct e xperimen ts to validate the perfo rmance of the proposed estimators. C. Notation: The element-wise matrix Hadamard produ ct is denoted by ⊙ and ( · ) ⊙ N denotes elem e nt-wise matrix exponent. The Kron ecker produ ct is indicated b y ⊗ , the tran spose o perator b y ( · ) T and ˆ ( · ) denotes an estimated value. 1 N ∈ R N × 1 is a vector o f ones, I N is an N × N identity matrix, 0 M ,N is an M × N matr ix of ze r os and k · k is the Euclidea n norm. For any vector a , diag ( a ) rep resents a diagonal matrix with a on the primary diag onal. For a diagonal matrix A , d iag ( A ) represents a vector o f approp riate len gth, co n taining the diago nal elements of the matrix A . Th e block diagon al matrix A = bdiag ( A 1 , A 2 , . . . , A N ) consists of m atrices A 1 , A 2 , . . . , A N along the diagonal a n d zeros else where. The first and second deriv ativ es ar e indica ted by ˙ ( · ) and ¨ ( · ) re spectiv ely , and more gen erally the m th order deriv ativ e is represen ted by ( · ) ( m ) . Unless oth erwise noted , ( · ) is used to ind icate parameters of the relati ve k inematic model. For matrices o f compatible d imensions, we will fre q uently u se the following p roper ties vec ( ABC ) =  C T ⊗ A  vec  B  , (1) vec ( A ) = J vec ( A T ) , (2) where J is an orth o gonal permutatio n m atrix. W e define an N dim ensional centerin g m atrix as P = I N − N − 1 1 N 1 T N . A brief list of frequ ently used notations are tab u lated in T ABLE I. Re vised manuscript 3 I I . T I M E - V A RY I N G D I S TA N C E S A N D N O D E K I N E M A T I C S W e begin by mod eling the r e la tio nship betwe e n the time-varying distances, the time-varying po sitions and the nod e kin e- matics. In Section I I-A, we expan d th e time-varying p osition using a T ay lor series, the coe fficients of wh ich yield th e absolute node kinematics. As an extensio n, we present a novel relative kinematics model in Section II-B. In Sections II-C and II-D, the relationship between the time-varying distances a n d the node kinematics is deriv ed. Using these definitions, we formalize the prob lem statem e nt in Section II-E. A. Absolute kinematics Consider a cluster of N mob ile nodes in a P dime n sional Euclid ean space ( N > P ), whose positions at time t are given by S ( t ) ∈ R P × N . For a small time interval ∆ t = t − t 0 around t 0 , we assume the time-varying position is continuou sly differentiable M times an d th at the M th derivati ve exists in the interior of this inter val. T h erefor e , the time-dep e ndent position vectors of the respecti ve nodes can be expanded using a T ay lor series, S ( t ) = S ( t ) | t = t 0 + ˙ S ( t ) | t = t 0 ( t − t 0 ) + ¨ S ( t ) | t = t 0 ( t − t 0 ) 2 + . . . (3) where ( S ( t ) , ˙ S ( t ) , ¨ S ( t ) , . . . ) a re the deriv ativ es o f the time-varying position vectors. Now let X , S ( t ) | t = t 0 be a P × N matrix containin g the in itial coordin ates of the mobile n odes at time t = t 0 . Furtherm ore, let th e instantaneo us velocities of the nodes i.e., the first-order deriv atives of the position vectors ˙ S ( t ) | t = t 0 be denoted by Y 1 ∈ R P × N , and in general higher-order deriv ativ es as Y m ∀ 1 ≤ m ≤ M . Then , the above equ ation simplifies to S ( t ) = X + M X m =1 ( m !) − 1 Y m ( t − t 0 ) m . (4) B. Relative kinematics The absolute instantaneo us positions at t = t 0 are an affine transformation of the relati ve positions, i.e., X = H 0 X + h 0 1 T N , (5) where X ∈ R P × N is the r elativ e position matrix upto a rotation and tran slatio n , H 0 ∈ R P × P is the unk nown rotation and h 0 ∈ R P × 1 is the unknown translation of the network [8]. No w , we extend this well-known relative position definition to the higher-order deri vatives . F or instance, the velocity of the nod es c an be written as Y 1 = H 1 ˜ Y 1 + h 1 1 T N , (6) where ˜ Y 1 represents the instantaneo us relati ve velocities of the n etwork at t = t 0 . The tr a nslational vector h 1 is the gro up velocity and H 1 is the u nique rotation matrix of the relativ e velocities [18]. Mor e g enerally , the m th o r der deriv ativ e is an affine model defined as Y m = H m ˜ Y m + h m 1 T N . (7) W e no w define the r elativ e time-varying position as S ( t ) = H T 0 S ( t ) P , and substituting the affine expre ssions (5) and (7) in (4) we have S ( t ) = XP + M X m =1 ( m !) − 1 H T 0 H m ˜ Y m P ( t − t 0 ) m , (8) where we exploit the property P1 N = 0 N to eliminate the translation vectors, and enforce the orthono rmality of the rotation matrix i.e., H T 0 H 0 = I N . Ob ser ve that th e tran slation vector h 0 does not affect the above equation . Secondly , for a mean in gful interpretatio n of the relative time-varying p osition, a refer e nce coor dinate system must be chosen e.g., H 0 = I P . In summ ary , without the loss of gen e r ality , we assume H 0 = I P and h 0 = 0 P . (9) and subsequen tly (8) simp lifies to S ( t ) = X + M X m =1 ( m !) − 1 Y m ( t − t 0 ) m , (10) where we use the following properties X = XP = XP , (11a) Y m = H m ˜ Y m = Y m P , (11b) S ( t ) = S ( t ) P . (11c) Re vised manuscript 4 Note that (10) repr esents the relative counterpa r t o f the absolute T aylor expansion (4), wher e the  X , Y 1 , Y 2 , . . . , Y M  denote the relative k inematics of the correspon d ing abso lute kinematics  X , Y 1 , Y 2 , . . . , Y M  . Our quest in this article is to estimate the relative and absolute kinem atic matrices, gi ven time-varying pairwise distance measur ements between the nodes. Consequently , the absolute position S ( t ) and relativ e position S ( t ) can then be estimated using (4 ) and (10) respectively . C. T ime-varying distances Similar to the nod e positions, the pairwise distances are also time-varying which we denote by the time-varying E u clidean distance matrix (EDM) D ( t ) , [ d ij ( t )] ∈ R N × N where d ij ( t ) is the pairwise Euclidean distan c e between th e no de pair ( i, j ) at time instant t . More explicitly  D ( t )  ⊙ 2 = ζ ( t ) 1 T N + 1 N ζ T ( t ) − 2 S T ( t ) S ( t ) , (12) where ζ ( t ) = diag  S T ( t ) S ( t )  . Observe th at D ( t ) is a non-linear fun ction of time t , e ven when the nodes a re in ind e p enden t linear m otion an d h ence D ( t ) is a continuo usly differentiable function in time. No w , b a sed on the tim e-varying EDM D ( t ) , we define the doub le centered matrix B ( t ) B ( t ) , − 0 . 5 P  D ( t )  ⊙ 2 P , (13a) and the time deriv atives of the dou ble centered matrix  ˙ B ( t ) , ¨ B ( t )  for upto M = 2 as, ˙ B ( t ) , − P  D ( t ) ⊙ ˙ D ( t )  P , (13b) ¨ B ( t ) , − P  D ( t ) ⊙ ¨ D ( t ) + ( ˙ D ( t )) ⊙ 2  P , (13c) where  ˙ D ( t ) , ¨ D ( t ) , . . .  are th e de riv atives of the time -varying EDM , wh ic h in d icate th e rad ial velocity and o th er h igher-order deriv ativ es. Now , let the EDM and the correspond ing deriv atives at t = t 0 be d enoted b y D ( t ) | t = t 0 , R = [ r ij ] , ˙ D ( t ) | t = t 0 , ˙ R = [ ˙ r ij ] , ¨ D ( t ) | t = t 0 , ¨ R = [ ¨ r ij ] , ∀{ i, j } ≤ N , then with an abuse of notation (13) becomes B (0) , B ( t ) | t = t 0 = − 0 . 5 PR ⊙ 2 P , (14a) B (1) , ˙ B ( t ) | t = t 0 = − P h R ⊙ ˙ R i P , (14b) B (2) , ¨ B ( t ) | t = t 0 = − P h R ⊙ ¨ R + ˙ R ⊙ 2 i P , (14c) and higher-order deriv atives can be defined along similar line s. In general, giv en the distance derivati ves at t 0 , i.e., the range parameters ( R , ˙ R , ¨ R , . . . ) , the do uble center ed ma tr ix B (0) and th e correspondin g higher-order deriv atives ( B (1) , B (2) , . . . ) can be co nstructed. In a mo bile network, the range par ameters may not be r eadily av ailab le, howe ver gi ven all the n odes are capable of two way ranging, the range parameter s can be estimated using dyn a mic ranging [14]. D. Model T o und e r stand th e relationship between the time- varying distances and the relati ve kinematics of the nodes, we sub stitute the definition of the EDM from (1 2) in (13a) and differentiate recu rsiv ely to obtain B ( t ) = S T ( t ) S ( t ) , (15a) ˙ B ( t ) = ˙ S T ( t ) S ( t ) + S T ( t ) ˙ S ( t ) , (1 5b) ¨ B ( t ) = S T ( t ) ¨ S ( t ) + ¨ S T ( t ) S ( t ) + 2 ˙ S T ( t ) ˙ S ( t ) , (15c) where we use the d e finition (11c) and introd uce ( ˙ S ( t ) , ¨ S ( t ) , . . . ) as the der iv atives of S ( t ) . Now , r earrang ing the term s and substituting the definition of S ( t ) at t = t 0 from (10) , we ha ve B 0 , B (0) = X T X , (16a) B 1 , B (1) = X T Y 1 + Y T 1 X , (16b) B 2 , B (2) − 2 Y T 1 Y 1 = X T Y 2 + Y T 2 X , (16c) where we in troduce the m atrices ( B 0 , B 1 , B 2 ) . The joint left and right centering using the centering matrix P in (13) ensu res that the phase center of the relative kinematic matrices ( Y 1 , Y 2 ) are at 0 P , similar to the relative po sition X . 1) Relative kinematics: Now , for M = 0 , comb ining (14a) and (16a), we ha ve B 0 = X T X = − 0 . 5 PR ⊙ 2 P , (17) Re vised manuscript 5 and more gener ally for a gi ven M ≥ 1 , (16 ) can be generalized to B M , B ( M ) − M − 1 X m =1  M − 1 m  Y T M − m Y m (18a) = X T Y M + Y T M X , (18b) where B ( M ) is the M th deriv ative of th e double centere d matr ix at t 0 , which is given by (14) an d Y M is the M th order relativ e kinematic matrix. Remark 1 : ( Measurement matrix B M ): W e make two critical observations on B M in (18a) . • Firstly , no te that B M is depen dent on the range parameters ( R , ˙ R , ¨ R , . . . ) via the definition of B ( M ) (14). • Secondly , B 0 , B (0) and B 1 , B (1) can be constructed only based o n the range parameter s (see (16 )). Howe ver for M ≥ 2 , B M in ad d ition to B ( M ) , add itionally relies o n the relati ve kinematic matrices o f order less than M . Hence, if the lower o rder kinema tics Y m ∀ 2 ≤ m < M are kn own, then the measu r ement matrix B M can be recon structed. 2) Absolute kinematics: In addition to the relative k inematics, (1 8b) can also b e reformula te d to estimate the absolu te kinematics Y M of the network . Recall from (11 b), that the relative kinematics of the M th order is Y M = Y M P under the assumption (9). Substituting this expression in (18b), we hav e B M = X T Y M P + PY T M X , (19) which is the absolute kinem atic model. 3) Model summary: I n summar y , if the r ange parameters ( R , ˙ R , ¨ R , . . . ) are av ailable, B ( M ) can be co nstructed from (1 4). Giv en B (0) , we aim to so lve for the relative position X using the equation (1 7), which we use to estimate the higher order kinematics. For M ≥ 1 , the measurement matrix B M can be constructed using B ( M ) and by substituting the lower order relativ e k inematic matrices Y m ∀ 2 ≤ m < M in (18 a). Finally , given the measuremen t matrix B M and an estimate o f X , our go al is to estimate the M th ord er relativ e kinem a tics Y M and the absolute kinem atics Y M for M ≥ 1 , u sing (18 b) and (19) respectively . W e n ow formulate the prob lem m ore concretely in the following section. E. Pr ob lem Sta te m e n t Problem statement: Gi ven the time-varying pairwise d istances D ( t ) between the N n odes in a P dimension a l Euclidean space, estimate the r elativ e kinematics ( X , Y 1 , Y 2 . . . ) and a b solute kinematics ( Y 1 , Y 2 . . . ) of the m obile network. Th ese estimates subsequen tly yield the relati ve (and absolu te) time-varying positions. Solution: W e propo se a tw o-step solution to the above estimation problem. S1) Dyn amic ranging an d r elative position : Giv en the time-varying distance measurem ents D ( t ) , we employ dynam ic ran ging to obtain the ran ge parameters ( R , ˙ R , ¨ R , . . . ) in Section III, un der the assumption that all the nodes are capable of commun icating with ea ch other . Secondly , we a lso estimate the initial relati ve position X using (17). S2) Kinema tics : Th e m easurement ma tr ix B M can be con structed using the estimated range parameters, and lower ord er kinematics (18a). Gi ven the relati ve position X and B M estimates, we solve for the relati ve kinema tics Y M (in Section IV), and the absolute kinematics Y M (in Section V), using (18 b) and (19) respecti vely . Finally , given the initial relative position and the node kinem atics, the time-varying absolute an d relative positions { S ( t ) , S ( t ) } can be estimated using (4) and (1 0 ) respecti vely . I I I . D Y N A M I C R A N G I N G A N D R E L A T I V E P O S I T I O N In this section , we a im to estimate the ran ge par ameters ( R , ˙ R , ¨ R , . . . ) , given two-way communicatio n b e tween the n odes in the mobile network. In Section III-A, we relate the time-varying propaga tion delay between the nodes and th e ran ge p arameters. Giv en this r elationship, w e present a Dynamic r anging mod el in Sectio n II I-B, and su b sequently pre sen t a closed form algorithm to estimate the range par ameters in Section III-C. Finally , we app ly the MDS algorithm to find the initial relative position of the nod e s in Section III-D. A. T ime -varying pr o pagation delay Consider a p a ir of mobile no des capable of comm u nicating with each othe r . Let τ ij ( t 0 ) , τ j i ( t 0 ) = c − 1 d ij ( t 0 ) be the propag ation delay of this comm unication between the n ode pair ( i, j ) at time instant t 0 , where d ij ( t 0 ) is the corresp onding pairwise d istance an d c is the spe e d of the electromagn etic wa ve in the m edium. Now , f or a small interval ∆ t = t − t 0 , we assume the relati ve distance to be a smooth ly varying polynom ia l of ti me which enables us to describe the p ropaga tio n delay τ ij ( t ) at t as an infinite T aylor series in th e neighbor hood of t 0 τ ij ( t ) = c − 1 d ij ( t ) = r ij + ˙ r ij ( t − t 0 ) + ¨ r ij ( t − t 0 ) 2 + . . . , (20) Re vised manuscript 6 Node j T j i, 1 T j i, 2 T j i, 3 T j i, 4 T j i,K Node i T ij, 1 T ij, 2 T ij, 3 T ij, 4 T ij,K Ti me t → Distance d ij ( t ) → Fig. 1: Dynamic ra nging: A Generalized T wo W ay Ranging (GTWR) scenario between a pair of mobile nodes, where t he nodes exchange K time stamps asymmetrically with each other [14]. The curved lines symbolize the non-linear motion of the mobile nodes with time. Unlike our prev ious models [19], [18] which considered only linear independe nt velocities of the nodes, in this article we consider non-linear motion of the nodes. where the T aylor coefficients are defined as  r ij , ˙ r ij , ¨ r ij , . . .  T = diag ( γ ) − 1  r ij , ˙ r ij , ¨ r ij , . . .  T , (21) and γ = c [0! , 1! , 2! , . . . ] T . Here, ( r ij , ˙ r ij , ¨ r ij , . . . ) are the deriv atives of th e tim e-varying p airwise distanc e d ij ( t ) esimtated at t = t 0 , which are the eleme nts of the ma tr ices ( R , ˙ R , ¨ R , . . . ) , pr esented earlier in Section II-C. The ph ysical significance of these co efficients is as follows. The p airwise distance at t 0 is r ij , which is conv ention a lly obtained fro m time o f arriv al measuremen ts. ˙ r ij is the radial velocity , typica lly observed from Dop pler shifts, an d the second-ord e r range parameter ¨ r ij is the rate of rad ial velocity b etween the node pair at t 0 . W e will now use this relation in a scen ario where mobile nodes are capable of two-way communicatio n . B. Data model Consider a Gen eralized T wo W ay Rang ing ( GTWR) scenario b etween a p air of m o bile nodes (Fig . 1 ), where the n odes commun icate asym metrically with each other, and r e cord K timestamp s on each node. Th e timestamps recor ded at the k th time instant ( k < K ) at node i and node j a r e giv en by T ij,k and T j i,k respectively . The nodes are mobile during these timestamp exchanges, and therefore th e propagation delay between the no des is un ique at ev ery time instant. W ith an abuse of n otation, let τ ij,k and d ij,k be th e propag ation delay and th e distanc e between the node pair ( i, j ) at th e k th time instant. Then assuming the distance is (ap prox) con stant during the pro pagation time of the message,the non-relativistic prop agation delay is τ ij,k = c − 1 d ij,k = | T ij,k − T j i,k | . No w , observe that the pairwise propagation delay for GTWR can also be written as (20), by replacing t with T ij,k (or T j i,k ). More concretely , the propagation d e lay τ ij is giv en as τ ij,k = | T ij,k − T j i,k | = r ij + ˙ r ij ( T ij,k − T 0 ) + ¨ r ij ( T ij,k − T 0 ) 2 + . . . , (22) where the range param eters are estimated at T 0 where T ij,k ≤ T 0 ≤ T ij,K . Aggregating all the K timestamps for eac h no d e pair ( i, j ) , and po pulating all measureme n ts f r om ¯ N , 0 . 5 N ( N − 1) unique pairwise links for a network of N nodes, we have V z }| {  I ¯ N ⊗ 1 K T T ⊙ 2 . . .  θ z }| {      r ˙ r ¨ r . . .      = τ , (23) where for an L th o r der po lynomial ap prox im ation, θ ∈ R ¯ N L × 1 is a vector o f unkn own coefficients. T he ¯ N dim ensional vector r = [ r ij ] ∀ 1 ≤ i ≤ N , j ≤ i con tains all the p airwise distan ces at t 0 , and vecto r s containin g the higher or der deriv atives ( ˙ r , ¨ r , . . . ) are similarly defined . The matrix V is a V andermo nde-like matrix de fined as V = [ I ¯ N ⊗ 1 K T T ⊙ 2 . . . ] ∈ R ¯ N K × ¯ N L , where T = bdiag ( t 12 , t 13 , . . . t 1 N , t 23 , . . . ) ∈ R ¯ N K × ¯ N and t ij = [ T ij, 1 − T 0 , T ij, 2 − T 0 , . . . , T ij,K − T 0 ] T ∈ R K × 1 contain all the time stamps. All the u nique pairwise propagation delays are collected in τ = [ τ T 12 , τ T 13 , . . . τ T 1 N , τ T 23 , . . . ] T ∈ R N K × 1 where τ ij = | t j i − t ij | . Our goal in the following sectio n, is to estimate the values  r ij , ˙ r ij , ¨ r ij , . . .  from (23), whic h will help us construc t the range matrices  R , ˙ R , ¨ R , . . .  . Re vised manuscript 7 C. Dynamic ranging algorithm In reality , the propag ation delay is erroneous and hence, more practically (23) is ˆ τ = V θ + η , (24) where ˆ τ is the noisy propa g ation delay , and th e noise parameter s p la g uing the data mod el ar e populated in η = [ η T 12 , η T 13 , . . . η T 1 N , η T 23 , . . . ] T ∈ R ¯ N K × 1 , wher e η ij = [ η ij, 1 , η ij, 2 , . . . , η ij,K ] is th e error un iq ue to the node pair ( i, j ) . I n practice, the noise ar e on the time markers T ij,k and subsequently on the V ander monde matrix, wh ich has been simplified under no minal assumptions to arrive at th e elegant m odel (24). The appro x imations in volved are d iscussed in Append ix A. Now , suppose the covariance of the no ise on the normal equations Σ , E  η η T  , (25) is known a n d in vertible, then the weighted least squares solution ˆ θ is obtained by minimizing the following l 2 norm, ˆ θ = arg min θ k Σ − 1 / 2 ( V θ − ˆ τ ) k 2 =  V T Σ − 1 V  − 1 V T Σ − 1 ˆ τ . (26) A valid solution is f easible if K ≥ L fo r each of the ¯ N pairw ise lin ks. Mo re gen erally , whe n L is un known, an or der rec u rsiv e least squar es can be e mployed to ob tain the ra n ge coefficients [18]. Given θ , estimates o f the range param eter matrices  ˆ R , ˆ ˙ R , ˆ ¨ R , . . .  can be constructed using (21) and subsequen tly , f rom (14) we have th e following estimates ˆ B (0) = − 0 . 5 P ˆ R ⊙ 2 P , (27a) ˆ B (1) = − P h ˆ R ⊙ ˆ ˙ R i P , (27b) ˆ B (2) = − P h ˆ R ⊙ ˆ ¨ R + ˆ ˙ R ⊙ 2 i P . (27c) D. Relative position Giv e th e initial pairwise d istances at t 0 i.e., R , the initial r elativ e p ositions X can b e determ in ed via M DS. Gi ven ˆ R , let ˆ B 0 be an estimate of B 0 , B (0) , obtained using (27a). A spectral decom position o f this matrix yields ˆ B 0 = V x Λ x V T x , where Λ x is an N d imensional diago nal matrix containin g the eigenv alu e s of the ˆ B 0 and V x the cor respond ing eigenvectors. An estimate of the relative po sition estimate using MDS is then given by ˆ X = arg min X k ˆ B 0 − X T X k s.t. rank ( X ) = P = Λ 1 / 2 x V T x , (28) where Λ x contains the fir st P nonz e ro eigenv alues fro m Λ x and V x is a subset of V x containing the cor r espondin g eigenv ectors [8]. I V . R E L A T I V E K I N E M AT I C S In the pre vious section, we estimated the range parameters given time-varying distance me a surements D ( t ) , which was the first step (S1) in our problem statement described in Section I I-E. Using these range parameters, we constructed the d ouble centered matrices  ˆ B (0) , ˆ B (1) , ˆ B (2) , . . .  (27) and estimated the relative position ˆ X using MDS (28). Given these estimates, we now aim to solve the unknown relative kin ematic matrices Y M using (18), as pro posed in (S2) of Section II-E. A. Linearized multidimension a l scaling (LMDS) Prior to in vestigating the genera l kinem atic mo del (18), we re visit a special case when th e nodes are mobile under linear indepen d ent m otion [1 8]. In such a scen ario, the accelera tion and other higher o rder de riv atives are absent i. e ., Y m = 0 , ∀ m ≥ 2 . Therefo r e, u nder constant velocity assumption, equations (16b) and (16c) simplify to B (1) = X T Y 1 + Y T 1 X , (29a) B (2) = 2 Y T 1 Y 1 , ( 29b) and for m ≥ 3 { B m , B ( m ) } defined in (1 8) d oes not exist [18, Appendix B]. Now sub stitutin g the definition of relative velocity from (11b) and exploiting the property H T 1 H 1 = I , we have B (1) = X T H 1 ˜ Y 1 + ˜ Y T 1 H T 1 X , (30a) B (2) = 2 ˜ Y T 1 ˜ Y 1 . (30b) Re vised manuscript 8 The LMDS algorith m to es timate the r elativ e velocity (upto a translation) is then a two s tep meth od as decribed below . 1) MDS-ba sed r ela tive velo city e stimator: Firstly , th e r elativ e velocity up to a rotation and tran slation is obtained by minimizing the strain functio n using (3 0b). Let ˆ B (2) be an estimate of B (2) from (2 7c), with an eigen value d ecompo sition ˆ B (2) , V y Λ y V T y , then the relative velocity estimate is given b y ˆ ˜ Y 1 = arg min ˜ Y 1 k ˆ B (2) − 2 ˜ Y T 1 ˜ Y 1 k s.t. rank ( ˜ Y 1 ) = P = Λ 1 / 2 y V T y , (31) where Λ y and V y contain the first P nonzero eigenv a lu es and co r respond ing eigenvectors of Λ y and V y respectively . 2) Estimating th e un k nown r o tation: The MDS-based solution (31) yield s the relativ e velocity up to a r otation an d translation , which is not suffi cient to re c o nstruct the time-varying relative position using (8). T o estimate the un iq ue rotatio n ma tr ix, we vectorize (30a), apply the transformation (1), and solve the following con strained cost function arg min H 1 k ˆ Φ vec ( H 1 ) − vec ( ˆ B (1) ) k 2 s.t H T 1 H 1 = I P , (32) where ˆ Φ = ( I N 2 + J )( ˆ ˜ Y T 1 ⊗ ˆ X T ) , { ˆ X , ˆ ˜ Y 1 } ar e estimates obtained fro m (28) and (31) respectively and, J is a pe rmutation matrix such that (2) ho lds. Thus, u nder line a r velocity assump tion, the relative velocity Y 1 = H 1 ˜ Y 1 up to a translation c a n be rec onstructed for a general P d imensional scenario using the estimators (3 1 ) and (32). It is worth notin g that the LMDS solution is feasible, only under the constan t velocity assum ption. I n gener al, the assumption on linear motion is not always valid a n d h ence we addr e ss the more general kinem atic motion in the following sections. B. Lyapunov-like equatio ns More generally , when the n odes are in non-linear motion, the kinematics Y m ∀ m ≥ 1 exists and must be estimated. T o solve for the r elativ e kin e matics in this scenario , we refer back to our relative kinematic model (18 ) . For any M ≥ 1 , the model (18b ) B M = X T Y M + Y T M X , (33) is the relati ve L yapunov-like equ ation [20], [21], where B M is the N − dimen sional measurement ma tr ix and Y M is the M th order kinematics to be estimated. As pointed o ut in Remark 1 in Section II-D, B M can be constructed by B ( M ) and lower order relativ e kinematics { Y m } M − 1 m =1 . The above equ ation is very similar , but not the same as the following equations, A H Y + Y A = B , A Y + Y A = 0 , A Y + YC = E , which are the (con tinuous) Lyapunov equation , commu tativity equa tion [22, c h apter 4] and Sylvester equation [2 3], [2 4] respectively , where the unkn own matrix Y has to be e stima te d , given A , B , C , E . The solutions to the se equations exist and are extensi vely in vestigated in control theory literature [25]. Howe ver the L y apunov- like equation (33) has recei ved relati vely less attention. The L yap unov-like equatio n ha s a straigh t forward solution for P = 1 . But, f or P ≥ 2 , althoug h a general solution was proposed by Braden [26], a unique solution to (3 3 ) does not e xist which we discu ss in App e ndix B. Now , vectorizing (33) and using (1), we aim to solve ˆ y M = arg min y M k ( I N 2 + J )( I N ⊗ X T ) y M − b M k 2 = arg min y M k A y M − b M k 2 , (34 ) where A = ( I N 2 + J )( I N ⊗ X T ) ∈ R N 2 × N P , (35a) y M = vec ( Y M ) ∈ R N P × 1 , (35b) b M = vec ( B M ) ∈ R N P × 1 , (35c) and J is an orthog onal permutation matrix (2). The m a tr ix ( I N ⊗ X T ) ∈ R N 2 × N P is f u ll colu mn rank, sinc e X is typically non-sing ular . Howe ver , the sum of p ermutatio n matrices ( I N 2 + J ) ∈ R N 2 × N 2 is always ran k deficient by at least  N 2  . Hence, the ma trix p rimary o bjective function A is not fu ll colu m n rank , but is rank de ficie n t by at least ¯ P , 0 . 5 P ( P − 1) , which is discussed in Appendix B. In (33), since the tran slatio n al vectors of both X and Y M are pro jec ted out u sin g the center ing matrix P , the ¯ P d epende n t columns in A indicate the rotational degrees of f reedom in a P -dim ensional Euclidean space. Re vised manuscript 9 C. Lyapunov-like least squares ( LLS) A unique solutio n to the L yapun ov-like equ ation is not feasible without sufficient constraints on the linear sy stem (34). Let ˆ A be an estimate of the A , obtained by substituting the estimated relative position ˆ X (28 ). Similarly , let ˆ b M be an estimate of b M obtained by substituting the ra n ge parameters and appropriate relati ve kinematic matrices u pto order M − 1 . Then the constrained L yapun ov-like least sq uares (LLS) so lu tion to estimating the relati ve kin ematic matrices is g iv en b y minimizing the cost function ˆ y M ,lls = arg min y M k ˆ Ay M − ˆ b M k 2 s.t. ¯ Cy M = ¯ d , (36) where ¯ C is a set of non-redu n dant constrain ts. The a bove o ptimization problem has a closed- form solution, giv en by solving the KKT equation s ( Append ix C). D. W eig h ted L yapu nov-like LS (WLLS) In reality , bo th A and b are plague d with err ors and hence th e solution to the cost fun c tion ( 36) is sub -optimal. Let ¯ W be an appro priate weighting m atrix on the L y apunov- like equation, then the W eigh ted L y apunov- like Least Squar es (WLLS) solution is obtained by minimizin g the cost function ˆ y M ,w lls = arg min y M k ¯ W 1 / 2 M  ˆ A y M − ˆ b M  k 2 s.t. ¯ Cy M = ¯ d , (37) which, similar to (3 6), can be solved using the constraine d KKT solutions (Ap pendix C). An appro priate choice of the weigh ting matrix ¯ W M will be discussed in Section VI-D. E. Choice of constraints: Relative immobility In th e absence of absolute location in formatio n, a unique solution is feasible if the relative motion of at lea st P nod es or features are inv a r iant (o r known) over a small time dur ation ∆ t . In an anchor le ss framework, a set of gi ven nodes would have equiv alen t relativ e kin ematics, if they are identical in motion upto a translation or if they are immobile for the small measuremen t time ∆ t . Such sit uation s could ar ise, for examp le, in underwater lo calization, when a few immobile no des could be fixed with unknown absolute location s, w h ich in turn could assist the r elativ e localization of the o th er nod es. For P = 2 , if the first P n odes are relatively imm obile for the small measurem e nt time, a v alid constraint for (36) and (37) is ¯ C 1 =  I 2 − I 2 0  , ¯ d 1 = 0 , (38) which ca n be readily extended for P > 2 and if require d , for a larger nu mber of immobile nodes. In essence, the relative immobility constraint redu c es th e parame te r space in pur su it of a unique solution for the ill-posed L y a punov-like equation. F . T ime-va rying r ela tive po sition In this section, we solved for the relativ e kinematics of motion, using the range parameters and relative po sition estimates. When the nodes are in linear motion, the first-or der relative kin ematics can be estimated using the LMDS algorith m ( 3 1, 32). More genera lly , for estimating the relati ve kinematics in a non-linear scenario, we solve th e L y a p unov-like equation (33) usin g constrained Least squares (36 , 37). Substitutin g these estimates in (10), an estimate of the relative time-varying position is ˆ S ( t ) = ˆ X + ˆ Y 1 ( t − t 0 ) + 0 . 5 ˆ Y 2 ( t − t 0 ) 2 + . . . (39) where ˆ X is a r e la tive po sition estimate fro m (2 8) an d { ˆ Y 1 , ˆ Y 2 , . . . } are the estimate s fro m (3 6) or (37). I n the fo llowing section, we aim to estimate the absolute kine m atics of the nodes and subsequently the time-varying a bsolute position. V . A B S O L U T E K I N E M A T I C S In this section, we solve for the absolute kinematics Y M , given B M and the relativ e position X . W e hav e from (19), X T Y M P + PY T M X = B M . (40) The above e q uation is similar , b ut not the same, to the genera lized (continuous-tim e) L yapu nov e q uation A T YC + C T Y A = B , where A , B , C are known square matrices [27]. W e now vectorize (40) and aim to minim ize the following co st fun c tio n ˆ y M = arg min y M k Ay M − b M k 2 , (41) Re vised manuscript 10 where A = ( I N 2 + J )( P ⊗ X T ) ∈ R N 2 × N P , (42a) y M = vec ( Y M ) ∈ R N P × 1 , (42b) and b M is given by (35 c). I n comparison to (34), the matrix ( I N ⊗ X T ) is replaced with ( P ⊗ X T ) in (42a). The rank of the centering ma tr ix P is N − 1 and since X is typically f ull row rank, the Kro necker pro duct is utmo st of rank N P − P . This rank-deficien cy of P is also reflected in the matrix A . Unlike A which has ¯ P dep endent c olomns, A is rank -deficient by  P +1 2  = ¯ P + P . The additional P d epende n t column s a r e perhaps not surprising, as they in dicate the lack of information on the translation al v ector, i.e., the gr oup center of the M th order kinematic matrix . A. Generalized L yapun ov-like least squ ar es (GLLS) In pursuit of a unique solution to the rank-deficient system (41), we propose a c onstrained generalized L y apunov- like least squares (GLLS) to estimate the absolu te kinematic matrices which is obtained by minimizing the cost function ˆ y M ,g lls = arg min y M k ˆ Ay M − ˆ b M k 2 s.t. Cy M = d , ( 43) where ˆ A and ˆ b M are estimates o f A and b M respectively . The matrix C is a set o f n on-red undan t constraints, which will be discussed in Section V -C. B. W eig hted generalized Lyapunov-like LS (WGLLS) The perfor mance of the es timator c a n b e imp r oved by weigh ting the cost function (43), i.e., ˆ y M ,w g lls = arg min y M k W 1 / 2 M  ˆ Ay M − ˆ b M  k 2 s.t. Cy M = d , (44) which yields the weighted generalized L ya p unov-like least sq u ares (WGLLS) solution (see Append ix C), where W M is an approp riate we ig hting matrix (see Section VI-D). C. Choice of con straints: Anchor-awar e network For an anchor ed scenario, if the M th o rder absolute kinem a tics of a few nodes a r e k nown, then the a bsolute velocity , acceleration and h igher ord er deriv a tives can be estimated. A straightforward minimal constra in t for the feasib le solution is then C 1 =  I ¯ P + P , 0  , (45 ) where without loss of gener ality , we a ssume the first ¯ P + P param eters are kno wn. D. T ime- varying absolu te position In (43, 44) , we solved for the absolute kinematics given the measurem ent matrix B M and the relative p osition, using constrained Least squares estimators. Given these estimates, w e ha ve from (4) ˆ S ( t ) = ˆ X + ˆ Y 1 ( t − t 0 ) + 0 . 5 ˆ Y 2 ( t − t 0 ) 2 + . . . , (46) where ˆ S ( t ) is an estimate of the tim e -varying ab solute p osition, ˆ X is an estimate of the relati ve position (28), an d { ˆ Y 1 , ˆ Y 2 , . . . } are the absolute kinem atic estimates obtained by solving (43) or (44). V I . C R A M ´ E R - R AO B O U N D S The Cra m ´ er-Rao lower Boun d (CRB) sets a lower b ound on the minimum achiev ab le variance of any unb iased estimato r . In this s ection, we d e riv e th e CRBs for the estimated par ameters based on the presented data models. I n the fo llowing section, we will use these boun ds to benchmark the performance of the propo sed estimato rs. A. Range parameters W e b egin by d eriving th e lower bou nds o n th e ran ge para m eters. L et ψ = [ r T , ˙ r T , ¨ r T , . . . ] T and let ˆ ψ be th e co rrespond ing estimate, then the CRB Σ ψ , E n ( ˆ ψ − ψ )( ˆ ψ − ψ ) T o on the range parame te r s for the lin ear model ( 24 ) is Σ ψ ≥ Γ  V T Σ − 1 V  − 1 Γ =      Σ r ∗ ∗ ∗ ∗ Σ ˙ r ∗ ∗ ∗ ∗ Σ ¨ r ∗ ∗ ∗ ∗ . . .      , (47 ) Re vised manuscript 11 where Σ ψ is the covariance of ψ and Σ is th e covariance of the no ise o n the timestamps d e fined in (25 ) . Her e, the co variance matrices { Σ r , Σ ˙ r , Σ ¨ r , . . . } are the lowest achiev able bo unds for the correspon ding rang e parameters { r , ˙ r , ¨ r , . . . } . The en tries not of interest a r e deno ted b y ∗ and Γ = diag ( γ ) ⊗ I ¯ N is a transformation ma tr ix, w h ere γ is g iv en by (2 1 ). It is w orth n o ting that our propo sed solution (26) achie ves th is lower bound for an appropriate L . B. Relative position The CRB on the relative p ositions y 0 , vec ( X ) is giv en by the in verse of the Fisher Informatio n Matrix (FIM) i.e., Σ x , E  ( ˆ y 0 − y 0 )( ˆ y 0 − y 0 ) T  ≥ F † x , (48) where ˆ y 0 is an estimate of th e unknown relative position y 0 , Σ x is th e co variance of ˆ y 0 [18] and the FIM F x ∈ R N P × N P is F x = J T x ¯ Σ − 1 r J x , (49) where ¯ Σ r , bdiag ( Σ r , Σ r ) , J x is the Jacobian [18, Appendix C] and Σ r is obtained from (47). In the absence of known anchor s in the network, the FIM is inheren tly nonlinear and hence we employ the Moore-Penro se pseud oinv erse in (4 8). C. Kinematics W e now derive th e lower bounds on th e v ariance of th e estimates of the r e la tive kinematics y M = v ec ( Y M ) an d absolute kinematics y M = vec ( Y M ) . The Gaussian noise vectors plaguin g the cost functions (34 ) and (41) are mo deled as ρ M ∼ N ( A y M − b M , Σ ρ,M ) , (50) ρ M ∼ N ( Ay M − b M , Σ ρ,M ) , (51) where ρ M , ρ M are N 2 dimensiona l noise vectors, an d the correspond ing covariance matrices are of the form Σ ρ,M , E  ρ M ρ T M  ≈ A y ,M ¯ Σ x A T y ,M + Σ b,M , ( 52a) Σ ρ,M , E  ρ M ρ T M  ≈ A y ,M ¯ Σ x A T y ,M + Σ b,M , (52b) where A y ,M = ( I N 2 + J )( I N ⊗ Y T M ) ∈ R N 2 × N P , (53a) A y ,M = ( I N 2 + J )( P ⊗ Y T M ) ∈ R N 2 × N P , (53b) and an expression for Σ b,M is derived in App endix D. 1) Uncon strained CRBs: The lowes t achievable variance by an unb iased estimator is gi ven by Σ y ,M , E n ( ˆ y M − y M )( ˆ y M − y M ) T o ≥ F † y ,M , ( 54a) Σ y ,M , E  ( ˆ y M − y M )( ˆ y M − y M ) T  ≥ F † y ,M , (54b) where the correspo n ding FI M s are given by F y ,M = A T Σ † ρ,M A , (55a) F y ,M = A T Σ † ρ,M A . (55b) It is worth noting that the Moore-Pen rose p seudoinverse is emp loyed since the FIM is rank- deficient, and consequ ently th e derived bo unds (54) are oracle-boun ds. 2) Constrained CRBs: When the FIM is rank-deficien t, a co nstrained CRB can be derived giv en differentiable a nd d eter- ministic con straints on the par ameters [28]. Let ¯ U , U be an orthono rmal basis for the nu ll space of th e con stra int matrices ¯ C , C , then the constraine d Cram ´ e r-Rao bound (CCRB) o n the M th or der kinematics are gi ven by Σ C y ,M , E n ( ˆ y M − y M )( ˆ y M − y M ) T o ≥ ¯ U  ¯ U T F y ,M ¯ U  − 1 ¯ U T , (56a) Σ C y ,M , E  ( ˆ y M − y M )( ˆ y M − y M ) T  ≥ U  U T F y ,M U  − 1 U T , (56b) where the FIMs are given b y (5 5). Re vised manuscript 12 10 13 16 20 25 32 40 50 63 79 100 10 -2 10 -1 10 0 P S f r a g r e p l a c e m e n t s No. of T wo-way comms. ( K ) RMSE Range [m] Range Rate: [m/s] Rate of Ran ge Rate: [m/s 2 ] RCRB -10 -8 -6 -4 -2 0 2 4 6 8 10 10 -2 10 -1 10 0 10 1 10 2 P S f r a g r e p l a c e m e n t s − 10 log 10 ( σ ) [dB meters] RMSE Range [m] Range Rate: [m/s] Rate of Ran ge Rate: [m/s 2 ] RCRB Fig. 2: Range parameters: V arying K : RMSEs (and RCRBs) of relative range parameters ( r , ˙ r , ¨ r ) for varying number of communications ( K ) between the N = 10 mobile nodes for σ = 0 . 1 meters. V arying σ : RMSEs (and RCRBs) of relative range parameters ( r , ˙ r , ¨ r ) for a network of N = 10 nodes exchan ging K = 10 timestamps, where the noise on the time markers ( σ ) is varied . Unlike our prev ious experimen ts [14], [18], we con sider acce leration in the current setup. X =  − 244 38 5 81 − 19 − 792 − 554 − 965 − 985 − 49 − 503 − 588 − 456 − 992 − 730 879 970 155 318 − 8 5 8 419  m (58a) Y 1 =  − 5 − 5 − 6 6 − 1 2 1 − 5 9 − 5 − 8 − 8 − 7 − 9 − 3 − 2 − 2 − 10 2 − 1  ms − 1 (58b) Y 2 =  − 0 . 17 − 0 . 17 0 . 22 − 0 . 07 0 . 21 − 0 . 15 0 . 55 − 0 . 72 − 0 . 49 − 0 . 34 0 . 42 0 . 42 0 . 98 0 . 73 0 . 48 0 . 08 − 0 . 43 − 0 . 14 0 . 5 6 0 . 91  ms − 2 (58c) D. Choice of weighting matrices ¯ W M , W M T o admit a BLUE solution, we use the in verse of the covariance matrice s Σ ρ,M , Σ ρ,M as weights to solve the regression problem s ( 37) and (44), i.e., ¯ W M , ˆ Σ † ρ,M =  ˆ A y ˆ ¯ Σ x ˆ A T y + ˆ Σ b,M  † , (57a) W M , ˆ Σ † ρ,M =  ˆ A y ˆ ¯ Σ x ˆ A T y + ˆ Σ b,M  † , (57b) where the estimates ˆ A y , ˆ A y are obtained by substituting ˆ Y M from LLS [(36) and (43)], in (53), ˆ ¯ Σ x is an estimate of (48) and ˆ Σ b,M is derived in Appen dix D from appropriate range parameter estimates. V I I . S I M U L A T I O N S In this section, we c onduc t experiments to validate the pr o posed d a ta m odel, and the solutions against their resp ectiv e de riv ed lower bound s. A network of N = 1 0 no des is considere d in P = 2 dimen sional space, with instantan eous po sition, velocity and acceler ation values arb itrarily chosen as in (5 8), such that the constraint (38) holds. All the n odes comm unicate with each other within a small time-interval of ∆ T = [ T ij,k , T j i,k ] = [ − 1 , 1] secon ds, wherein the transmit time markers are cho sen to be linearly spaced Without loss of generality , we are interested in the in stantaneou s kinematics of the nod es at time instant t 0 = T 0 = 0 . W e assume that all the pairwise co mmunica tions are ind ependen t o f each other, i.e ., Σ = σ 2 I ¯ N K . T h e metric used to ev alu a te the perf ormance of the range parameters is the root mean square error (RMSE), given b y RMSE ( z ) = N − 1 z v u u t N − 1 exp N exp X n =1 k ˆ z ( i ) − z k 2 , (59) where ˆ z ( i ) is th e estimate o f the unkn own vector z ∈ R N z × 1 related to the i th r un o f N exp = 500 Monte Carlo runs. T o ev alu a te th e estimates of the relative and absolute kinem atic matrices, we use z = vec ( U ) , where U is th e matrix u nder ev alu a tio n. T o q u alify these estimates, the square root of the C ram ´ er-Rao Boun d (RCRB) is plotted along with the respectiv e RMSE. It is worth n oting that the th eoretical lower bound s for th e range p arameters (47), and sub sequently the boun ds for relativ e position (49) and node kinematics (54, 56) are de p enden t on the cov ar iance of the noise on time markers i.e., Σ . Re vised manuscript 13 10 13 16 20 25 32 40 50 63 79 100 10 -3 10 -2 10 -1 P S f r a g r e p l a c e m e n t s (a) Relative position: X [m] No. of T wo-w ay comms. ( K ) RMSE MDS RCRB-Unconstrained 10 13 16 20 25 32 40 50 63 79 100 10 -2 10 -1 10 0 P S f r a g r e p l a c e m e n t s (b) Relative velocity: Y 1 [ms − 1 ] No. o f T wo-way comms. ( K ) LLS WLLS RCRB-Constrained RCRB-Unconstrained 10 13 16 20 25 32 40 50 63 79 100 10 -2 10 -1 10 0 10 1 P S f r a g r e p l a c e m e n t s (c) Relative acceleration Y 2 [ms − 2 ] No. of T wo-way comms. ( K ) LLS WLLS RCRB-Constrained RCRB-Unconstrained -10 -8 -6 -4 -2 0 2 4 6 8 10 10 -2 10 -1 10 0 10 1 P S f r a g r e p l a c e m e n t s (d) Relative position: X [m] − 10 log 10 ( σ ) [dB meters] RMSE MDS RCRB-Unconstrained -10 -8 -6 -4 -2 0 2 4 6 8 10 10 -2 10 0 10 2 P S f r a g r e p l a c e m e n t s (e) Relati ve velocity: Y 1 [ms − 1 ] − 10 log 10 ( σ ) [dB meters] LLS WLLS RCRB-Constrained RCRB-Unconstrained -10 -8 -6 -4 -2 0 2 4 6 8 10 10 -2 10 0 10 2 P S f r a g r e p l a c e m e n t s (f) Relative acceleration Y 2 [ms − 2 ] − 10 log 10 ( σ ) [dB meters] LLS WLLS RCRB-Constrained RCRB-Unconstrained Fig. 3: Rel ative Kinematics: V arying K : RMSEs (and R CRBs) of (a) Relati ve position ( X ), (b) Relativ e velocity ( Y 1 ) and (c) Relati ve acceleration ( Y 2 ) for v arying numb er of communications ( K ) between the N = 10 mobile nodes for σ = 0 . 1 meters. V arying σ : RMSEs (and R CRBs) of (d) Relative position ( X ), (e) R elativ e velocity ( Y 1 ) and (f) Relative acceleration ( Y 2 u ), for a network of N = 10 excha nging K = 10 timestamps, where the Noise on the time mark ers ( σ ) is varied. For all the proposed estimato r s in Section s VII A-C , we con duct two types o f e xperim ents. Firstly , for (a) varying num ber of pairwise commun ications K from 0 to 100 , with constant noise of σ = 0 . 1 m, and secondly for (b) v arying SNR from [ − 10 , 10 ] dB meter with a fixed K = 1 0 time-stamp exchan ges. The noise considered on the time-markers is typ ical of TWR based fixed localization experiments [29]. A. Range parameters W e employ th e dynamic ranging algorithm (2 6) for L = 3 , to estimate th e desired range co e fficients from the time-varying propag ation delays. In co mparison to our previous exper iments [14], [18], we addition ally co nsider acceler ation in the current simulation. Fig. 2 shows the RMSE and RCRB of the first 3 ran ge co efficients, for b o th varying K and varying SNR, wh ere we observe th at the RMSEs achieve the c o rrespon ding derived RCRBs asymptotically . Ob serve th at in th e Mon te carlo experim ents, we consider the noise on the time makers, whereas the lower bounds are deri ved on the data model with approximated noise (24). Hence, the RMSEs achieving th e correpo nding RCRB s validates our noise approximatio n d iscussed in Ap pendix A for the given exp e rimental setup. For the lin ear model ( 2 4), the p r oposed solution is the minimu m v ariance u nbiased estimator under Gaussian n oise assumption. In th is simulation , witho ut loss o f gener a lity , we assume that the o rder of ap proxim ation L is known. Alternatively , iterati ve solutions such as iMGLS [1 4] can b e employed to estimate L . For a d etailed discussion o n the effect o f L on the distance estimation , particularly for an asynchron ous network , see [14]. B. Relative kinematics The e stima ted r elativ e ran ge param eters yie ld the desired r elativ e kin ematics ma tr ices. Fig. 3 shows the RMSEs (and RCRBs) of all the relative kinematic estimates. Th e MDS-b ased r elativ e p osition estimates presented in Fig. 3 (a) and Fig. 3(d) , pe r form well against the derived or acle-bou nd, which was also o b served in [14]. In case of the relative velocity and acceleration, we assume the minimal con straint ¯ C 1 for analysis. No te that the unconstrained ora cle-boun ds are lo wer as compared to the CCRB, for a fixed SNR an d in creasing K. T he WLLS solution outperf orms the LLS solutions for both velocity an d acceleration estimation, and asympto tically achie ve th e deriv ed respectiv e CCR Bs. T o compare the p erform ance of the propo sed relati ve velocity estimator aga in st th e MDS-based relati ve velocity estimation (31), we perfo rm ano ther exp e r iment. Th e MDS-b ased alg orithm for relati ve velocity estimation assumes the n odes a re in linear motion. Hence, we set Y 2 = 0 P,N in (58) and re-implement the dyn amic ranging algorithm for L = 2 an d plot the standard deviation of the estimates in Fig. 4. Un der the constant velocity assumption, the CCRB is compa r able to the or acle-bou n d. The pro posed WLLS solution o utperfo rms the MDS-b ased estimator , especially for h igher SNR and lower nu mber of pair-wise Re vised manuscript 14 10 13 16 20 25 32 40 50 63 79 100 10 −3 10 −2 10 −1 10 0 10 1 P S f r a g r e p l a c e m e n t s No. of T wo-w ay comms. ( K ) RMSE LMDS WLLS RCRB-Constrained RCRB-Unconstrained −10 −8 −6 −4 −2 0 2 4 6 8 10 10 −2 10 −1 10 0 10 1 10 2 P S f r a g r e p l a c e m e n t s − 10 log 10 ( σ ) [dB meters] RMSE LMDS WLLS RCRB-Constrained RCRB-Unconstrained Fig. 4: Comparison of relativ e velocity esti mators: RMSE s (and RCR Bs) of relativ e range parameters Y 1 for varying number of communications ( K ) for σ = 0 . 1 meters (t op) and varying σ (bottom) between the N = 10 mo bile nod es. 10 13 16 20 25 32 40 50 63 79 100 10 -3 10 -2 10 -1 10 0 P S f r a g r e p l a c e m e n t s (a) Absolute v elocity: Y 1 [ms − 1 ] No. of T wo-way comms. ( K ) GLLS WGLLS RCRB:Anchored RCRB:Anchor-Free 10 13 16 20 25 32 40 50 63 79 100 10 -2 10 -1 10 0 10 1 P S f r a g r e p l a c e m e n t s (b) Absolute acc eleration Y 2 [ms − 2 ] No. of T wo-w ay comms. ( K ) GLLS WGLLS RCRB:Anchored RCRB::Anchor-Free -10 -8 -6 -4 -2 0 2 4 6 8 10 10 -2 10 0 10 2 P S f r a g r e p l a c e m e n t s (c) Absolute v elocity: Y 1 [ms − 1 ] − 10 log 10 ( σ ) [dB meters] GLLS WGLLS RCRB:Anchored RCRB:Anchor-Free -10 -8 -6 -4 -2 0 2 4 6 8 10 10 -2 10 0 10 2 P S f r a g r e p l a c e m e n t s (d) Absolute acc eleration Y 2 [ms − 2 ] − 10 log 10 ( σ ) [dB meters] GLLS WGLLS RCRB:Anchored RCRB:Anchor-Free Fig. 5: Absolute Kinematics: V arying K : RMSEs (and RCR B s) of (a) Absolute velocity ( Y 1 ) and (b) Absolute acceleration ( Y 2 ) for v arying number of co mmunications ( K ) between the N = 10 mobile nodes for σ = 0 . 1 meters. V arying σ : RMSEs (and R CRBs) of (c) Absolute velocity ( Y 1 ) and (d) Absolute acceleration ( Y 2 ), for a network of N = 10 nodes exchanging K = 10 timestamps, where the noise on the t i me markers ( σ ) is v aried. commun ications. This is per haps n ot surprising, since the MDS-based estimator relies on all the R , ˙ R , ¨ R where the noise variance on these regression coefficients typic a lly incr ease with the rang e-order fo r a T ay lor basis (see Fig . 2). In comparison, the WLLS solution is dep endent only on range R and range rates ˙ R . C. Absolute kinematics Fig. 5 sh ows the RMSEs an d the corr espondin g RCRB s of the absolute velocity Y 1 and acceleration Y 2 . W e assume constraint (45) to solve the proposed GLLS (43) and WGLLS (44) algorithms. The propo sed estimato rs are seen to conv erge asymptotically to the d eriv ed CCRBs, while the CCR B itself is an order higher th an the theoretical or acle-bou nd. Th e perfor mance of the absolute kin ematics is very similar to tha t of th e relativ e kinematics (see Fig. 3), which is due to the fact that the FIMs in both scenarios are dominate d by the singu lar v alues of the relative po sition matrix. D. Relative and absolu te time-varying positions The estimatio n of the no de kinema tics enable us to re c o nstruct the time-varying relati ve positions S ( t ) and time-varying absolute po sitions S ( t ) , f rom (3 9) and (46) respectively . W e cond uct experiments to study the effect of the prop osed estimators on the tim e-varying p ositions. The RMSE p lo t fo r the ab so lute and relative time-varying p o sitions ar ound the region of interest at t 0 = 0 are shown in Fig. 6, where the n umber of com munication s a re varied as K = [50 , 100 , 50 0] with a Gaussian noise Re vised manuscript 15 -1 -0.5 0 0.5 1 10 -2 10 -1 10 0 10 1 K=50 K=100 K=500 σ =1 meter P S f r a g r e p l a c e m e n t s Relativ e T ime-varying Position: S ( t ) T ime ( t ) RMSE K = 5 0 K = 1 0 0 K = 5 0 0 σ = 1 m e t e r -1 -0.5 0 0.5 1 10 -2 10 -1 10 0 10 1 K=50 K=100 K=500 σ =1 meter P S f r a g r e p l a c e m e n t s Absolute Time-v arying Position: S ( t ) T ime ( t ) RMSE K = 5 0 K = 1 0 0 K = 5 0 0 σ = 1 m e t e r Fig. 6: P osition ov er time: RMSEs of r el at ive position S ( t ) and ab solute p osition S ( t ) over time for K = [50 , 100 , 500] communications between a clu ster of N = 10 mob ile nod es, with σ = 1 meter . 2 3 4 5 6 10 -2 10 -1 P S f r a g r e p l a c e m e n t s (a) Relative velocity: Y 1 [ms − 1 ] No. of relati vely i mmobile nodes GLLS WGLLS RCRB:Constrained 2 3 4 5 6 10 -1 10 0 P S f r a g r e p l a c e m e n t s (b) Relati ve acceleration Y 2 [ms − 2 ] No. of relati vely i mmobile nodes GLLS WGLLS RCRB:Constrained 2 3 4 5 6 10 -2 10 -1 P S f r a g r e p l a c e m e n t s (c) Absolute v elocity: Y 1 [ms − 1 ] No. o f known node kine matics GLLS WGLLS RCRB:Anchored R C R B : A n c h o r - F r e e 2 3 4 5 6 10 -1 10 0 P S f r a g r e p l a c e m e n t s (d) Absolute acc eleration Y 2 [ms − 2 ] No. of k no wn node kinematics GLLS WGLLS RCRB:Anchored R C R B : A n c h o r - F r e e Fig. 7 : E ffec t of increasing constraints: Relativ e ki nematics: RMSEs (and RCRB) of (a) Relativ e velocity ( Y 1 ) and (b) Relativ e acceleration ( Y 2 ) for varying number o f relati vely immobile nodes. Absolute kinematics: RMSEs (and RCRBs) of (c) Absolute velocity ( Y 1 ) and (d) Absolute acceleration ( Y 2 ) for v arying numb er of kno wn node kinematics on the distance of σ = 1 meter . For K = 500 , the RMSE estimate of bo th the relati ve and absolute position around t 0 shows an imp rovement by an order mag nitude in compariso n to the noise on the d istance me asurement, for the g iv en experime n tal setup. This g ain is primarily contributed during dyna mic ranging, where K da ta poin ts are averaged using the T ay lor basis which yields a factor √ K improvement on the estimate of the ran g e parameters. Secondly , the performan ce deterio r ates as we move away f rom t 0 , which is a ty p ical characteristic of the T ay lor approximatio n. However , if Doppler measurements are av ailab le for radial velocities and other higher ord er d e riv atives, then the standar d deviation of the estimators can be furth e r reduced . E. Choice of constraints In the pre vious sections, we e valuated th e pro posed algorithms under m in imal constrain ts. Now , we perform experiments to understand the ef fect of incorpo rating additio n al co nstraints (or references) on the perfo r mance o f the proposed estima to rs. These ad ditional co nstraints im plicitly r educe the parameter sub space, and con sequently affect the overall RMSE of the p ropo sed estimators. In order to understand this variation, we set N z = 1 in o ur performance m etric (59) fo r th e following simu lations. T o estimate the r elativ e kinematics in a 2 dimension al scenario, a unique solution is feasible if at least 2 nodes are relatively immobile (see Appendix B). I f m o re n odes a r e immobile, the n the constrain ts in (38) can be readily extended to incorp orate this supplemen tary information. Similarly , in case of absolute velocity and acceler ation estimation, a m inimum of at lea st 2 nod e kinematics must b e known. Therefore , in the following experim ents we vary the number o f known kinematics ( or immobile Re vised manuscript 16 nodes) from 2 to 6 , for a fixed numb er of two way com munication s K = 100 with σ = 0 . 1 meters. Fig. 7 shows the results of the GLLS an d WGLLS alg o rithms for estimating the absolute and r elativ e k in ematics, a lo ng with the r espective C CRBs. Not surprisin g ly , we observe an imp rovement in the performance of th e algo rithms with the additional co nstraints. In add ition, unlike the GLLS estimator , the WGLLS estimator asymptotically achieves the respec tive CCRBs. V I I I . C O N C L U S I O N S Understand ing the relative kinematics o f an anchorless network of mobile nod es is paramount for r eference- free localization technolog ies of the futu re. W e presented a novel data model wh ich relates the time-varying distance m e asurements to the M th order relative kin e m atics for an an chorless network of m obile no des. The der ived data model takes the form of a L y apunov-like equation, which und er certain constraints, can be recur sively solved fo r estimating the r elativ e velocity , acc e leration an d h igher order deriv a tives. Closed f orm constrained estimato rs, such as th e LS and WLS are proposed, which are also the BLUE for the gi ven data model. Cram ´ er-Rao lower bo u nds are derived for th e ne w data mode l and the p erform a nce of the pr oposed algorithm s is v alidated using simu la tio ns. Althou g h ou r focus is on r elativ e localization , th e proposed model an d solution s can be broadly applied to understand feature variations in Euclidean space, with applica tions in general exp lo ratory data a n alysis. In our future work, we are keen in addressing tw o research challenges. Firstly , our focus in this article has been o n finding unique solutions to time-d e riv atives of the r e lati ve position matrix. T o this en d, unbiased constrained estimators are p roposed to so lve the under-determined L y apunov-like equation . However , more generally , regularize d algo rithms can be employed, such as Ridge regression [ 30], subset selection [ 3 1] or Lasso [32], without the need fo r equality constrain ts on the cost function . The estimates o f such unconstrain ed alg orithms can be cor r obora te d against the unconstrained Cram ´ er-Rao bou n d derived in this article. Furth ermore, the algo r ithms a r e inh erently ce ntralized in nature, which could be distrib uted for resource constrained implementation. Finally , the p r oposed fram ew ork is particularly he lp ful f or cold -start scena r ios when there is no apriori in formatio n on th e position or hig her ord er kine m atics. In p ractice, given the cold- start solution on relative velocity and higher ord er kinem a tics, a state-space model readily e merges for dyn a mic tracking of the relative positions over time, which can be elegantly solved u sing adaptive filters. A P P E N D I X A A P P ROX I M AT E N O I S E M O D E L T o estima te the ra n ge parameter s f rom time-varying propaga tion delays, we pr esented the d ynamic rang ing model in (24), with additive Gaussian noise i.e . , V θ = τ + η , (60) where V is the V and ermond e-like matrix, θ contain s th e unknown range coefficients, τ contains all th e propa g ation dela y s and η is the noise vector plaug ing the p ropagatio n d elays. In practise, the noise is o n the time markers and subsequently on the V andermo nde matrix. Howe ver , u nder cer tain no minal assumptio n s, the above model is v alid, which we discuss in this section. W e begin with the noiseless pairwise time-varying dynamic ranging model, which we recollect from (22) as below r ij + ˙ r ij ∆ T k + ¨ r ij ∆ T 2 k + . . . = | T ij,k − T j i,k | = τ ij,k , (61) where we in troduce ∆ T k = ( T ij,k − T 0 ) for no tational simplicity . In r eality , there is noise plagu ing the time m arkers and h ence we have, r ij + ˙ r ij (∆ T k + η i,k ) + ¨ r ij (∆ T k + η i,k ) 2 + . . . = τ ij,k + η ij,k , (62) where { η i,k , η j,k } are the no ise ter ms on the time mar kers at n ode i and node j respectively , and η ij,k = η i,k − η j,k is the pairwise noise error of the no d e pair ( i, j ) . Ex pandin g the po lynomial and rearran g ing the terms, we ha ve r ij + ˙ r ij ∆ T k + ¨ r ij ∆ T 2 k + . . . + ¯ η i,k = τ ij,k + η ij,k . (63 ) Here ¯ η i,k is the cumulative noise err or from the T aylor appro ximation, which is e xpressed as ¯ η i,k = η i,k  ˙ r ij,k + 2¨ r ij,k ∆ T k + . . .  + η 2 i,k  ¨ r ij,k + . . .  + . . . ≈ 0 , (64) and appro ximated to 0 . This appro x imation is v alid under two assumption s. Firstly , we assume that the time stamps are measured with high SNR, i.e., we consider standard de viations of ≤ 10 − 7 seconds on the time stamps, which is necessary to achieve meter le vel accuracies is co nventional two-way ranging b ased localization solutions [ 33], [34]. As a c o nsequen ce, we ignore the second order noise term η 2 i,k , and o ther high er order no ise terms in (64). Secondly , observe from definitio n (21) that the coef ficients { ˙ r , ¨ r, . . . } are scaled by c − 1 , where c = 3 × 10 8 m/s for free space. Therefor e, the T aylor coefficients are significantly small and subsequently , the term ( ˙ r ij,k + 2¨ r ij,k ∆ T k + . . . ) is negligible fo r a measurement period of upto a f ew seconds. This is a p ragmatic assump tion, since we are on ly interested in the instantan eous relativ e kinematics of the no des around a small time interval. In summar y , for small measure ment per iods in hig h SNR scen a rios, the n oise pa rameter ¯ η i,k ≈ 0 , and under these assumptio n s (24) holds. Re vised manuscript 17 A P P E N D I X B U N D E R D E T E R M I N E D L Y A P U N OV - L I K E E Q U A T I O N Theor em 1 (Underdetermined Lyapunov-like eq uation) : Given X ∈ R P × N and B ∈ R N × N for N > P , the L yapun ov-like equation X T Y + Y T X = B , (65) is rank-de ficient by at least ¯ P =  P 2  . Pr oof: Let th e singu lar value decomposition of X be X = U x  Λ x 0  V T x , (66) where Λ x ∈ R P × P is a diago nal matrix con taining the singular values and U x ∈ R P × P , and V x ∈ R N × N are the correspo nding singular vectors. Then, (65) is  Λ x 0  T ˜ Y + ˜ Y T  Λ x 0  = ˜ B , (67) where ˜ B =  ˜ B 11 ˜ B 12 ˜ B T 12 ˜ B 22  = V T x B x V x , (68) ˜ Y =  ˜ Y 1 ˜ Y 2  = U T x YV x , (69 ) where ˜ Y 1 ∈ R P × P , ˜ Y 2 ∈ R P × N − P and ˜ B 22 = 0 for the equation to be consistent. A solution to the system ( 6 5) is obtained by solving for ˜ Y the set of equation s, Λ x ˜ Y 1 + ˜ Y T 1 Λ x = ˜ B 11 , (70) Λ x ˜ Y 2 = ˜ B 12 . (71) An estimate for ˜ Y 2 is straightforward and is gi ven by ˆ ˜ Y 2 = Λ − 1 x ˜ B 12 . Let ˜ Λ x , ˜ Y 1 and ˜ B 11 be partitioned into  σ 1 0 0 Λ x, 1  ,  y 11 ˜ y 12 ˜ y 21 ˜ Y 1 , 1  ,  ˜ b 11 ˜ b 12 ˜ b T 12 ˜ B 11 , 1  , (72) then (70) is equivalent to solving y 11 = ˜ b 11 / 2 σ 1 , (7 3 ) σ 1 ˜ y 12 + ˜ y T 21 Λ x, 1 = ˜ b 12 , (74) Λ x, 1 ˜ Y 1 , 1 + ˜ Y T 1 , 1 Λ x, 1 = ˜ B 11 , 1 . ( 7 5) Note th at the solutio n to y 11 in (7 3) is straigh tforward, howe ver the so lution to off-diago nal term s ˜ y 12 , ˜ y 21 is u nderd e termined. Furthermo re, since (75 ) is in for m similar to the (70) , ˜ Y 1 , 1 can b e estimated recu rsiv ely [35]. Th us, the diagonal ter m s of the P dimensional matrix ˜ Y 1 , 1 can be estimated , ho wever to resolve th e ambiguity o f the off-diago nal ter ms atleast ¯ P =  P 2  constraints are required . A P P E N D I X C K A RU S H - K U H N - T U C K E R ( K K T ) E Q U A T I O N S A solution to minimize the equ ality constrained function of the form min y k Ay − b   2 s.t. Cy = d , (76) is obtained by solving the Karu sh-Kuhn-T ucker (KKT ) system of equations,  ˆ y ˆ λ  =  2 A T A C T C 0 N 2 ,N 2  − 1  2 A T b d  , (77 ) where ˆ y is an estimate of the unknown pa rameter y and ˆ λ collects the correspondin g La grange multipliers. The problem has a feasible solution provide d  A C  is full column rank [36]. Re vised manuscript 18 A P P E N D I X D E X P R E S S I O N F O R Σ b,M W e presen t an explicit expression for the covariance matrix Σ b,M , which is obtained by ignor in g higher order noise terms i.e., for sufficiently large SNR. For M = 1 , i.e., rela tive velocity , we ha ve Σ b, 1 ≈ ˜ P  Ψ r ¯ Σ ˙ r Ψ r + Ψ ˙ r ¯ Σ r Ψ ˙ r  ˜ P , (78) and for M = 2 , i.e., relative acce leration, we have Σ b, 2 ≈ ˜ P  Ψ r ¯ Σ ¨ r Ψ r + Ψ ¨ r ¯ Σ r Ψ ¨ r + 4 Ψ ˙ r ¯ Σ ˙ r Ψ ˙ r  ˜ P + 4 Ψ y ¯ Σ ˙ x Ψ y , (79) where we ˜ P , P ⊗ P , Ψ r , d iag  vec  R )  , Ψ ˙ r , d iag  vec  ˙ R )  and Ψ ¨ r , d iag  vec  ¨ R )  . The matrix Ψ y = A y , 1 for absolute kinematics an d Ψ y = A y , 1 for relati ve k inematics. Observe that the diagonal elements of th e range p arameters R , ˙ R , ¨ R , . . . contain zeros an d conseq uentially the matrices Ψ r , Ψ ˙ r , Ψ ¨ r , . . . are singular . Hence the covariance matrix Σ b,M is in general ran k d eficient. F urther more, A y in (53b) is rank deficient by definition and subsequ ently Σ ρ (52) is ill-co nditioned and theref o re, we u se the Moor e-Penrose pseu do-inverse in (55) and (57). An expression for higher order M > 2 can be similarly derived. R E F E R E N C E S [1] R. T . Rajan, Relative Space- T ime Kinematics of an Anchorless Network . Doctoral dissertation , Gilde print, The Netherlands, 2016. [Online]. A vaila ble: http:/ /doi.org/ 10.4233/uuid:0bcfc55b-be81-4326-855c-3a97ba126521 [2] M. Koehl er , T . Rabino witz, J. Hirdes, M. St ones, G. I. Carpente r , B. E. Fries, J . N. Morris, and R. N. Jones, “Measuring depression in nursing home resident s with the mds and gds: an observ ational psyc hometric study , ” BMC geriatrics , vol. 5, no. 1, p. 1, 2005. [3] C.-C. Ho, K . F . MacDorman, and Z. D. Pramono, “Human emotion and the uncanny vall ey: a glm, mds, and isomap analysis of robot video ratings, ” in Proc eedings of the 3rd ACM/IEEE internat ional confer ence on Human r obot intera ction . A CM, 2008, pp. 169–176. [4] B. Dil, S. Dulman, and P . Havinga, “d, ” in W ire less Sensor Networks . Spr inger , 2006, pp. 164–179. [5] D. K. Agra fiotis, D. N. Rassokhi n, and V . S. Lobano v , “Multidime nsional s caling and visuali zation of l arge molec ular similarity table s, ” Journal of Computati onal Chemistry , vol. 22, no. 5, pp. 488–500, 2001. [6] F . Rehm, F . Klawo nn, and R. Kruse, “Mds polar: a new approach for dimension reduction to visualize high dimensional data, ” in Advances in Intellig ent Data Analysis V I . Spri nger , 2 005, pp. 316–327. [7] O. C. Jenkins and M. J. Matari ´ c, “ A spatio-temporal extension to isomap nonlinear dimension reduction, ” in Pro ceedings of the twenty-first internatio nal confer ence on Machine learning . A CM, 2004, p. 56. [8] I. Borg and P . J. F . Groenen, Modern Multidimen sional Sc aling: Theory and Application s (Springer Series in Statisti cs) , 2nd ed. Spr inger , 8 2005. [9] R. T . Rajan, A.-J. Boonstra, M. Bentum, M. Klein-W olt, F . Belien, M. Arts, N. Saks, and A. -J. V een, “Space -based aperture array for ul tra-long wa velen gth radio astronomy , ” Experimental Astro nomy , vol. 41, no. 1, pp. 271–306, 2 2016. [Online]. A va ilable : ht tp://dx .doi.org/10.1007 /s10686-015-9486-6 [10] V . Chandrasek har , W . K. Seah, Y . S. Choo, and H. V . Ee, “Localiza tion in underwat er se nsor network s: surv ey and challeng es, ” in Proc eedings of the 1st ACM international workshop on Underwater networks . A CM, 2006, pp. 33–40. [11] Z. Y ang, C. W u, and Y . Liu, “Loca ting in finge rprint space: wirel ess indoor locali zation with little human interv ention, ” in Proc eedings of t he 18th annual internatio nal confere nce on Mobile computing and networking . AC M, 2012, pp. 269–280. [12] F . W . Y oung, Multidimensio nal scali ng: History , theory , a nd application s . Psychology Press, 2013. [13] D. Halliday , R. Resnick, and J. W alk er , Fundame ntals of physics exte nded . John W ile y & Sons, 2010, vol. 1. [14] R. T . Rajan and A.-J. v an der V een, “Joint ranging and synchro nizatio n for an anchorless network of mobile nodes, ” IE EE T ransactions on Signal Pr ocessing , , vol. 63, no. 8, pp. 1925–1940, 4 2015. [15] S. M. Kay , Fundamentals of statistical signal pr ocessing: estimation theory . Upper Saddle Ri ver , NJ, USA: Prentice -Hall, Inc., 1993. [16] H.-W . W ei, R. Peng, Q. W an, Z.-X. Chen, and S.-F . Y e, “Multidimensio nal Scali ng Ana lysis for Passiv e Mo ving T arget Localizati on W ith TDOA and FDOA Measurements, ” IE EE T ransactions on Signal P r ocessing , vol. 58, no. 3, pp. 1677 –1688, 3 2010. [17] H. Jamali-Rad a nd G. Leus, “Dynamic mult idimensional sca ling for low-co mplexit y mobile network tr acking, ” Signal Proce ssing, IEEE T ransactions on , vol. 60, no. 8, pp. 4485–4491, 2012. [18] R. T . Raj an, G. Leus, and A.-J. v an der V een, “Joint relati ve positio n and velocit y estimati on for an anchorl ess netwo rk of mobile nodes, ” Signal Pr ocessing , vol. 115, no. 0, pp. 66 – 78, 10 2015. [Online]. A vai lable: http://www .sciencedirec t.com/science/article/pii/S0165168415000894 [19] R. T . Rajan and A.-J. va n der V een, “Joint motion estimation and clock synchronizat ion for a wireless network of mobile nodes, ” in IEEE International Confer ence on A coustic s, Speec h and Signal Processi ng (ICASSP) , March 2012, pp. 2845–2848. [20] J. H. Hodges, “Some m atrix equation s over a finite field, ” Annali di Matematica Pura ed A pplicata , vol. 44, no. 1, pp. 245–250, 1957. [21] C.-Y . Chiang, E . K. -W . Chu, and W .-W . Lin, “On the ⋆ − sylvest er equation ax ± x ⋆ b ⋆ = c , ” Applied Mathematics and Computation , vol. 218, no. 17, pp. 8393–8407, 2012. [22] R. Horn and C. R. Johnson, T opics in matrix analysis . Cambridge Uni versit y Press, 1991. [23] R. H. Bartels and G.W .Ste wart, “Solutio n of the matrix equation ax + xb = c , ” Communication s of the ACM , vol. 15, no. 9, pp. 820–826, 1972. [24] G. H. Golub, S. Nash, and C. V an Loan, “ A hessenberg-schur method for the problem ax+ xb= c, ” IEEE T ransactions on Automatic Contr ol , vol. 24, no. 6, pp. 909–913, 1979. [25] R. Bhatia and P . Rosenthal, “How and why to solv e the operato r equat ion AX − X B = Y , ” B ulleti n of the London Math ematical Soci ety , vol . 29, no. 1, pp. 1–21, 1997. [26] H. Braden, “The equati ons A T X ± X T A = B , ” SIAM Journal on Matrix A nalysis and Applicati ons , vol. 20, no. 2, pp. 295–302, 1998. [27] T . Penzl, “Numerical solution of generali zed lyapunov equations, ” Adva nces in Computati onal Mathematics , vol. 8, no. 1-2, pp. 33–48, 1998. [28] P . Stoica and B. C. Ng, “On the Cramer-R ao Bound under parametric constraints, ” IE EE Signal P r ocessing Letters , vol. 5, no. 7, pp. 177–179, 1998. [29] N. Patw ari, A. Hero, M. Perkins, N. C orreal, and R. O’Dea, “R elati ve loca tion e stimation in wi reless sensor networks, ” IEEE Tr ansactions on Signal Pr ocessing , vol. 51, no. 8, pp. 2137–2148, 2003. [30] G. H. Gol ub, P . C. Hansen, and D. P . O’L eary , “ Tikho nov regul arizati on and total least squares, ” SIAM Jou rnal on Matrix Analysis and A pplicat ions , vol. 21, no. 1, pp. 185–194, 1999. [31] C. L. Lawson and R. J. Hanson, Solving least squares pr oblems . SIAM, 1974, vol. 161. Re vised manuscript 19 [32] R. Tibsh irani, “Regression shrinkage and selectio n via the lasso, ” J ournal of the Royal Stati stical Societ y . Series B (Methodolog ical) , pp. 267–288, 1996. [33] N. Patw ari, J. Ash, S. Kyperoun tas, I. Hero, A. O . , R. Moses, and N. Correal, “Locating the nodes: Cooperati ve localiz ation in wireless sensor networks, ” IEEE Signal Proc essing Magazine , vol. 22, no. 4, pp. 54 – 69, 7 2005. [34] A. Simonetto and G. Leus, “Distribute d maximum likeli hood sensor network localizat ion. ” IEE E T rans. Signal Pro cessing , vol. 62, no. 6, pp. 1424–1437, 2014. [35] K.-w . E . Chu, “Sy mmetric solu tions of linear mat rix equ ations by m atrix decompositi ons, ” Linear Alg ebra and its Appli cations , v ol. 119, pp. 35–50, 1989. [36] S. Boyd and L. V andenberghe , Conv ex Optimization . Cambridge Univ ersity Press, Mar . 2004. aX -1 -0.5 0 0.5 1 aY 10 -2 10 -1 10 0 10 1 aT W G LLS : K = 50 W G LLS : K = 100 W G LLS : K = 100 σ = 1 m e t e r rX -1 -0.5 0 0.5 1 rY 10 -2 10 -1 10 0 10 1 rT WLLS: K=50 WLLS: K=100 WLLS: K=100 σ =1 meter −20 −10 −0 10 20 10 −4 10 −2 10 0 10 2 10 4 zSX zSY zST zSLa zSLb zSLc zSLd −20 −10 −0 10 20 10 −4 10 −2 10 0 10 2 ySX ySY yST ySLa ySLb ySLc ySLd