Abstract:
An approach to adaptively positioning a set of mobile routers to provide communication services to a set of clients makes use of estimated direction profiles of communication between routers and clients. The approach does not rely on a Euclidean model in which communication characteristics (e.g., signal strength, data rate, etc.) depend on distance between communicating nodes, and does not necessarily require sampling of communication characteristics in unproductive directions in order to move the routers to preferable locations.

Description:
CROSS-REFERENCE TO RELATED APPLICATIONS 
       [0001]    This application claims the benefit of U.S. Provisional Application No. 62/069,496 filed Oct. 28, 2014, which is incorporated herein by reference. 
     
    
     STATEMENT AS TO FEDERALLY SPONSORED RESEARCH 
       [0002]    This invention was made with Government support under Contract No. FA8721-05-C-0002 awarded by the U.S. Air Force and under Grant No. W911NF-08-2-0004 awarded by the Army Research Laboratory. The Government has certain rights in the invention. 
     
    
     BACKGROUND 
       [0003]    This invention relates to adaptive communication for mobile router systems. 
         [0004]    Providing radio communication services to a distributed set of clients can be challenging in the face of these clients moving in an environment and/or requiring varying communication capacity over a wireless medium, for example, for radio frequency communication over the air. One approach is to use a set of movable “routers” and to position those routers to satisfy the client requirements, or at least to approach full satisfaction of those requirements in a best possible manner according to an overall cost function that depends on the degrees to which the clients&#39; requirements are satisfied. 
         [0005]    In this description, the references to “routers” should be understood broadly as devices that route communication, without implication that they perform any particular type of processing or that they operate at any particular layer of a communication protocol stack. For example, in certain example that addresses wireless radio frequency data communication, the described routers function as WiFi access points (or repeaters) functioning at layer 2, and in other embodiments, the also router perform layer 3 functions. However, it should be understood that no particular processing of communication is essential for a device providing communication services to be considered to be a “router”. Similarly “client” should be understood broadly to be a communication node in communication with the router. For example, a client may be a computer, smartphone, or other fixed location or mobile device. 
         [0006]    Examples in which positioning mobile routers can be important is in swarm robotic systems, which perform many complex tasks through coordination, such as cooperative search of an environment, consensus, rendezvous, and formation control. For example, Google&#39;s Project Loon, Facebook&#39;s Connectivity Lab, and similar projects envision using a network of controllable routers to provide wireless communication infrastructure in remote areas of the world. Beyond simply maintaining connectivity, reliable communication may include supporting heterogeneous and possibly time-varying communication rates amongst different pairs of agents. For example, some agents may need to use the network for transmitting video while others may simply wish to transmit status information. 
         [0007]    The general problem of optimal placement of routers to service a set of spatially distributed clients has been address in previous work in Multi-Robot Coordination. Past work on this topic includes two classes of approaches. 
         [0008]    The Euclidean Disk Model approach employs Euclidean disk assumptions where signal quality is assumed to be deterministic and mapped perfectly to the Euclidean distance between the communication nodes. A Euclidean metric allows for quadratic cost for the edges of the network and enables a geometric treatment of an otherwise complex problem. However, this approach can be limited because, in reality, signal strength suffers from large variations over small displacements. 
         [0009]    More recently, efforts have focused on giving the communication quality over each link in the network a more realistic treatment by sampling the signal strength and building closed-loop controllers using this feedback. Stochastic Sampling Methods includes methods that either supplement theoretical models for signal strength with a stochastic component based on the collected samples, or use the collected samples to design stochastic gradient controllers However, in general, these approaches suffer inefficiencies associated with sampling, and associated router movement, needed to estimate the gradient, and suffer from local minimum problems inherent with such gradient approaches. 
         [0010]    In other domains, direction of arrival information has been estimated for wireless clients using a variation of a Synthetic Aperture Approach, for example, as described in International (PCT) Patent Pub. WO2015/100237, titled “LOCALIZATION WITH NON-SYNCHRONOUS EMISSION AND MULTIPATH TRANSMISSION,” which is incorporated herein by reference. However, this or other approaches to acquiring direction or arrival information have not been applied to the problem of positioning mobile routers. 
       SUMMARY 
       [0011]    In a general aspect, an approach to adaptively positioning a set of mobile routers to provide communication services to a set of clients does not rely on a Euclidean model in which communication characteristics (e.g., signal strength, data rate, etc.) depend on distance between communicating nodes, as is assumed in general in Euclidean Disk Model approaches. Furthermore, one or more embodiments of the approach do not require sampling of communication characteristics in unproductive directions in order to move the routers to preferable locations as is done in certain Stochastic Sampling approaches. 
         [0012]    In another aspect, in general, a method for adaptively locating a set of mobile routers provides wireless communication services to a plurality of clients. The method includes, at a first mobile router of the set of mobile routers, for each client of a set of clients of the plurality of clients that are in wireless communication with the first router, (a) determining a direction profile characterizing direction-dependence of one or more signal paths between the first router and the client, and determining using the direction profile directional information comprising at least one or a communication direction between the router and the client and a spread of the direction profile, (b) computing a deficiency in communication service associate with the client and the first router, and (c) providing the directional information and the deficiency in communication for determining a location for the first router. The determined location for the first router is then used to adapt the location of said router. In some examples, the router is caused to travel toward the determined location. This method can be performed at each mobile router of the set of mobile routers, and performed repeatedly to repeatedly update the location. 
         [0013]    In another aspect, in general, software stored on a non-transitory machine readable medium includes instructions that when executed on a processor of a mobile router cause the processor to perform the method. 
         [0014]    In another aspect, in general, software stored on non-transitory machine-readable media cause multiple processors, including processors of the mobile routers, to act together to adaptively updated locations of the set of mobile routers. 
         [0015]    In another aspect, in general, a mobile router is configured to perform adaptively update its location according to a method specified above. In some examples, the mobile router comprises at least one of a robot, an unmanned vehicle, and aerial vehicle. 
         [0016]    Other features and advantages of the invention are apparent from the following description, and from the claims. 
     
    
     
       DESCRIPTION OF DRAWINGS 
         [0017]      FIG. 1  is a schematic diagram illustrating mobile routers providing communication services to spatially distributed clients. 
           [0018]      FIGS. 2-4  are flowcharts illustrating a first, a second, and a third approach, respectively, to updating router locations. 
           [0019]      FIG. 5  is a diagram illustrating channel acquisition during travel of a router. 
           [0020]      FIGS. 6A-B  are timing diagrams illustrating channel estimation using reciprocity. 
           [0021]      FIG. 7  is a flowchart illustrating an approach to estimating direction and direction spread. 
           [0022]      FIG. 8  is a block diagram of a router node. 
       
    
    
     DESCRIPTION 
       [0023]    Referring to  FIG. 1 , in an example of application of one or more embodiments described below, a set of mobile routers  110 A-C provide communication services to a set of clients  120 A-E. In the example illustrated in  FIG. 1 , the mobile routers are aerial vehicles (e.g., multi-rotor aircraft (e.g., “quadcopters”), balloons, fixed wing aircraft, low earth orbit satellites, etc.), which communicate among one another along radio links  150  and communicate with clients along radio links  140 . In the illustration a representative obstruction  160  (e.g., a building) inhibits communication between certain routers and clients (e.g., between router  120 C and client  120 D). As a result of such obstructions, and for other reasons, communication capacity between router-client pairs does not depend solely on Euclidean distance between them. In this example, at any time, each client is associated with a particular router (e.g., client  120 C is associated with router  110 B), in this examples, other routers may communicate with a client (e.g., router  110 A with client  120 C over a communication path  142 ). In operation, locations of the routers  120 A-C, and association of clients with routers, is controlled to improve overall communication service to the clients according to a criterion discussed below. This control reacts to changes in client locations, changes in capacity requirements by clients, and/or changes in link capacities, for example, resulting from changing environmental conditions. 
         [0024]    It should be understood that in the example shown in  FIG. 1  the routers operate in three dimensions. However, equivalent approaches are applicable in two dimensions. Furthermore, for the sake of exposition, two-dimensional embodiments are presented below with the understanding that the manner of modification to three dimensions is evident from this Description. Also, the approaches are not limited to radio frequency techniques. For example, freespace optical communication, and acoustic communication (e.g., in the air or underwater), can be addressed using the techniques described herein. Also the nature of the communication passing between the clients and routers is not critical. For example, the approach can applied to voice (e.g., for mobile telephones), multimedia, or other data communication. 
         [0025]    In some embodiments, the overall quality of communication service is based on the service that can be provided to a particular client j by a particular router i. (To aid in the exposition below, in general, routers are index by the letter i and clients are indexed by the letter j). For a pair of nodes (i, j) the following real quantities relate to the service that can be provided between the nodes:
       q j &gt;0 represents the demanded service (e.g., required data rate, quality) by client j. As introduced above, time dependence is not explicity indicated.       
 
         [0027]    ρ ij &gt;0 represents the available service to the client j from the the router i, for example, representing the Effective Signal to Noise Ratio (ESNR), which has a direct mapping to data rate.
       max((q j −ρ ij ), 0)≧0 is the service deficiency if client j is serviced by router i.   0≦max((q j −ρ ij )/q j , 0)≦1 is the fractional service deficiency if client j is serviced by the the router i.   α j 0 is an “importance” of client j.   0≦w ij =max(α j (q j −ρ ij )/q j , 0)≦α j  is the importance weighted fractional service deficiency       
 
         [0032]    Given the importance weighted fractional service deficiencies, an overall “max-min” cost is defined as
       G=max j (min i  w ij )       
 
         [0034]    As introduced above, a goal is to reduce or minimize G by controlling the locations of the routers:
       c i  ∈           d  (in d—dimensional space) is position of the i th router.   C=(c 1 , c 2 , . . . , c K ) is the tuple of locations of K robots.
 
We assume that the locations of the routers are known a central controller of the system. As with other quantities, the time-dependence of these locations is in general omitted unless necessary to the understanding of the notation. Similarly, the locations of the clients is represented as
   p i  ∈           d  (in d—dimensional space) is position of the j th client.   P=(p 1 , p 2 , . . . , p L ) is the tuple of locations of L clients.
 
In some embodiments, the locations of the clients are known a central controller of the system. In other embodiments, locations of the clients are inferred as discussed more fully below. As with other quantities, the time-dependence of these locations is in general omitted but should be understood.
       
 
         [0039]    Another property of the system in a particular configuration relates to the directional characteristics communication between a particular client j by a particular router j. In general, in free space with no obstruction, reflections, or other factors that affect signal quality or propagation (e.g., attenuating walls, etc), the router would expect the signal from the client to be arriving in a direction θ dir  in the vector direction of p j −c i . However, in the presence other factors a router may observe a direction-dependent signal strength, which may or may not have a peak in the θ dir  direction, and that may also have multiple peaks, for example, due to multipath effects in the environment, with the direct path not necessarily corresponding to the strongest peak. As discussed below, for each a router-client pair (i, j) the following, the system maintains a direction estimate of
       f ij :[−π, π]→         , which is a function that maps direction from the client to the router to a positive real quantity. In some instances, this function is referred to below as the “direction profile”, and a normalized version is scaled to unit area such that {circumflex over (f)} ij (θ)=(1/F)f ij (θ) where F=∫ −π   π f ij (θ)dθ.
 
In various embodiments as described below, each such directionality function can be characterized by
   θ ij   max  which is an estimate of the direction to the client, for example, at the maximum of f ij (θ); and   0≦ρ ij ≦1, which is a measure spread of f ij (θ), or the certainty of the estimate θ ij   max , as the direction to the client
 
Approaches to computing f ij  and the characterizing quantities θ max  and ρ are described fully below.
       
 
         [0043]    A number of related embodiments are described below to solve the problem for updating the desired locations of the routers. Based on information known at time t, including the current positions of the routers C(t), the system determines desired new locations for the routers, C*. These approaches are outlined below, and explained more fully later in this document. 
         [0044]    Referring to  FIG. 2 , a first approach to updating the router locations uses the Euclidean distances weighted by the importance weighted fractional service deficiency, w ij , terms, which are determined in a service deficiency computation  210 . These computations may be distributed among the routers, for example, with each router i computing terms involving that router. It should be understood that the computations shown in  FIG. 2  may be distributed, or may be centralized, for example, at a node in communication with the routers (e.g., one of the routers or another type of node on the communication backbone of the routers). A router computation  220  then determines the a new tuple of locations as 
         [0000]        C* =arg min C {max p     j     ∈P(t)  min c     i     ∈C ( w   ij   2   ∥p   j   −c   i ∥ 2 )}  (1)
 
         [0000]    in which the minimization over C is over all tuples C=(c 1 , c 2 , . . . , c K ) of joint router locations, and the minimization over c i  ∈ C is effectively the assignment of the router to service the j th  client. In at least some implementations, this computation is performed at a single node to which the required quantities (e.g., the locations of the routers and clients, and the service deficiencies) have been communicated. 
         [0045]    The solution to the optimization of equation (1) is a variant of the well known and studied “K-centers” problem (also referred to as the “Facility Location” problem), which in its optimal form is known to be an NP-hard problem. Efficient solutions, including in the case of large numbers of clients, are described in D. Feldman, S. Gil, R. Knepper, B. Julian, and D. Rus. “K-robots clustering of moving sensors using coresets”, in  ICRA,  2013, and in S. Gil, D. Feldman, and D. Rus. “Communication coverage for independently moving robots,” in  IROS,  2012, which are incorporated herein by reference. 
         [0046]    Note that an alternative to solving equation (1) in an optimal manner, a gradient approach may be used in which the gradient of the term in the brackets { } is computed with respect to each c i , and the routers are then moved in the most negative gradient direction. 
         [0047]    Referring to  FIG. 3 , second approach uses the spread/certainties j th  of the estimated directions in a direction profile computation  310 , in order to compute (d×d) weighting matrices determined in a direction weighting computation  320  and then computing the updated router locations using a modified form of the computation  220  of  FIG. 2  as a location computation  320 : 
         [0000]        C*= arg min C  {max p     j     ∈P(t)  min c     i     ∈C (w ij   2 ( p   j   −c   i ) T    M   ij ( p   j   −c   i ))}  (2)
 
         [0048]    Referring to  FIG. 4 , a third approach shares much of the computation of the second approach shown in  FIG. 3 . In this approach, each router C=arg min C  {max p     j     ∈P(T)  min c     i     ∈C (w ij   2 (p j −c i ) T  M ij (p j −c i ))} infers locations of each of the clients j in a client position inference computation  415 , with the inferred locations represented as
       p′ ij , which is an estimated or proxy location of client j from the perspective of router i, in some embodiments determined from f ij  and w ij .       
 
         [0050]    The updated locations are then determined using the router location computation  320  (as in  FIG. 3 ) as 
         [0000]        C*=arg min   C  {max j∈[1, . . . , L]  min c     i     ∈C ( w   ij   2 ( p′   ij   −c   i ) T    M   ij ( p′   ij   −c   i ))}  (3)
 
         [0051]    As introduced above, the second and the third approaches make used of direction spread/certainties ρ ij , and in the case of the third approach, also make use of the direction estimates θ ij   max . 
         [0052]    As introduced above, alternative versions of these three approaches make use of computation of the gradient of the bracketed term ∇ c     i    { } and directing the router i in the most negative gradient direction. In the case of the approach corresponding to equation (2), the most negative gradient direction for a router i will result in 
         [0000]        c*   j   =c   i ( t )+α2 w   ij   M   ij ( p   j   −c   i ( t ))
 
         [0000]    for a particular client j, and an optional positive value α, which determines a length of the position updates of the routers. 
         [0053]    A number of approaches for computing a direction profile f ij (θ) for a particular router-client (i, j) pair, for use in determining C* at time t make use of multiple transmissions between that router-client pair prior to time t. In some such embodiments, these transmissions occur during the travel of router from location c i (t−Δt) to location c i (t). In general, the system determines complex channel estimates h ij (t) (i.e., representing magnitude and phase) at an increasing sequence to times t 0 =t−Δt, t 1 , t 2 , . . . , t m =t (which are in general different for each (i, j) pair). Using an approach described in detail in International (PCT) Patent Pub. WO2015/100237, titled “LOCALIZATION WITH NON-SYNCHRONOUS EMISSION AND MULTIPATH TRANSMISSION,” the direction profile is formed as 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       f 
                        
                       
                         ( 
                         
                           θ 
                           + 
                           
                             θ 
                             i 
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                          
                         
                           
                             ∑ 
                             
                               
                                 τ 
                                 = 
                                 
                                   t 
                                   0 
                                 
                               
                               , 
                               
                                   
                               
                                
                               … 
                                
                               
                                   
                               
                               , 
                               
                                   
                               
                                
                               
                                 t 
                                 m 
                               
                             
                           
                            
                           
                             
                               h 
                                
                               
                                 ( 
                                 τ 
                                 ) 
                               
                             
                              
                             
                                
                               
                                 
                                   - 
                                   j 
                                 
                                  
                                 
                                   
                                     2 
                                      
                                     
                                         
                                     
                                      
                                     π 
                                   
                                   λ 
                                 
                                  
                                 
                                   d 
                                    
                                   
                                     ( 
                                     τ 
                                     ) 
                                   
                                 
                                  
                                 cos 
                                  
                                 
                                   
                                       
                                   
                                    
                                   
                                       
                                   
                                 
                                  
                                 θ 
                               
                             
                           
                         
                          
                       
                       2 
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
         [0000]    where assuming that the path from c i (t−Δt) to c i (t) is a straight line at an absolute angle θ i =∠(c i (t)−c i (t−Δt)), d(τ) is the scalar distance along that line. Note that it is not essential that the router travel in a straight line, and vector displacements {right arrow over (d)}(τ) can be used instead of scalar displacements in an expression {right arrow over (d)}(τ)·{circumflex over (v)} θ  instead of d(τ)cos(θ), where {circumflex over (v)} θ  is a unit vector in the θ direction. Also note that in situations in which the router has multiple antennas (e.g., for a Multiple-Input Multiple-Output, MIMO, communication approach), each antenna contributes separately to the direction profile. 
         [0054]    As introduced in Patent Pub. WO2015/100237, the routers and clients are independent in that they have local oscillators that, in general, exhibit frequency offsets. Unfortunately, even a small offset Δ f  in the frequency of the oscillators introduces a time varying phase to the wireless channel. For instance, let h(t 0 ), h(t 1 ), . . . , h(t m ) be the actual wireless channel from the robotic client to the robotic router at times t 0 , t 1 , . . . , t m . The channel observed by the router (to with a constant unit magnitude term) from its client ĥ(t 0 ), ĥ(t 1 ), . . . , ĥ(t m ) are given by: 
         [0000]        ĥ ( t   0 )= h ( t   0 ),  ĥ ( t   1 )= h ( t   1 ) e   −2πΔ     f     (t     1     −t     0     )   , . . . , ĥ ( t   m )= h ( t   m ) e   −2πΔ     f     (t     m     −t     0     ) . 
         [0055]    One approach to addressing this time-varying phase is presented in Patent Pub. WO2015/100237 through the use of a stationary antenna in conjunctions with a moving antenna, each demodulated by a same local oscillator. In the use cases introduced above, in general, a moving router cannot take advantage of such a stationary antenna. Note that in general, the phase offsets will not be precisely the result of a constant frequency difference Δf, and the terms 2πΔf(t m −t 0 ) are replaced with φ m  without any change in the procedure. 
         [0056]    In another approach to addressing the time varying phase is rather than having a fixed antenna, each router has a rotating antenna, and a fixed antenna that is maintained in a fixed relative position to center of rotation of the rotating antenna. In the case that the relative position is in or fixed relative to the direction of travel of the router, then the prior approach can be used directly. Similarly, two or more moving antennas on the router (e.g., two rotating antennas 180 degrees offset in rotation) can be used to acquire transmissions from the clients to determine directional information while compensating for frequency and phase offset between the router and the clients. 
         [0057]    Direct application of the techniques presented in WO2015/100237 that use moving antennas on the routers may increase the complexity, cost, weight, etc., of the router, and also do not directly take advantage of the motion the router itself A new approach permits the router to have a single antenna (or multiple antennas is a configuration that has a relatively fixed orientation) and to take advantage of the motion of the router itself, while also addressing the frequency (and phase) offsets described above. It should be understood that although this approach is particularly applicable to the case of controlling mobile routers as described herein, the approach has much wider applicability, for example, to situations described in Patent Pub. WO2015/100237, including in situations that address multipath and identification of most direct directions of arrival. 
         [0058]    Referring to  FIG. 5 , in an example of one update period for one router  110  moving in a direction θ i  and one client  120 , the router and client have m transmission exchanges, as introduced above. In this example, at each time t i , the router makes a transmission  510  to the client, and the client replies with a transmission  520 , for example, in the form of an acknowledgment of the transmission from the router. The client performs a channel estimation of the channel from the router to the client, for example, based on a known complex signal A being modulated than transmitted to the client where it is received and demodulated as a complex signal B. The channel from the router to the client is estimated by the client to be h=B/A. Similarly, the transmission  520  from the client to the router, the client sends a known signal C, which is demodulated as D, with the router estimating the reverse channel as h r =D/C. Although the magnitude of these channels may differ, for example, because of different signal powers emitted by the router and by the client, one would expect that the phases would be the same, or at least differ by a constant that depends on the transmit and receive chains. However, the local oscillators (LOs) of the router and client are not assumed synchronized. In this example, we assume that the LO of the router  110  generates as signal 
         [0000]      e j2πft    
         [0000]    while the LO of the client generates a signal 
         [0000]      e j(2π(f+Δf)t+ψ)    
         [0000]    for a frequency offset Δf, and a phase offset ψ, which is not assume constant between the pairs of transmissions  510 - 520 . 
         [0059]    Referring to  FIG. 6A , in the case of assuming that there is only an unknown phase ψ≠0 but no frequency offset so that Δf=0, a packet  610  with the known signal A  612  (e.g., a “pilot tone”) is received as packet  620  and demodulated as a signal B  622  as 
         [0000]    
       
      
       B=|h 
       f 
       |Ae 
       −j(2πfΔT+ψ) 
       =|h 
       f 
       |Ae 
       −j2πfΔT 
       
         e 
       
       −jψ 
      
     
         [0000]    where |h f | is the magnitude of the channel, and ΔT is the propagation time between the router  110  and the client  120 . Therefore, the client estimates the channel from the router to the client as 
         [0000]    
       
      
       h 
       f 
       =|h 
       f|e 
       −j2πfΔT 
       e 
       −ψ 
      
     
         [0000]    which as can be seen differs from a desired channel estimate by e −jψ . 
         [0060]    Within a fixed (or alternatively known) time δt, client sends a packet  630 , with a known signal C  632 , to the router, which receives it as packet  640 , and demodulates the received symbol D  642 . We assume in this example, that the client  120  communicates its estimate h f , or one or more quantities that permit the router to compute h f  (e.g., by transmitting B or a function of B). The router then makes an estimate of the reverse channel as 
         [0000]    
       
      
       h 
       r 
       =|h 
       r 
       |e 
       −j2πfΔT 
       e 
       +jψ 
      
     
         [0000]    which as can be seen differs from a desired channel estimate by e +jψ . 
         [0061]    In this example, the router uses the computed reverse channel, and an assumption that |h r |=γ|h f | and the client&#39;s estimate of h f  to compute a channel estimate as 
         [0000]        h =√{square root over ( h   f   h   r /γ)}=√{square root over (( B/A )( D/C )/γ)}{square root over (( B/A )( D/C )/γ)}=√{square root over (| h   f   ∥h   r /γ)}e −j2πfΔT  
 
         [0000]    This channel estimate is used as h(t) in the sum of the m terms to form f(θ) as set forth above (see equation (4)). For simplicity, the quantity γ is dropped (or assumed equal to 1) because it does not enter into further computations. 
         [0062]    Referring to  FIG. 6B , if we assume that there is not only a phase difference but also a frequency offset Δf, then as illustrated, computation of the channel at the router yields an estimate that is similar, with an addition of a phase term that depends on the time between receiving of the packet  620  and transmission of the return packet  630 , as 
         [0000]        h =√{square root over ( h   f   h   r /γ)}=√{square root over (( B/A )( D/C )/γ)}{square root over (( B/A )( D/C )/γ)}=√{square root over (| h   r   ∥h   r |γ)}e −j2πfΔT   e   +δtΔf/2  
 
         [0063]    This channel estimate is used as h(t) in the sum of the m terms to form f(θ) as set forth above (see equation (4)). To the extent that δt and Δf are relatively constant over the m exchanges between the router and the client, this phase term does not affect the resulting directional profile f(θ). 
         [0064]    As an alternative to computing f(θ) as described above, a “squared” version can be computed as 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       f 
                        
                       
                         ( 
                         
                           θ 
                           + 
                           
                             θ 
                             i 
                           
                         
                         ) 
                       
                     
                     = 
                     
                       
                          
                         
                           
                             ∑ 
                             
                               
                                 τ 
                                 = 
                                 
                                   t 
                                   0 
                                 
                               
                               , 
                               
                                   
                               
                                
                               … 
                                
                               
                                   
                               
                               , 
                               
                                   
                               
                                
                               
                                 t 
                                 m 
                               
                             
                           
                            
                           
                             
                               
                                 h 
                                 2 
                               
                                
                               
                                 ( 
                                 τ 
                                 ) 
                               
                             
                              
                             
                                
                               
                                 
                                   - 
                                   j 
                                 
                                  
                                 
                                   
                                     4 
                                      
                                     
                                         
                                     
                                      
                                     π 
                                   
                                   λ 
                                 
                                  
                                 
                                   d 
                                    
                                   
                                     ( 
                                     τ 
                                     ) 
                                   
                                 
                                  
                                 cos 
                                  
                                 
                                     
                                 
                                  
                                 θ 
                               
                             
                           
                         
                          
                       
                       2 
                     
                   
                   , 
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
         [0000]    thereby avoiding the necessity of computing the square root above. 
         [0065]    In some implementations, the profile f(θ) is evaluated only at a representative set of D directions (e.g., θ k =2πk/D) by accumulating the complex sums at each of these directions, and after the m transmissions, determining the D magnitudes. 
         [0066]    It should be understood that this way of using reciprocity is only one example. For instance, using m interchanges between a router and a client, the client may send C=B and the router the uses √{square root over (D/A)} to compute the estimated channel. In other examples, a router does not necessarily have to have a direct interchange with a client to estimate the channel using reciprocity. For example, a second router may overhear the transmission of C (i.e., receiving D′) and knowledge (e.g., prior knowledge, or by receiving them from the first router) of the quantities A, B, C, D as well as an exchange of messages with the first router to determine a phase correction for a channel estimated using D′/C. 
         [0067]    As introduced above, for example as illustrated in  FIGS. 3 and 4 , for each router-client pair (or at least a subset of pairs that are in communication with each other—the other pairs are assumed to be inaccessible to one another) a direction and spread/certainty of the direction from the router to the client is determined in a direction computation  310 . Referring to  FIG. 7 , and implementation of this computation is described in more detail. 
         [0068]    As discussed above, each router i and for a number of clients j, the router collects channel estimates h ij  at a number of points along its trajectory. For a particular (i, j) pair, the router collects m channel estimates h ij  at corresponding distances d i  along its trajectory, for example, as illustrated in  FIG. 5 . The channel estimates are used for form a direction profile. In some embodiments, the approach described above is used to compute f ij (θ) at a discretized set of angles. 
         [0069]    In the implementation shown in  FIG. 7 , each direction profile evaluated to determining the maximum power direction in a direction computation  730 . In one embodiment, this direction, referred to as θ max  is the angle at which f(θ) reaches its maximum. It should be understood that there are alternative ways of determining this maximum direction. For example, f(θ) may be normalized to form a circular distribution, and the first order moment of the distribution use to determine the maximum direction. As another alternative, the direction may be selected as a local maximum of the direction profile, and a most direct path may be selected by identifying maxima associated with indirect (i.e., reflected) paths using techniques described in WO2015/100237. 
         [0070]    Each profile is also evaluated to determine a quantity characterizing a spread or certainty of the determined maximum direction, with is denoted using the variable σ. One approach to computing this quantity uses the following sequence of calculations: 
         [0000]    
       
         
           
             
               
                 
                   
                     
                       
                         1 
                         ) 
                       
                        
                       
                           
                       
                        
                       
                         
                           f 
                           ~ 
                         
                          
                         
                           ( 
                           θ 
                           ) 
                         
                       
                     
                     = 
                     
                       
                         
                           ( 
                           
                             1 
                             / 
                             F 
                           
                           ) 
                         
                          
                         
                           f 
                            
                           
                             ( 
                             θ 
                             ) 
                           
                         
                          
                         
                             
                         
                          
                         where 
                          
                         
                             
                         
                          
                         F 
                       
                       = 
                       
                         
                           ∫ 
                           θ 
                         
                          
                         
                           
                             f 
                              
                             
                               ( 
                               θ 
                               ) 
                             
                           
                            
                           
                              
                             θ 
                           
                         
                       
                     
                   
                    
                   
                     
 
                   
                    
                   2 
                 
                 ) 
               
                
               
                   
               
                
               σ 
             
             = 
             
               
                 ∫ 
                 θ 
               
                
               
                 
                   
                     ( 
                     
                       θ 
                       - 
                       
                         θ 
                         max 
                       
                     
                     ) 
                   
                   2 
                 
                  
                 
                   
                     f 
                     ~ 
                   
                    
                   
                     ( 
                     θ 
                     ) 
                   
                 
                  
                 
                   
                      
                     θ 
                   
                   / 
                   
                     1 
                     
                       2 
                        
                       
                           
                       
                        
                       π 
                     
                   
                 
                  
                 
                   
                     ∫ 
                     θ 
                   
                    
                   
                     
                       
                         ( 
                         
                           θ 
                           - 
                           
                             θ 
                             max 
                           
                         
                         ) 
                       
                       2 
                     
                      
                     
                        
                       θ 
                     
                   
                 
               
             
           
         
       
     
         [0000]    Alternative to the continuous integration, equivalent discretized summations may be used. Note that σ can range from close to zero for a direction profile with a single very sharp peak, to one for a uniform direction profile, to quantities greater than one for example, in the case of multiple peaks (e.g., as a result of multi paths). 
         [0071]    There are alternative approaches to computing the spread or certainty. For example, a second moment of {tilde over (f)}(θ), a circular distribution, can be used. Other alternatives may consider the “mass” for the direction profile in the vicinity of the maximum direction, or equivalent percentile ranges. 
         [0072]    As introduced above, in some embodiments, the spread quantities σ are used to form weighted distances. Recall that in the first approach (i.e., equation (1)), a Euclidean distance between a router at location c i  and a client a location p j  is computed as 
         [0000]      ∥ c   i   −p   j ∥ 2 =( c   i   −p   j ) T ( c   i   −j )
 
         [0000]    while in the second and third approaches a weighted (Mahalanobis Distance) approach quantity 
         [0000]      (c i −p j ) T  M ij (c i −p i )
 
         [0000]    is used in place of the Euclidean distance. 
         [0073]    One approach to computing M (omitting the ij subscript) uses θ max  and σ. A unit vector {circumflex over (v)} max  is in the direction between the router and the client in the (c i −p j ) direction. A unit vector {circumflex over (v)} ⊥   max  is perpendicular to {circumflex over (v)} max . The weighting matrix M is selected to not scale distances in directions perpendicular to {circumflex over (v)} max , and to scale Euclidean distance by 1/σ 2  in the {circumflex over (v)} max  direction. This type of scaling is achieved with the computation 
         [0000]      M=QΛQ T  
 
         [0000]    where 
         [0000]    
       
         
           
             Q 
             = 
             
               
                 
                   [ 
                   
                     
                       
                         v 
                         ^ 
                       
                       max 
                     
                     , 
                     
                       
                         v 
                         ^ 
                       
                       ⊥ 
                       max 
                     
                   
                   ] 
                 
                  
                 
                     
                 
                  
                 and 
                  
                 
                     
                 
                  
                 Λ 
               
               = 
               
                 diag 
                 [ 
                 
                   
                     1 
                     
                       σ 
                       2 
                     
                   
                   , 
                   1 
                 
                 ] 
               
             
           
         
       
     
         [0074]    Note that if σ=1, no directional weighting is performed. When the direction profile is peaked and σ&lt;1, the optimization preferably moves routers towards the such clients. When the direction profile is multi-peaked and σ&gt;1, the optimization de-weights distances to those clients. 
         [0075]    As introduced above, in the thirds approach, which is illustrated in  FIG. 4 , locations of clients are inferred in a computation block  415 . One approach to such inference is performed separately at each router and makes use of the following vector computation: 
         [0000]        p′   ij   =c   i ( t )+γ w   ij    {circumflex over (v)}   ij   max  
 
         [0000]    which is a point at a distance γw ij  away from the router in the maximum direction to the client. Recall that w ij  is larger if there is a greater deficiency in the service to the client, and therefore the optimization will preferentially move the router toward that client. The constant γ is chosen to be greater to cause greater changes in desired locations for the routers. 
         [0076]    As introduced above, the approach can be implemented as a two-dimensional or a three-dimensional solution. Two-dimensional solutions are appropriate, for example, for a ground based situation for example, in locating robot-based routers in a factory with many obstructions and moving robot and/or human-based clients. Three-dimensional clients are appropriate, for example, in locating aerial routers to service ground-based or aerial clients. Although specific equations are presented, in general, for a two-dimensional implementation, analogous three-dimensional forms can be used. For example, as described in WO2015/100237, three-dimensional direction profiles, which depend on two scalar angular quantities, are formed. In the three-dimensional case, in some embodiments multiple antennas and/or non-linear (e.g., arced) trajectories are used to mitigate ambiguity in the maximum angle in the direction profile. 
         [0077]    Referring to  FIG. 8 , in an embodiment of a router, a processor  810  controls operation of the router executing instructions and using data in a memory  830 . The router uses a radio  820  to communicate with the clients  120  and determines the directional information via the antenna  822  (or optionally via multiple antennas). The determination of the desired position may be performed at a central server  810 , in communication with the router, or alternatively may be performed in whole or in part in a distributed manner using the processor in the router. As the router travels, its location is tracked, for example, by the processor accepting measurements from an inertial sensor  850  (e.g., accelerometers, gyroscope, etc.) and/or from a GPS receiver  840 . In some examples, in addition to tracking clients, the routers track each other using the same or similar techniques as described above. 
         [0078]    Implementations of the approaches described above may use software stored on machine-readable media. The software can include instructions for causing processors, for example, processors in the routers and in some embodiments in the client nodes. In some implementations, some of the processing, for example channel estimation using reciprocity, is implemented in a processor in a software-based radio. In some implementations, the router location updating is performed at a central server, while in other implementations, each router computes part of the router repositioning optimization in a distributed manner. Some embodiments implement at least part of the functions (e.g., channel estimation) in hardware, for example, using Application Specific Integrated Circuits (ASICs). 
         [0079]    It is to be understood that the foregoing description is intended to illustrate and not to limit the scope of the invention, which is defined by the scope of the appended claims. Other embodiments are within the scope of the following claims.