Relative kinematics of an anchorless network

Delft University of Technology

Relative kinematics of an anchorless network

Rajan, Raj Thilak; Leus, Geert; van der Veen, Alle Jan

DOI

10.1016/j.sigpro.2018.11.005

Publication date

2019

Document Version

Final published version

Published in

Signal Processing

Citation (APA)

Rajan, R. T., Leus, G., & van der Veen, A. J. (2019). Relative kinematics of an anchorless network. Signal

Processing, 157, 266-279. https://doi.org/10.1016/j.sigpro.2018.11.005

Important note

To cite this publication, please use the final published version (if applicable).

Please check the document version above.

Copyright

Other than for strictly personal use, it is not permitted to download, forward or distribute the text or part of it, without the consent of the author(s) and/or copyright holder(s), unless the work is under an open content license such as Creative Commons. Takedown policy

Please contact us and provide details if you believe this document breaches copyrights. We will remove access to the work immediately and investigate your claim.

This work is downloaded from Delft University of Technology.

Green Open Access added to TU Delft Institutional Repository

‘You share, we take care!’ – Taverne project

ContentslistsavailableatScienceDirect

Signal

Processing

journalhomepage:www.elsevier.com/locate/sigpro

Relative

kinematics

of

an

anchorless

network

R

Raj

Thilak

Rajan

a ,∗

_,

_Geert

_Leus

a

_,

_Alle-Jan

_van

_der

_Veen

a

TU Delft, Delft, The Netherlands

a

r

t

i

c

l

e

i

n

f

o

Article history: Received 27 July 2016 Revised 8 November 2018 Accepted 10 November 2018 Available online 23 November 2018

Keywords: Lyapunov-like equation Relative velocity Relative acceleration Multidimensional scaling Time-varying distance

a

b

s

t

r

a

c

t

The estimationofthecoordinatesofnodes theirproximity(or distance) measurements,isaprincipal challengeinnumerousfields.Conventionally,whenlocalizingastaticnetworkofimmobilenodes, non-lineardimensionalityreductiontechniquesareappliedonthemeasureddistancestoobtaintherelative coordinatesuptoarotationand translation.Inthisarticle,weconsider ananchorlessnetworkof mo-bile nodes,wherethedistance measurementsbetweenthemobile nodes aretime-varying.Insuchan anchorlessframework,wheretheabsoluteknowledgeofanynodeposition,motionorreferenceframeis absent,weaimtoestimatetherelativepositionsusingthemeasuredtime-varyingdistances.Tothisend, wederiveadatamodelwhichrelatesthetime-varyingdistancestothetime-varyingrelativepositions ofananchorlessnetwork.Giventhisdatamodel,weestimatetherelative(position,velocity)andhigher orderderivatives,whicharecollectivelytermedastherelativekinematicsoftheanchorlessnetwork.The derived datamodelis inherently ill-posed,howeverunder certainimmobility constraints,wepropose closed-formsolutionstorecursivelyestimatetherelativekinematics.Forthesakeofcompleteness, we alsoestimatetheabsolutenodekinematics,givenreferenceanchors.Theoreticalboundsarederived,and simulationsareconductedtobenchmarktheperformanceofproposedsolutions.

© 2018PublishedbyElsevierB.V.

1. Introduction

The estimation of the relative coordinates of N points (or nodes)in a P-dimensional Euclideanspaceusing proximity mea-surements(orpairwisedistances)isafundamentalproblem span-ning a broad range of applications. These applications include, butarenot limitedto,psychometricanalysis[2] ,perceptual map-ping [3] , range-based anchorless localization [4] , combinatorial-chemistry[5] ,polar-based navigation[6] ,sensorarray calibration

[7]and in generalexploratorydata analysis[8] .Inanchorless lo-calizationscenariosforinstance,nodesheavilyrelyonco-operative estimationof relativecoordinates. Such anchorless networks nat-urally arise when nodes are inaccessible or only intermittently monitored,as is the case in space-basedsatellite arrays [9] , un-derwater networks[10] or indoor wireless sensor networks[11] . Insuch reference-freescenarios, the proximityinformation, often measuredaspairwisedistancesbetweenthenodes,formakey in-putinestimatingtherelative coordinatesofnodes.Theserelative coordinatesare typicallyestimated usingnon-linear dimensional-ityreductionalgorithms(suchasmultidimensionalscaling(MDS)), whichhave beenstudied rigorouslyover the pastdecades [8,12] .

R _{A part of this work is published in the doctoral dissertation [1] .} ∗ _{Corresponding author.}

E-mail addresses: rtrajan@ieee.org (R.T. Rajan), g.j.t.leus@tudelft.nl (G. Leus),

a.j.vanderveen@tudelft.nl (A.-J. van der Veen).

However,considerablylessattentionhasbeendirectedtowards an-chorlessmobilescenarios.

Ourprimary focusin thisarticle is onan anchorless network of mobile nodes, where we use the term anchorless to indicate noabsoluteknowledgeofthenodepositions,motionorreference frame. Furthermore, since the nodes are mobile, both the node positions and the pairwise distance measurements between the nodesaretime-varyinginnature.Ourmotiveistorelatethe time-varyingpairwisedistancemeasurementstotime-derivativesofthe nodecoordinates.Forananchorlessnetwork,theseincludethe rel-ative position, relative velocity, relative acceleration and higher-orderderivativeswhichwe cumulativelyrefer toasrelative kine-matics in thisarticle. It is worth noting that the universally ac-cepteddefinitionofrelativekinematicsinherentlyreliesonthe in-formationin theabsolutereferenceframe.Forexample,the non-relativisticrelativevelocitybetweentwo objectsisrightlydefined asthedifferencebetweentheirrespectiveabsolutevelocityvectors

[13] .Inananchorlessframeworkhowever,anaturalquestionarises on whetherthe relative kinematics can be estimated, given only time-varyingdistancemeasurements.Ergo,wewishtounderstand therelationshipbetweenthetime-varyingdistancemeasurements andtherelativekinematicsofmobilenodes,whichistheprime fo-cusofthisarticle.Theestimatedrelativekinematicscanbereadily usedtofindthetime-varyingrelativepositionsofthenodes.

https://doi.org/10.1016/j.sigpro.2018.11.005

1.1. Previouswork

A key challenge in our pursuit is that both the time-varying distanceandthe time-varyingrelative positionsare non-linearin nature.Inparticular,theEuclideandistancebetweenapairof mo-bile nodesisalmostalways anon-linear functionoftime,even if the nodes are in linear independent motion[14] . Therefore,it is perhapsnotsurprisingthattraditionalmethodstosolvesuch chal-lenges have been to employ state-space based approaches, with the assistance of known anchors [15] . The initial position of the nodesis estimatedusingMDS-likealgorithms, whichusethe Eu-clideandistancematrix(EDM)atasingletimeinstanttoestimate therelativenodepositions.Giventhisinitialestimate,therelative positionsaretrackedoveraperiodoftimewithDoppler measure-ments andknownanchors[16] ,orviasubspacetrackingmethods

[17] .Unfortunately,Dopplermeasurementsandanchorinformation are notalwaysavailable. Secondly,subspacetrackingisapplicable only for smallperturbations inmotion andtherefore offers little insightonthekinematicsofthemotionitself.

Inour previousstudy,we proposed atwo-step solutionto es-timaterelative velocitiesofthe nodesfromtime-varyingdistance measurements [18] . Firstly, the derivatives of the time-varying distances were estimated by solving a Vandermonde-like system of linear equations. The estimated regression coefficients (called rangeparameters)jointlyyield therelativevelocities andthe rel-ativepositions,usingMDS-likealgorithms.However,theproposed solutionisvalidonlyforlinearmotion,whichisnotalways prac-tical.Furthermore,thepreviouslyproposedMDS-basedrelative ve-locityestimatorheavily reliesonthesecond-ordertime-derivative of distance, and under Gaussian noise assumptions, it performs worse than the relative position estimator. Thus, designing more optimalestimators forthe relativevelocity isone ofthekey mo-tivations for the pursuit of a generalized framework presented in this article.Moreover, understanding the higher order relative kinematicsofmotioninEuclideanspaceviatime-varyingdistance measurements iscrucialfornext-generationlocalization technolo-gies.

1.2. Contributions

Ourkeycontributionsaresummarizedasfollows.

1. We derive a generalized relative kinematics model for a net-work of mobile nodes, relating the derivatives of the time-varyingdistancemeasurementsbetweentherespectivepairsof mobilenodestotheirindividualrelativekinematics.Unlikeour previouswork[18] ,wherewelimitedourstudytorelative po-sitionandrelativevelocity,theproposedmodelinthisarticleis moregenerallyapplicableforrelativeposition,velocity, acceler-ationandhigher-orderkinematics.

2. Weproposealgorithmstoestimatetherelativekinematics, un-derrelativeimmobilityconditionsofafewnodes.Theproposed algorithms are novel for relative acceleration estimation, and simulations reveal that the proposed relative velocity estima-torsoutperformourpreviousMDS-likealgorithm[18] . 3. Forthesakeofcompletion,inthepresenceofanchor

informa-tion, we show that the absolutekinematics ofthe nodes can alsobeestimatedusingthederivedmodel.

4. Giventherelative(andabsolute)kinematicestimatesuptothe desiredorder,weshow thatthetime-varyingrelative(and ab-solute positions) of the nodes can be subsequently obtained. Simulations show that the proposed kinematics-based time-varying position estimation, offers significant improvement in positionaccuracyaroundthetime-periodofinterest.

1.3.Overview

WepresentthedatamodelinSection 2 ,whichrelatesthe time-varying distances to the kinematics of the mobile nodes. More concretely, this relationship is established via the derivatives of thetime-varyingdistance (calledrange parameters),which is es-timatedinSection 3 usingdynamicranging.InSection 4 weshow that the relationship between the rangeparameters andthe rel-ative kinematics takes the form of a Lyapunov-like set of equa-tions, which is inherently ill-posed. In pursuit of unique solu-tions, we propose leastsquares algorithms, which can be solved undercertain assumptions. In Section 5 , we alsopropose similar algorithms for estimating the absolute kinematics of the nodes, givenknownreferenceparametersinthecluster.Tobenchmarkthe performanceofourestimators,wederive constrainedCramér-Rao bounds(CRBs),underaGaussiannoiseassumptiononthedata.An optimalchoice oftheweighting matrix ensures theproposed es-timatoristhebestlinearunbiasedestimator(BLUE)forthegiven datamodel.Inaddition,unconstrainedoracle boundsarealso de-rivedinSection 6 ,asabenchmarkfornextgenerationestimators. InSection 7 ,weconductexperimentstovalidatetheperformance oftheproposedestimators.

1.4.Notation:

The element-wise matrixHadamard product is denoted by  and(· )N _{denotes element-wisematrixexponent.The Kronecker}

product is indicated by , the transpose operator by (· )T _and

ˆ

(

·

)

denotes an estimated value. A vector of ones is denoted by

1N∈RN×1,IN isanN× Nidentitymatrix,0M,NisanM× Nmatrix

ofzerosand



_·



is theEuclideannorm.Foranyvector a,diag(a) isadiagonalmatrixcontainingtheelementsofaalongthe diag-onal.The blockdiagonalmatrixA=bdiag

(

A1,A2,...,AN

)

consists

ofmatricesA1,A2,...,AN alongthediagonalandzeroselsewhere.

Thefirstandsecondderivativesareindicatedby

(

·˙

)

and

(

¨·

)

respec-tively,andmoregenerallythemthorderderivative isrepresented by (_{· )}(m)_{. Unless otherwisenoted, (· )is used to indicate}

param-etersof therelative kinematicmodel.For matricesofcompatible dimensions,wewillfrequentlyusethefollowingproperties

vec

(

ABC

)

=

(

CTA

)

vec

(

B

)

, (1)

vec

(

A

)

=Jvec

(

AT

)

, (2)

whereJisan orthogonalpermutationmatrix.Wedefine an N di-mensionalcenteringmatrixasP=IN− N−11N1TN.Forasetofn

el-ements, the number of k-combinations is given by the binomial coefficient,whichisdefinedas



n k



=n

(

n− 1_k

₍

_k

)

· · ·

(

n− k+1

)

− 1

)

· · · 1 . (3)

AlistoffrequentlyusednotationsisgiveninTable 1 .

2. Time-varyingdistancesandnodekinematics

We begin by modeling the relationship between the time-varying distances, the time-varyingpositions and the node kine-matics.InSection 2.1 ,weexpand thetime-varyingpositionusing a Taylorseries, thecoefficients ofwhich yield the absolutenode kinematics.Asanextension, wepresentanovelrelative kinemat-icsmodelinSection 2.2 .InSections 2.3 and2.4 ,the relationship betweenthetime-varyingdistancesandthenodekinematicsis de-rived.Usingthesedefinitions,weformalizetheproblemstatement inSection 2.5 .

Table 1 Notations.

Notation Description

P Number of dimensions

N Number of nodes ( N > P ) D(t) ∈ R N×N _{Euclidean distance matrix at time t}

S(t) ∈ R P×N _{Absolute positions at time t}

S(t) ∈ R P×N _{Relative positions at time t}

X ∈ R P×N _{Absolute instantaneous positions at time t 0}

X ∈ R P×N _{Relative instantaneous positions at time t 0}

Ym ∈ R P×N _{m th order absolute kinematics at t 0}

Y_{m ∈ R} P×N _{m th order relative kinematics at t 0}

Hm ∈ R P×P Rotation matrix of the m th order kinematics hm ∈ R P×1 _{Translational vector of the m th order kinematics}

2.1.Absolutekinematics

Consider a cluster of N mobile nodes in a P-dimensional Eu-clideanspace(N>P),whosepositionsattimetaregivenbyS

(

t

)

∈ RP×N_{.Forasmalltime interval}

_

_{t=t− t}

0 aroundt0, weassume

that the time-varying position is continuously differentiable and thatthederivativeexistsintheinteriorofthisinterval.Therefore, the time-dependentposition vectors ofthe respective nodes can beexpandedusingaTaylorseries,

S

(

t

)

=S

(

t

)

|

t=t0+S˙

(

t

)

|

t=t0

(

t− t0

)

+0.5¨S

(

t

)

|

t=t0

(

t− t0

)

2_{+... (4)}

where

(

S

(

t

)

,_S˙

₍

_t

₎

_,¨S

₍

_t

₎

_{,. . .}

₎

_{arethe derivativesofthetime-varying} positionvectors.NowletXS

(

t

)

|

t=t0 beaP× Nmatrixcontaining theinitialcoordinatesofthemobilenodesattimet=t0.

Further-more,let the instantaneous velocities of the nodesi.e., the first-order derivatives of the position vectors S˙

(

t

)

|

t=t0 be denoted by

Y1∈RP×N,andingeneralthehigher-orderderivativesasYm

∀

m≥

1.Then,theaboveequationsimplifiesto

S

(

t

)

=X+∞ m=1

(

m!

)

−1Ym

(

t− t0

)

m. (5)

2.2.Relativekinematics

The absolute instantaneous positions at t₌t0 are an affine

transformationoftherelativepositions,i.e.,

X=H0X+h01TN, (6)

whereX∈RP×N_{istherelativepositionmatrixuptoarotationand}

translation,H0∈RP×P is the unknown rotation and h0∈RP×1 is

theunknowntranslationofthenetwork[8] .Now,weextendthis well-knownrelativepositiondefinitiontothehigher-order deriva-tives.Forinstance,thevelocityofthenodescanbewrittenas

Y1=H1Y˜1+h11TN, (7)

where Y˜1 represents the instantaneous relative velocities of the

networkatt₌t0. The translationalvector h1 isthe group

veloc-ityandH1 istheunique rotationmatrix oftherelative velocities

[18] . More generally,the mth order derivative is an affine model definedas

Ym=HmY˜m+hm1TN. (8)

We now define the relative time-varying position as S

(

t

)

=

HT

0S

(

t

)

P, and substituting the affine expressions (6) and (8) in

(5) wehave S

(

t

)

=HT 0XP+ ∞  m=1

(

m!

)

−1HT 0HmY˜mP

(

t− t0

)

m, (9)

wherewe exploitthepropertyP1N=0N to eliminatethe

transla-tionvectors,andenforcetheorthonormalityoftherotationmatrix

i.e.,HT

0H0=IN.Observethatthetranslationvectorh0 doesnot

af-fect theaboveequation. Secondly,fora meaningfulinterpretation oftherelativetime-varyingposition,areferencecoordinatesystem mustbechosene.g.,H0=IP.Insummary,withoutlossof

general-ity,weassume

H0=IP and h0=0P. (10)

andsubsequently(9) simplifiesto

S

(

t

)

=X+∞ m=1

(

m!

)

−1Y_m

(

t− t0

)

m, (11)

whereYm is therelative kinematics matrixof themthorder

de-fineduptoarotation.Inderiving(11) ,weusethefollowing prop-erties

X=XP=XP, (12a)

Ym=HmY˜m=YmP, (12b)

S

(

t

)

=S

(

t

)

P. (12c)

Note that (11) represents the relative counterpart of the ab-solute Taylor expansion (5) , where the

(

X,Y₁,Y₂,...

)

denote the relative kinematics ofthe corresponding absolute kinematics

(

X,Y1,Y2,...

)

.Ourquestinthisarticleistoestimate therelative

andabsolutekinematicmatrices,giventime-varyingpairwise dis-tancemeasurements betweenthenodes. Consequently,the abso-lute positionS(t) andrelative positionS(t) canthen beestimated using(5) and(11) respectively.

2.3. Time-varyingdistances

Similar to the node positions, the pairwise distances are also time-varyingwhichwe denotebythetime-varyingEuclidean dis-tancematrix(EDM) D

(

t

)

[di j

(

t

)

]∈RN×N wheredij(t) isthe

pair-wiseEuclideandistancebetweenthenodepair(i,j)attimeinstant t.Moreexplicitly

(

D

(

t

))

2=

ζ

(

t

)

1T

N+1N

ζ

T

(

t

)

− 2ST

(

t

)

S

(

t

)

, (13)

where

ζ

(

t

)

₌diag

(

ST

₍

_t

₎

_S

₍

_t

₎₎

_{. Observe that D(t) is a non-linear}

function oftime t,even when thenodes are inindependent lin-earmotionandhenceD(t)isacontinuouslydifferentiablefunction intime.Now,basedonthetime-varyingEDM D(t),we definethe doublecenteredmatrixB(t)

B

(

t

)

−0.5P



D

(

t

)



2P, (14a)

andthetimederivativesofthedoublecenteredmatrixas,

˙

B

(

t

)

−P



D

(

t

)

D˙

(

t

)



P, (14b) ¨B

(

t

)

−P



D

(

t

)

D¨

(

t

)

+

(

D˙

(

t

))

2



P, (14c)

where

(

D˙

(

t

)

,_D¨

₍

_t

₎

_{,. . .}

₎

_{arethederivativesofthetime-varyingEDM,} which indicate the radial velocity and other higher-order deriva-tives.Now,lettheEDMandthecorrespondingderivativesatt=t0

bedenotedby D

(

t

)

|

t=t₀R=[ri j],D˙

(

t

)

|

t=t₀R˙ =[r˙i j],D¨

(

t

)

|

t=t₀

¨R=[¨ri j],

∀

{

i,j

}

≤ N,thenwithanabuseofnotation(14) becomes B(0)B

(

t

)

|

t=t0 =−0.5PR 2_{P, (15a)} B(1)_B˙

₍

_t

₎

_|

_t=t 0 =−P



R_R˙



_{P, (15b)} B(2)¨B

₍

t

)

|

t=t0 =−P



R¨R+_R˙2



_{P, (15c)}

and higher-order derivatives can be defined along similar lines. In general, giventhedistance derivativesatt0,i.e., the range

pa-rameters

(

R,_R˙_{,¨R,. . .}

₎

_{, the double centered matrix B}(0) _{and the}

corresponding higher-orderderivatives

(

B(1)_,B(2)_,...

)

canbe con-structed. In a mobile network, the rangeparameters may not be available, however given all the nodes are capable of two-way ranging, the range parameters can be estimated using dynamic ranging[14] .

2.4. Model

To understand therelationship between the time-varying dis-tancesandtherelativekinematicsofthenodes,wesubstitutethe definition of theEDM from(13) in (14a) anddifferentiate recur-sivelytoobtain

B

(

t

)

=ST

₍

_t

₎

_S

₍

_t

₎

_{, (16a)}

˙

B

(

t

)

=_S˙T

₍

_t

₎

_S

₍

_t

₎

_+ST

₍

_t

₎

_S˙

₍

_t

₎

_{, (16b)} ¨B

(

t

)

=ST

(

t

)

¨S

(

t

)

+¨ST

(

t

)

S

(

t

)

+2S˙T

(

t

)

S˙

(

t

)

, (16c)

whereweusethedefinition(12c) andintroduce

(

S˙

(

t

)

,¨S

(

t

)

,. . .

)

as thederivativesofS(t).Now,rearrangingthetermsandsubstituting thedefinitionofS(t)att=t0 from(11) ,wehave

B0B(0)=XTX, (17a)

B1B(1)=XTY1+YT1X, (17b)

B2B(2)− 2YT1Y1=XTY2+YT2X, (17c)

where we introduce the matrices(B0, B1,B2). The joint left and

right centeringusing thecentering matrix P in(14) ensures that the phase centerofthe relative kinematicmatrices(Y1, Y2) is at

0P,similartothedefinitionoftherelativepositionX.

2.4.1. Relativekinematics

Now,combining(15a) and(17a) ,wehave

B0=XTX=−0.5PR2P, (18)

andmoregenerallyforagivenM≥ 1,(17) canbegeneralizedto

BM B(M)− M_−1 m=1



M− 1 m



YTM−mYm (19a) =XTYM+YTMX, (19b)

whereB(M) _{istheMthderivative ofthedoublecenteredmatrixat}

t0, whichisgivenby (15) andYM is theMthorderrelative

kine-maticmatrix.

Remark 1. (Measurement matrix BM): We make two critical

ob-servationsonBM in(19a) .

• Firstly, note that BM is dependent on the range parameters

(

R,_R˙_,¨R,...

₎

_{viathedefinitionofB}(M) _{(15) .}

• Secondly,B0B(0) and B1B(1) can be constructed onlybased

ontherangeparameters (see(17) ).HoweverforM≥ 2,BM not

onlydependsonB(M)_{,butalsoadditionallyreliesontherelative}

kinematic matrices of order lessthan M. Hence, if the lower orderkinematicsYm

∀

2≤ m<Mareknown,thenthe

measure-mentmatrixBMcanbereconstructed.

2.4.2. Absolutekinematics

Inaddition tothe relativekinematics, (19b) can alsobe refor-mulated to estimate the absolute kinematics YM of the network.

Recallfrom(12b) ,that therelativekinematicsoftheMthorderis

Y_M=YMPundertheassumption(10) .Substitutingthisexpression

in(19b) ,wehave

BM=XTYMP+PYTMX, (20)

whichistheabsolutekinematicmodel. 2.4.3. Modelsummary

Insummary,iftherangeparameters

(

R,R˙,¨R,. . .

)

are available,

B(M) _{canbeconstructedfrom(15) .GivenB}(0)_{,weaimtosolvefor}

therelativepositionXusingtheEq. (18) ,whichweusetoestimate the higher order kinematics. For M≥ 1, the measurement matrix

BM can be constructed using B(M) and by substituting the lower

orderrelativekinematicmatricesYm

∀

2≤ m≤ M in(19a) .Finally,

giventhemeasurementmatrix,BMandanestimateofX,ourgoal

is to estimate the Mthorder relative kinematics YM and the

ab-solutekinematics YM forM≥ 1,using(19b) and(20) ,respectively.

Wenow formulate theproblemmoreconcretely inthe following section.

2.5.Problemstatement

Problemstatement: Giventhetime-varyingpairwisedistances

D(t)betweentheNnodesinaP dimensionalEuclideanspace, es-timatetherelativekinematics(X,Y1,Y2...)andabsolute

kinemat-ics(Y1,Y2...)ofthemobilenetwork.Theseestimatessubsequently

yieldtherelative(andabsolute)time-varyingpositions.

Solution:Weproposeatwo-stepsolutiontotheabove estima-tionproblem.

S1) Dynamicrangingandrelativeposition:Giventhetime-varying distancemeasurementsD(t),weemploydynamicrangingto obtaintherangeparameters(R_,R˙_,¨R,...)inSection 3 ,under theassumptionthatallthenodesarecapableof communi-catingwitheachother.Secondly,wealsoestimatetheinitial relativepositionXusing(18) .

S2) Kinematics:ThemeasurementmatrixBMcanbeconstructed

usingtheestimatedrangeparameters,andlowerorder kine-matics(19a) .GiventherelativepositionXandBMestimates,

wesolve forthe relativekinematics YM (in Section 4 ), and

theabsolute kinematicsYM (in Section 5 ), using (19b) and (20) respectively.

Finally,giventheinitialrelativepositionandthenode kinemat-ics,thetime-varyingabsoluteandrelativepositions{S(t),S(t)}can beestimatedusing(5) and(11) respectively.

3. Dynamicrangingandrelativeposition

In this section, we aim to estimate the range parameters

(

R,_R˙_,¨R,...

₎

_{, given two-way communication between the nodes} inthemobile network.In Section 3.1 ,we relate thetime-varying propagationdelay betweenthe nodes and the rangeparameters. Given this relationship, we present a dynamicranging model in

Section 3.2 ,and subsequentlypresenta closedformalgorithm to estimatetherangeparametersinSection 3.3 .Finally,weapplythe MDSalgorithmtofindtheinitialrelativepositionofthenodesin

Section 3.4 .

3.1. Time-varyingpropagationdelay

Considerapairofmobilenodescapableofcommunicatingwith eachother.Let

τ

i j

(

t0

)



τ

ji

(

t0

)

=c−1di j

(

t0

)

bethepropagation

de-layofthiscommunicationbetweenthenode pair(i,j)attime in-stantt0,where dij(t0) isthe corresponding pairwise distanceand

Fig. 1. Dynamic ranging: A generalized two-way ranging (GTWR) scenario between a pair of mobile nodes, where the nodes exchange K time stamps asymmetrically with each other [14] . The curved lines symbolize the non-linear motion of the mo- bile nodes with time. Unlike our previous models [18,19] which considered only linear motion of the nodes, in this article we consider non-linear motion of the nodes.

cis thespeed ofthe electromagneticwave in themedium. Now, forasmallinterval



t=t− t0,weassumetherelativedistanceto

beasmoothlyvaryingpolynomialoftimewhichenablesusto de-scribethepropagationdelay

τ

ij(t) attasan infiniteTaylorseries

intheneighborhoodoft0

τ

i j

(

t

)

=c−1di j

(

t

)

=ri j+r˙i j

(

t− t0

)

+¨ri j

(

t− t0

)

2+..., (21)

wheretheTaylorcoefficientsaredefinedas

ri j,r˙i j,¨ri j,...

T =diag

(

γ

)

−1

ri j,r˙i j,¨ri j,...

T , (22) and

γ

=c[0!,1!,2!,...]T_{. Here,}

₍

_r

i j,r˙i j,¨ri j,...

)

are the derivatives

of the time-varying pairwise distance dij(t) esimtated at t=t0,

which are the elements of the matrices

(

R,_R˙_,¨R,...

₎

_{, presented} earlierinSection 2.3 .Thephysicalsignificanceofthesecoefficients isasfollows.Thepairwisedistanceatt0isrij,whichis

convention-allyobtainedfromtime of arrivalmeasurements. r˙i j is theradial

velocity, typically observed from Doppler shifts,and the second-orderrangeparameter¨ri j istherateofradialvelocitybetweenthe

nodepairatt0.Wewillnowusethisrelationinascenariowhere

mobilenodesarecapableoftwo-waycommunication.

3.2.Datamodel

Consider a generalized two-way ranging scenario between a pairofmobilenodes(Fig. 1 ),wherethenodescommunicate asym-metricallywitheachother,andrecordKtimestampsoneachnode. Thetimestampsrecordedatthekthtime instant(k≤ K)atnodei andnode jare givenbyTij,k andTji,k respectively.Thenodesare

mobileduringthesetimestampexchanges,andthereforethe prop-agationdelaybetweenthenodesisunique atevery timeinstant. Withan abuse of notation, let

τ

ij,k anddij,k be the propagation

delayandthedistancebetweenthenodepair(i,j)atthekthtime instant.Thenassumingthedistanceis(approx)constantduringthe propagationtimeof themessage,thenon-relativistic propagation delayis

τ

_{i j,k=}c−1d_{i j,k=}

|

T_{i j,k− Tji,k}

|

. Now, observe that the pair-wisepropagationdelayforGTWR canalso bewritten as(21) ,by replacing t with Tij,k (or Tji,k). More concretely, the propagation

delay

τ

_ijisgivenas

τ

i j,k=

|

Ti j,k− Tji,k

|

=ri j+r˙i j

(

Ti j,k− t0

)

+¨ri j

(

Ti j,k− t0

)

2+..., (23)

where the range parameters are estimated at t0 where

Tij,k≤ t0≤ Tij,K.

Aggregating allthe Ktimestampsforeach nodepair(i,j),and populating all measurements from N¯0.5N

(

N− 1

)

unique

pair-wiselinksforanetworkofNnodes,wehave

V







IN¯1K T T2 ...

θ







⎡

⎢

⎢

⎣

r ˙ r ¨r . . .

⎤

⎥

⎥

⎦

=

τ

, (24)

where for an Lth order polynomial approximation,

θ

∈RN¯L×1

is a vector of unknown coefficients. The N¯ dimensional vec-tor r = [rij],

∀

1≤ i≤ N, j≤ i contains all the pairwise

dis-tances at t0, and vectors containing the higher-order derivatives

(

r˙_,¨r_{,. . .}

)

are similarly defined. The matrix V is a Vandermonde-like matrix defined as V=[IN¯ 1K T T2 ...]∈RN¯K× ¯NL,

where T=bdiag

(

t12,t13,...t1N,t23,...

)

∈RN¯K× ¯N and ti j=

[T_{i j,}1− t0,Ti j,2− t0,...,Ti j,K− t0]T∈RK×1 contain all the time

stamps. All the unique pairwise propagationdelays are collected in

τ

=[

τ

T

12,

τ

T13,...

τ

T1N,

τ

T23,...]T∈RNK×1 where

τ

i j=

|

tji− ti j

|

.

Our goal in the following section, is to estimate the values [r_{i j,}r˙_{i j,}¨r_{i j,...}] from (24) , which will help usconstruct the range matrices(R,_R˙_,¨R,...).

3.3. Dynamicrangingalgorithm

Inreality,thepropagationdelayiserroneousandhence,more practically(24) is

ˆ

τ

=V

θ

+

η

, (25)

where

τ

ˆ is the noisy propagation delay, and the noise param-eters plaguing the data model are populated in

η

= [

η

T

12,

η

T13, ...

η

T 1N,

η

T23,...]T∈R ¯ NK×1_{, where}

_η

i j=[

η

i j,1,

η

i j,2,...,

η

i j,K] is the

erroruniquetothenodepair(i,j).Inpractice,thenoiseisonthe time markersTij,kandsubsequentlyontheVandermondematrix,

which has been simplified under nominal assumptions to arrive at themodel (25) . The approximationsinvolved are discussed in

Appendix-A .

Now,suppose thecovarianceofthenoiseonthenormal equa-tions



E

{

ηη

T

_}

_{, (26)}

isknownandinvertible,thentheweightedleastsquaressolution ˆ

θ

isobtainedbyminimizingthefollowingl2 norm,

ˆ

θ

=argmin

θ





−1/2

₍

_V

_θ

_{− ˆ}

_τ

₎

_

2

=

(

VT

_

−1_V

₎

−1_VT

_

−1

_τ

_{ˆ, (27)}

which is a feasible solution, if K_{≥ L} for each of the N¯ pairwise links.Moregenerally,whenLisunknown,anorderrecursiveleast squares can be employed to obtain the range coefficients [18] . Given

θ

,estimatesoftherangeparametermatrices(Rˆ,_Rˆ˙_,ˆ_¨R,...)can beconstructedusing(22) andsubsequently,from(15) wehavethe followingestimates ˆ B(0)=−0.5PRˆ2P, (28a) ˆ B(1)=−P



RˆRˆ˙



P, (28b) ˆ B(2)=−P



Rˆˆ¨R +Rˆ˙2



P. (28c) 3.4. Relativeposition

Givetheinitialpairwisedistancesatt0i.e.,R,theinitialrelative

positionsXcanbedeterminedviaMDS.GivenRˆ,letBˆ0bean

ofthismatrixyields Bˆ0=Vx



xVTx,where



x isanNdimensional

diagonal matrix containingthe eigenvalues of theBˆ0 andVx the

correspondingeigenvectors.Anestimateoftherelativeposition es-timateusingMDSisthengivenby

ˆ X=argmin X



ˆ B0− XTX



s.t.rank

(

X

)

=P =



1/2 x VTx, (29)

where



x containsthefirstPnonzeroeigenvaluesfrom



xandVx

isasubsetofVxcontainingthecorrespondingeigenvectors[8] .

4. Relativekinematics

In the previous section, we estimated the range parameters given time-varying distance measurements D(t), which was the first step(S1) inourproblemstatementdescribedinSection 2.5 . Usingtheserangeparameters,weconstructedthedoublecentered matrices

(

Bˆ(0)_,Bˆ(1)_,Bˆ(2)_,...

)

(28) and estimated the relative po-sition Xˆ using MDS (29) .Given these estimates, we now aim to solve the unknown relative kinematic matricesYM using (19) , as

proposedin(S2)ofSection 2.5 .

4.1. Linearizedmultidimensionalscaling(LMDS)

Priortoinvestigatingthegeneralkinematicmodel(19) ,we re-visit aspecialcasewhenthenodesaremobile underlinear inde-pendentmotion[18] .Insuchascenario,theaccelerationandother higherorderderivativesareabsenti.e., Ym=0,

∀

m≥ 2.Therefore,

underaconstantvelocityassumption,(17b) and(17c) simplifyto

B(1)=XTY1+YT1X, (30a)

B(2)_=2YT

1Y1, (30b)

and for m≥ 3 {Bm, B(m)} definedin (19) doesnot exist [18, Ap- pendix B] .Nowsubstitutingthedefinitionofrelativevelocityfrom

(12b) and exploiting the property HT

1H1=I, we have

B(1)=XTH1Y˜1+Y˜

T

1HT1X, (31a)

B(2)=2Y˜T₁Y˜₁. (31b)

The LMDS algorithmto estimate therelative velocity (upto a translation)isthenatwostepmethodasdecribedbelow. 4.1.1. MDS-Basedrelativevelocityestimator

Firstly, the relativevelocity up to arotation andtranslationis obtainedbyminimizingthestrainfunctionusing(31b) .LetBˆ(2)be an estimate ofB(2) _{from(28c) , withan eigenvaluedecomposition}

ˆ

B(2)_Vy



yVTy,thentherelativevelocityestimateisgivenby ˆ ˜ Y₁ =argmin ˜ Y1



Bˆ(2)− 2Y˜T₁Y˜₁



s.t.rank

(

Y˜₁

)

=P =



1/2 y VTy, (32)

where



yandVycontainthefirstPnonzeroeigenvaluesand

cor-respondingeigenvectorsof



yandVyrespectively.

4.1.2. Estimatingtheunknownrotation

TheMDS-basedsolution(32) yieldstherelativevelocityuptoa rotationandtranslation,whichisnotsufficienttoreconstructthe time-varyingrelativepositionusing(9) .Toestimatetheunique ro-tationmatrix,wevectorize(31a) ,applythetransformation(1) ,and solvethefollowingconstrainedcostfunction

argmin H1



ˆ_vec

₍

_H1

₎

_{− vec}

₍

Bˆ(1)

)



2 _{s.t H}T 1H1=IP, (33) where

ˆ ₌

₍

I_N2+J

)(

Yˆ˜T1 ˆX T

)

,

{

Xˆ_,Yˆ˜₁

}

areestimatesobtainedfrom

(29) and(32) respectivelyand,Jisapermutationmatrixsuchthat

(2) holds.

Thus, under a linear motion assumption, the relative velocity

Y₁=H1Y˜1 up to a translation can bereconstructed fora general

P-dimensional scenario using the estimators (32) and (33) . It is worth noting that the LMDS solution is feasible, only under the constantvelocityassumption.Ingeneral,theassumptiononlinear motionisnotalwaysvalidandhenceweaddressthemoregeneral kinematicmotioninthefollowingsections.

4.2.Lyapunov-likeequations

Moregenerally,when thenodesare innon-linear motion,the kinematics Ym,

∀

m≥ 1exist and mustbe estimated. Tosolve for

therelativekinematicsinthisscenario,wereferbacktoour rela-tivekinematicmodel(19) .ForanyM_{≥ 1,}themodel(19b) BM=XTYM+Y

T

MX, (34)

is the relative Lyapunov-like equation [20,21] , where BM is the

N−dimensional measurement matrix and YM is the Mth order

kinematics to be estimated. As pointed out in Remark 1 in

Section 2.4 , BM can be constructed by B(M) andlower order

rel-ativekinematics

{

Ym

}

_mM₌₁−1.Theabove equationisvery similar,but

notthesameasthefollowingequations,

AH_Y+YA=B, AY+YA=0,

AY+YC=E,

whicharethe(continuous)Lyapunovequation,commutativity equa-tion [22 , chapter 4] and Sylvester equation [23,24] respectively, wheretheunknown matrix Yhasto beestimated, given A,B,C,

E. The solutions to these equations exist and dummyTXdummy-areextensivelyinvestigatedincontroltheoryliterature[25] . How-evertheLyapunov-like Eq. (34) has received relativelyless atten-tion. The Lyapunov-like equation has a straight forward solution forP=1.But,forP≥ 2,althoughageneralsolutionwasproposed byBraden[26] ,auniquesolutionto(34) doesnotexistwhichwe discussinAppendix-B .

Now,vectorizing(34) andusing(1) ,weaimtosolve

ˆ y_M=argmin y_M



(

IN 2+J

)(

I_NXT

)

y M− bM



2 =argmin y_M



Ay_M− bM



2, (35) where A=

(

IN2+J

)(

INXT

)

∈RN 2_×NP , (36a) y_M=vec

(

Y_M

)

∈RNP×1_{, (36b)} bM=vec

(

BM

)

∈RNP×1, (36c)

andJ is an orthogonal permutation matrix (2) .The matrix

(

IN

XT

)

∈RN2_×NP

isfullcolumnrank,sinceXistypicallynon-singular. However,thesumofpermutationmatrices

(

I_N2+J

)

∈RN

2_×N2 is al-waysrankdeficientbyatleast



N₂



.Hence,thematrixprimary ob-jectivefunctionAisnotfull columnrank,butisrankdeficientby atleastP¯0.5P

(

P− 1

)

,whichisdiscussedinAppendix B .In(34) , sincethetranslationalvectorsofbothXandYMareprojectedout

usingthe centering matrix P, the P¯ dependent columnsin A in-dicate the rotational degrees of freedom in a P-dimensional Eu-clideanspace.

4.3.Lyapunov-likeleastsquares(LLS)

AuniquesolutiontotheLyapunov-likeequationisnotfeasible withoutsufficient constraintson thelinearsystem (35) . LetAˆ be an estimate of A,obtained by substituting the estimatedrelative position Xˆ (29) . Similarly, let ˆbM be an estimate of bM obtained

bysubstitutingtherangeparametersandappropriaterelative kine-maticmatricesuptoorderM− 1.Thentheconstrained Lyapunov-likeleastsquares(LLS)solutiontoestimate therelative kinematic matricesisgivenbyminimizingthecostfunction

ˆ y_M ,lls=argmin_y M



_Aˆ_y M− ˆbM



2 _{s.t. ¯Cy} M=¯d, (37)

where¯Cisasetofnon-redundantconstraints.Theabove optimiza-tionproblemhasaclosed-formsolution,givenbysolvingtheKKT equations[27 ,Section10.1.1].

4.4.Weightedlyapunov-likeLS(WLLS)

Inreality,bothAandbareplaguedwitherrorsandhencethe solutiontothecost function(37) issub-optimal.LetW¯ bean ap-propriate weighting matrix on the Lyapunov-like equation, then the weighted Lyapunov-like least squares (WLLS) solution is ob-tainedbyminimizingthecostfunction

ˆ

y_M

,wlls=argmin_y M



W¯1_M/2

(

Aˆy_M− ˆbM

)



2 s.t. ¯Cy_M= ¯d, (38)

which,similarto(37) ,canbesolvedusingtheconstrainedKKT so-lutions[27 ,Section10.1.1].Anappropriatechoiceoftheweighting matrixW¯M willbediscussedinSection 6.4 .

4.5.Choiceofconstraints:Relativeimmobility

Intheabsenceofabsolutelocationinformation,aunique solu-tionisfeasibleiftherelativemotionofatleastPnodesorfeatures areinvariant(orknown)overasmalltimeduration



t.Inan an-chorless framework, a set of given nodes would have equivalent relativekinematics,iftheyareidenticalinmotionuptoa transla-tionoriftheyare immobileforthesmallmeasurement time



t. Such situations could arise, for example, in underwater localiza-tion,when a few immobile nodes could be fixed withunknown absolutelocations,whichinturncouldassisttherelative localiza-tionoftheothernodes.ForP=2,ifthefirstPnodesarerelatively immobilefor thesmall measurement time, a validconstraint for

(37) and(38) is ¯C1=

I2 −I2 0

, ¯d1=0, (39)

whichcanbeextendedforP>2andifrequired,foralarger num-ber of immobile nodes. In essence, the relative immobility con-straintreducestheparameterspaceinpursuitofauniquesolution fortheill-posedLyapunov-likeequation.

4.6.Time-varyingrelativeposition

Inthissection,wesolvedfortherelativekinematicsofmotion, usingtherangeparameters andrelative positionestimates.When thenodesareinlinear motion,the first-orderrelativekinematics canbe estimatedusing theLMDS algorithm(32,33 ). More gener-ally,forestimatingtherelativekinematicsinanon-linearscenario, wesolvetheLyapunov-likeEq. (34) usingconstrainedleastsquares (37,38 ).Substitutingtheseestimatesin(11) ,anestimateofthe rel-ativetime-varyingpositionis

ˆ

S

(

t

)

=_Xˆ_+Yˆ₁

₍

_{t− t0}

₎

_+0.5Yˆ₂

₍

_{t− t0}

₎

2_{+... (40)}

whereXˆ isa relativepositionestimate from(29) and

{

Yˆ1,Yˆ2,...

}

arethe estimates from(37) or (38) . In thefollowing section, we aimtoestimate the absolutekinematics ofthenodesand subse-quentlythetime-varyingabsoluteposition.

5. Absolutekinematics

Inthissection,we solvefortheabsolutekinematicsYM,given

BMandtherelativepositionX.Wehavefrom(20) ,

XTYMP+PYTMX=BM. (41)

Theaboveequationissimilar,butnotthesame,tothegeneralized (continuous-time)Lyapunovequation

AT_YC+CT_YA=B,

whereA,B,Careknown squarematrices[28] .We nowvectorize

(41) andaimtominimizethefollowingcostfunction

ˆ yM=argmin yM



AyM− bM



2, (42) where A=

(

IN2+J

)(

PXT

)

∈RN 2_×NP , (43a) yM=vec

(

YM

)

∈RNP×1, (43b)

andbMisgivenby(36c) .Incomparisonto(35) ,thematrix(INXT)

is replaced with(PXT_{) in (43a) . The rankof the centering}

ma-trix P is N− 1 and since X is typically full row rank, the Kro-necker product is utmost of rank NP− P. This rank-deficiency of

P isalso reflected inthe matrixA.UnlikeA whichhas P¯ depen-dentcolomns, A isrank-deficientby



P+1₂



=P¯+P. Theadditional P dependentcolumnsareperhaps notsurprising,asthey indicate thelack ofinformationonthe translationalvector, i.e.,the group centeroftheMthorderkinematicmatrix.

5.1. Generalizedlyapunov-likeleastsquares(GLLS)

In pursuit of a unique solution to the rank-deficient system

(42) , we propose a constrained generalized Lyapunov-like least squares(GLLS)toestimatetheabsolutekinematicmatriceswhich isobtainedbyminimizingthecostfunction

ˆ

yM,glls=argmin

yM



AˆyM− ˆbM



2 s.t.CyM=d, (44)

whereAˆ andbˆMareestimatesofAandbMrespectively.Thematrix

Cisasetofnon-redundantconstraints,whichwillbediscussedin

Section 5.3 .

5.2. Weightedgeneralizedlyapunov-likeLS(WGLLS)

Theperformance oftheestimatorcanbe improvedby weight-ingthecostfunction(44) ,i.e.,

ˆ

yM,wglls=argmin

yM



W1_M/2

(

AˆyM− ˆbM

)



2 s.t.CyM=d, (45)

whichyieldstheweightedgeneralizedLyapunov-likeleastsquares (WGLLS)solution[27 ,Section10.1.1],whereWMisan appropriate

weightingmatrix(seeSection 6.4 ).

5.3. Choiceofconstraints:Anchor-awarenetwork

Foran anchoredscenario,iftheMthorderabsolutekinematics ofafewnodesareknown,thentheabsolutevelocity,acceleration and higher-order derivativescan be estimated. A straightforward minimalconstraintforthefeasiblesolutionisthen

C1=

IP¯+P, 0

, (46)

wherewithoutlossofgenerality,weassumethefirstP¯+P param-etersareknown.

5.4. Time-varyingabsoluteposition

In (44,45 ), we solved for the absolute kinematics given the measurement matrix BM and the relative position, using

con-strained leastsquaresestimators. Giventheseestimates,we have from(5)

ˆ

S

(

t

)

=_Xˆ_+Yˆ₁

₍

_{t− t0}

₎

_{+ 0.5Y}ˆ₂

₍

_{t− t0}

₎

2_{+..., (47)}

whereSˆ

(

t

)

isan estimateofthetime-varyingabsoluteposition,Xˆ

isanestimateoftherelativeposition(29) ,and

{

Yˆ1,Yˆ2,...

}

arethe

absolutekinematicestimatesobtainedbysolving(44) or(45) .

6. Cramér-Raobounds

TheCramér-Raolower bound(CRB)setsalower boundonthe minimum achievable variance of any unbiased estimator. In this section,wederivetheCRBsfortheestimatedparametersbasedon the presented datamodels. In the followingsection, we willuse theseboundstobenchmarktheperformanceoftheproposed esti-mators.

6.1. Rangeparameters

Webeginbystatingthelowerboundsfortherangeparameters basedon(25) .Let

ψ

=[rT_,r˙T_,¨rT_,...]T_{,thenthecovarianceofthe}

rangeparameters

ψ

andthe correspondingestimate

ψ

ˆ i.e.,



_ψ E



(

_ψ

ˆ₋

_ψ

₎₍

_ψ

ˆ ₋

_ψ

₎

T



_{,is bounded by}



ψ ≥

(

VT



−1V

)

−1



=

⎡

⎢

⎢

⎣



r ∗ ∗ ∗ ∗



r˙ ∗ ∗ ∗ ∗



¨r ∗ ∗ ∗ ∗ ...

⎤

⎥

⎥

⎦

, (48)

where



is the covariance of the noise on the timestamps de-fined in (26) . Here, the covariance matrices

{



r,



r˙,



¨r,...

}

are

thelowestachievableboundsforthecorrespondingrange param-eters

{

r,r˙,¨r_{,. . .}

}

.The entriesnot ofinterestaredenoted by∗ and



=diag

(

γ

)

_{ I}N¯ is atransformation matrix,where

γ

is givenby

(22) . It isworth notingthat our proposed solution(27) achieves thislowerboundforanappropriateL.

6.2. Relativeposition

TheCRBontherelativepositionsy0vec(X)isgivenbythe

in-verseoftheFisherInformationMatrix(FIM)i.e.,



xE



(

yˆ0− y0

)(

yˆ0− y0

)

T



≥ F† x, (49)

whereyˆ0isanestimateoftheunknownrelativepositiony0,



xis

thecovarianceofyˆ0[18] andtheFIMFx∈RNP×NPis Fx=JTx



¯

−1

r Jx, (50)

where

_

¯_rbdiag

(



r,



r

)

,Jx istheJacobian[18, Appendix C] and



r isobtainedfrom(48) .Intheabsenceofknownanchors inthe

network,theFIMisinherentlynonlinearandhenceweemploythe Moore-Penrosepseudoinversein(49) .

6.3. Kinematics

We now derivethe lower boundsonthe variance ofthe esti-matesoftherelative kinematicsy_M=vec

(

Y_M

)

andabsolute kine-maticsyM=vec

(

YM

)

.TheGaussiannoisevectorsplaguingthecost

functions(35) and(42) aremodeledas

ρ

_M∼ N

(

Ay_M− bM,



_ρ,M

)

, (51)

ρ

M ∼ N

(

Ay_M− bM,



ρ,M

)

, (52)

where

ρ

M,

ρ

M are N2 dimensional noise vectors, and the

corre-spondingcovariancematricesareoftheform



_ρ,ME

{

ρ

M

ρ

T M

}

≈ Ay,M



¯xA T y,M+



b,M, (53a)



ρ,ME

{

ρ

M

ρ

TM

}

≈ Ay,M



¯xATy,M+



b,M, (53b) where Ay,M =

(

IN2+J

)(

I_NYT_M

)

∈RN 2_×NP , (54a) Ay,M =

(

IN2+J

)(

PYT_M

)

∈RN 2_×NP , (54b)

andanexpressionfor



b,MisderivedinAppendix C .

6.3.1. UnconstrainedCRBs

The lowest achievable variance by an unbiased estimator is givenby



y,ME



(

yˆ_M− y_M

)(

yˆ_M− y_M

)

T



_{≥ F}† y,M, (55a)



y,ME



(

yˆM− yM

)(

yˆM− yM

)

T



≥ F† y,M, (55b)

wherethecorrespondingFIMsaregivenby

F_y,M=AT



†_ρ,MA, (56a)

F_y,M=AT

_

†

ρ,MA. (56b)

ItisworthnotingthattheMoore-Penrosepseudoinverseis em-ployed since the FIM is rank-deficient, andconsequently the de-rivedbounds(55) areoracle-bounds.

6.3.2. ConstrainedCRBs

Whenthe FIMis rank-deficient,a constrainedCRBcan be de-rivedgivendifferentiable anddeterministicconstraintsonthe pa-rameters[29] .LetU¯,Ubean orthonormalbasis forthenullspace ofthe constraint matrices ¯C,C, then the constrained Cramér-Rao bound(CCRB)ontheMthorderkinematicsaregivenby



C y,ME



(

yˆ_M− y_M

)(

yˆ_M− y_M

)

T



≥ ¯U

(

U¯T_F y,MU¯

)

−1U¯T, (57a)



C y,ME



(

yˆM− yM

)(

yˆM− yM

)

T



≥ U

(

UT_F y,MU

)

−1UT, (57b)

wheretheFIMsaregivenby(56) . 6.4.ChoiceofweightingmatricesW¯M,WM

ToadmitaBLUEsolution,weusetheinverseofthecovariance matrices



_ρ,M,



ρ,Masweightstosolvetheregressionproblems (38) and(45) ,i.e., ¯ WM



ˆ † ρ,M=

(

Aˆy



ˆ¯xAˆ T y+



ˆb,M

)

†, (58a) WM



ˆ † ρ,M=

(

Aˆy



ˆ¯xAˆTy+



ˆb,M

)

†, (58b)

wherethe estimatesAˆ_y,Aˆyare obtainedby substituting YˆM from

LLS[(37) and (44) ], in (54) ,



¯ˆx is an estimate of (49) and



ˆb,M

isderived inAppendix C from appropriate rangeparameter esti-mates.