{'fraction_non_alphanumeric': 0.03510633955178185, 'fraction_numerical': 0.006980446585296158, 'mean_word_length': 4.234615384615385, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Optimization with Zeroth-Order Oracles in Formation

In this paper, we consider the optimisation of time varying functions by a network of agents with no gradient information. The proposed a novel method to estimate the gradient at each agent's position using only neighbour information. The gradient estimation is coupled with a formation controller, to minimise gradient estimation error and prevent agent collisions. Convergence results for the algorithm are provided for functions which satisfy the Polyak-Lojasiewicz inequality. Simulations and numerical results are provided to support the theoretical results.


I. INTRODUCTION
In time varying optimisation tasks, the goal is to optimise a sequence of problems where each new objective is a variation of the previous. The time varying function can represent the position of a moving source, with measurements capturing signal strength. We assume that only measurements, with no higher order information, are available at each iteration. To compensate for the lack of gradient information, we consider a cooperating formation of agents, sharing information to minimise the time-varying function. In general, this may be posed as a sequence of independent optimization problems [1]. One could simply treat every new cost function as an entirely new optimization problem, although this may be computationally infeasible. Additionally, solving for the optimum at each iteration is unnecessary if it is sufficient to remain within some neighbourhood of the optimum at every iteration. If there is a limit on the variation of objective parameters between iterations, the solution of the previous iteration can be updated to approach of the solution of the current iteration.
At each iteration, some amount of information about the changing function f k must be measured. Here we adopt the term p-th order oracle [2] to describe the type of available information. If p = 0 then the zeroth-order oracle only makes available the current function value f k (x k ), and not any gradient or higher order information. Gradient descent makes use of a first-order oracle, Newton's method a secondorder oracle, etc. We derive an iterative approach to track the optima of a changing cost function using zeroth-order oracles, and minimal assumptions on the behavior of f k . This is similar to finite difference stochastic approximation (FDSA), except in this case the function f k can only be observed at the locations of the agents, rather than user E. Michael, T. A. Wood, C. Manzie, and I. Shames  chosen sample points. Therefore, the agents receive information from wherever their neighbours are to compute an approximate descent direction at each iteration. As the accuracy of the agent's gradient estimate is dependent on the geometry of its neighbours, we incorporate a formation control strategy to ensure the gradient estimation is accurate.
In the area of online optimization, we will give a short review of the time varying optimization problems, but largely focus on the gradient free solutions which are more relevant to our formulation. Time varying optimization problems are well studied, frequently under the name Online Convex Optimization or OCO [3] [4]. A predictive/corrective method for OCO is presented in [5], using gradient information and line search methods. Online convex optimisation with constraints is addressed by [9], with regret bounds and convergence results. These approaches use gradient information which we assume is unavailable in this formulation. The term bandit feedback is also used to describe this problem coupled with a zeroth-order oracle, as it conforms to a multi-armed bandit problem with convex costs [2]. Regret bounds assuming compactness and convexity are derived in [6], and a similar technique is used in [7] with a multi-point estimate at each iteration for bandit feedback problems. A similar technique but using only a stochastic two point sampling each iteration is derived in [8]. These results utilise random or user chosen function sampling at each iteration, are entirely centralised, and assume convexity of the unknown cost functions. A network of zeroth-order oracles localizing the source of a static scalar field is examined in [10], with existence and convergence results, by assuming the existence of controllers with given properties.
In this paper, we present a novel algorithm which combines the information from a network of zeroth oracles to optimise a time-varying cost function. We assume the agents follow single integrator dynamics, and construct a gradient estimate which uses only local information. As well, we provide a novel method of bounding the gradient estimation error, which has an interesting geometric interpretation. As such, we incorporate formation control, along with the gradient descent, to minimise the gradient estimation error. Both gradient estimation and formation control laws require only local information, leading to an entirely distributed approach. Additionally, we allow for a time-varying objective function, and the assumptions on the time-varying objective functions are only the Lipschitz continuity of the gradient and the Polyak-Łojasiewicz inequality. These assumptions are weaker than many which are used to provide the linear convergence of gradient descent algorithms [12].
The paper is organised as follows. Section II is devoted to basic assumptions on the time-varying function and agent dynamics. Section III covers the approximation of the gradient given only zeroth-order information from an agent and its neighbours and derives an error bound on the gradient approximation. Section IV introduces and unifies formation control with the minimization. Finally, simulation and conclusions are covered in Section V.

II. PROBLEM FORMULATION
Consider a network of n agents where x i k ∈ R d denotes the position of the i-th agent for i ∈ {1, ..., n} at iteration k in dimension d. Let G = (V, E) be the underlying graph of the network with the vertex set V = {1, ..., n} and the edge set E ⊆ V × V. The edge set E captures the communication topology of the network, i.e. agent i receives information from agent j if (i, j) ∈ E. Denote the neighbour set of each In this paper, we consider the edges to be bidirectional, i.e., if (i, j) ∈ E then (j, i) ∈ E. We begin with an assumption on the network structure.
Assumption 1. Assume that the network is undirected, connected, and that If the network is disconnected, then each connected subnetwork would display the same behavior as is presented in this paper. The assumption that the neighbor set has cardinality greater than or equal to the dimension is necessary for the algorithm presented, and as the authors primarily envision physical applications (2 or 3 dimensional), is not seen as a restrictive assumption.
The agents are modeled as single integrators, with dynamics At each time instance k, each agent i can measure y i k = f k (x i k ). Let X * k denote the set of minimisers of the timevarying function f k . The following assumptions hold for the functions being minimised f k .

Assumption 2. (Differentiability and Lipschitz Gradient):
The function f k : R d → R is continuously differentiable. The gradient is Lipschitz with constant L, that is there exists a positive scalar L such that, We allow for the cost function f k to change at each iteration, however we make two assumptions on the changing cost functions.
The problem of interest is given below. Problem 1. Let X * k denote the set of minimisers of the timevarying function f k . For a network of n agents modeled in (1), under Assumptions 1-4, find inputs u i k for all agents , the set of measurements available to agent i at iteration k, and a positive constant M such that

III. ZEROTH ORDER NETWORK
If, at each time step k, each agent i was able to query an oracle and receive ∇f k (x i k ), then a standard gradient descent method could be used to reach the set of minimisers. Motivated by this, we construct an approximate gradient oracle which combines the set of measurements from the agent and its neighbours to produce a descent direction at each iteration k.
Consider the directional derivative along the path from agent i to agent j ∈ N i We construct an estimate of the gradient where we are using the shortening x ji k = x j k − x i k to represent the difference vector and v ji k = x ji k /||x ji k || as the unit vector in the difference's direction. We use u, v to denote the standard inner product when superscripts make u T v cumbersome. Note that if the function f was linear, (2) would be the exact directional derivative with ji k = 0. Using the estimate v ji k , Λ i (x k ) and Assumption 2, the error term ji k is bounded by However, computing an approximation of −∇f k (x i k ) to use as a descent direction with bound-able error requires more than just information in the x ji k direction. We must use more than 1 neighbour to construct the approximation Λ i (x k ) of the full gradient ∇f k (x i k ). At each time step, agent i computes, as the approximate gradient, Note that if the sum of outer products on the left of (5) is not of appropriate rank, it cannot be inverted to estimate the gradient. Additionally, if any adjacent agents coincide, then the gradient estimate Λ i (x k ) cannot be computed. Both the rank requirement and the requirement that no agents coincide will be addressed using formation control strategies in Section IV. (5) using only x i k and x j k , j ∈ N i . In order to ensure that it is always possible to construct a gradient estimate, the formation control covered in Section IV will preventing neighbours from being co-linear or co-planar.
For a function f k and a set of agents satisfying Assumptions 1-4, the vector Λ i (x k ) as defined in (5) satisfies where, in R 2 , the error δ i k is defined to be and we have usedv ji to indicate a vector which is orthogonal to v ji .
Proof. Recall the error bound on a single directional derivative state in (4). Let a ji k := L 2 ||x ji k || and d ji k :=

Rearranging the error bound into a set of inequalities yields
In R 2 representing the space of all possible gradients, these two inequalities enclose a band of R 2 of width 2a ji k bordered by two parallel lines perpendicular to v ji k . Consider the set of inequalities from an additional neighbour l ∈ N i Then as long as v li k and v ji k are not parallel they enclose a finite area parallelogram P ijl ∈ R 2 , and ∇f k (x i k ) ∈ P ijl . Note that assuming v li k and v ji k are not parallel is equivalent to Assumption 4. The gradient ∇f k (x i k ) is inside the parallelogram P ijl because it satisfies (8) and (9). From the original definition of the gradient estimate Λ i (x k ) in (2) we have so Λ i (x k ) is inside P ijl well. Therefore, the error ||Λ i (x k )− ∇f k (x i k )|| is bounded by diameter of the smallest ball containing P ijl . The diagonals of P ijl have lengths To upper bound the error, the longer diagonal is used. For an agent i with neighbours j and l, the longest diagonal thus has length O Fig. 1. The set of feasible gradients which satisfy (8) and (9).
which is the bound used in the Theorem.
An example of the parallelogram P ijl is shown in Figure 1, This bound does not take into account where within the parallelogram the gradient estimate Λ i (x k ) falls, which may decrease the distance to the farthest point by up to a factor of 2. If the bound at each iteration k is of interest, it is straightforward to check which corners of the parallelogram the gradient estimate is farthest from. If only the two neighbours j, l ∈ N i which are used to calculate the bound in (7) are used to calculate the estimate (5), then it is straightforward to see the estimate is the center of the parallelogram and the bound is conservative by a factor of 2.
An almost identical bounding procedure for the error is possible in R n , with each neighbour specifying a pair of parallel hyper plane constraints, which given n neighbours result in an n-parallelotope. The approximate gradient formulation (5) is the same for any dimension.
Repeating this technique of partitioning the space of possible gradients with information from additional neighbours, the bound can be tightened. However, using additional agents significantly increases the computational burden, as the resulting polytope of possible gradients will have uncertain structure, and maximizing the norm under linear constraints is itself an NP-hard problem. Additionally, empirically the error bound from the complete set of neighbours largely seems to be determined by the pair of neighbours from Theorem 1. Finally, computing the bound involving only 2 neighbours is independent of the function measurements y i k , y j k , y l k . With additional neighbours forming a polytope with more facets, this computational convenience is lost.
Using the gradient approximation from (5), the agents are able to find an approximate descent direction which has an error bounded by (7). To minimise the bound (7), thereby creating a better gradient approximation, we combine formation control with the decentralised minimization.

IV. OPTIMIZATION IN FORMATION
Examining the parallelogram P ijl in Fig. 1, there are two intuitive methods to minimise the diameter of the smallest bounding ball. We can bring the parallel edges closer together, "thinning" the parallelogram, and ensure that the two bands are orthogonal, "squaring" the parallelogram. These two methods correspond to maximizing the inner product in the denominator in (7), "squaring" the parallelogram, and minimizing the distances term, "thinning" the parallelogram. The former can be achieved by keeping the vectors v li k and v ji k orthogonal, i.e. ensuring that v li k ,v ji k = 1 and the latter by keeping the agents as close as possible while maintaining a non-collision guarantee. Finally, it is critical to prevent violation of Assumption 4, where v li k ,v ji k = 0, which geometrically corresponds to both bands in Fig. 1 being parallel to each other. We use decentralised navigation functions [13], [14] to maintain a desirable formation while minimizing f k . Definition 1. Let φ i : R nid → R + be the navigation potential function for agent i where n i = |N i |+1, with the following properties: 1) The function φ i is continuously differentiable on R d .
2) The function φ i has a unique minimum, only attained when the agents are in the desired formation configuration.

3)
The function φ i is Morse (critical points are nondegenerate). 4) The function can be computed decentrally, i.e., each agent i can compute φ i (x k ) using only x i k and x j k , j ∈ N i .
Note that these navigation potential functions exclude distance based approaches such as in [15], as they are not Morse and we cannot, as of yet, prove the global convergence properties derived here. Using the decentralised navigation functions φ i (x k ), and information available locally to each agent i, the agents are able to decrease a common global potential function Remark 2. To evaluate φ i (x k ), agent i needs access only to x i k and x j k , j ∈ N i . No information about the position of all other agents is required. Consequently, i can compute ∇ i φ(x k ) using x i k and x j k , j ∈ N i . We make the following assumption throughout the remainder of the paper.
Assumption 5. The global potential function φ(x k ) is continuously differentiable, and the gradient is Lipschitz with constant L φ .
Defining the control input u i k := −α i k p i k , the dynamics are where α i k > 0 is a design constant. Define p i k for agent i at k as where Λ i (x k ) is the estimate of the gradient from (5) and λ i k ∈ [0, 1] allows the agents to "focus" on the primary goal of minimizing f k while maintaining formation. The rules for deciding the weight λ i k and constant α i k are laid out in Theorem 2.
. . , n} be positive constants. For a set of agents with dynamics as in (14), with step direction (15), define the weighting parameter λ i k λ i k := min(1, where σ is a class K function. Defineᾱ i k to bē where c is a constant. Let the design constant α i k be in the interval where L is the Lipschitz constant for the gradient of f k . Then the system is stable, and ∃k 0 such that ∀k ≥ k 0 the global potential function is bounded φ( Proof. See Appendix A. Remark 3. In light of (16) and (18), and Remarks 1 and 2, it can be seen that each agent i can compute u i k using only x i k and x j k , j ∈ N i . Let D i 0 be the projection of set {x | φ i (x) ≤ Φ i + c} onto the subspace defined by x i and x j , j ∈ N i . The boundedness of φ(x k ) corresponds to trajectories converging to D i 0 . We should choose the formation potential functions φ i and design constants Φ i , c such that having two collinear neighbours (in 2-D) or three coplanar neighbours (in 3-D) is impossible. Thus, we ensure that the matrix in (5) is fullrank.
We may also leverage the convergence of the potential function to bound the gradient estimate error, as the formation fixes the geometry of the estimation. Using the new step direction, the the modified gradient error term is where δ i k is the upper bound on the error of the original gradient estimate from (7). Define ρ i be the radius of the smallest ball which contains D i 0 and is centred at x i . By the Lipshitz gradient property of φ(x), we then have To achieve the desired cooperative tasks agent i executes the following steps at each k, (i) the gradient of f k at x i k is estimated using (5); (ii) ∇ i φ(x k ) is computed; (iii) the value of λ i k is chosen via (16); (iv) the state is updated through (15) and an appropriate choice of α i k . We conclude this section by commenting on the overall performance of the agents in tracking the minimiser(s) of f k . To this end, note that the directions p k at each iteration are still only approximations of the true gradients. The formation of zeroth-order agents cooperating is thus equivalent to individual agents querying a δ−first order oracle at each iteration k. The definition of a δ−first order oracle is given in Definition 2. Here we show that a δ-first order oracle is sufficient to converge to a neighbourhood of the minimisers X * k , using (22) to construct an error bound on the δ-first order Oracleδ With the bounds introduced in (22), we may also define a constant α for each agent, which satisfies all of the required properties for Theorem 2.
Note that theδ used in Proposition 1 includes the formation control term L φ ρ i , because it is an additional error in the gradient estimate, although it benefits the network as a whole.

Proposition 1.
If α i is chosen such that |(1 − α i s)|< 1, then an agent using the δ-first order oracle will reach an M = η0+η * 2α i s 2 + (δ i ) 2 4s 2 neighbourhood of the optimiser X * k as the time steps k → ∞.
Proof. See Appendix B.

V. SIMULATIONS
In this section we will implement, illustrate, and analyze the method described in the previous sections. We use a formation potential adapted from [13], where each agent uses the following potential function In the formation potential function given in (25), fully explored in [13], the numerator is a quadratic attraction potential to the desired difference between agents i and j. The function in the denominator β(x k ) is described as a "collision function", which is nominally equal to 1 but quickly vanishes as the agents reach a prescribed safety distance of each other or an obstacle. The decentralised formation control from [13], [14] is shown to almost always converge, except from a set of initial conditions with measure zero. We have chosen the desired displacements c ij to form a hexagon with side lengths s = 4. The gradient error bound (7) is If the agents were in a perfect hexagon formation, 2sL would be their error bound, for s the side length and L the Lipschitz constant. However, the formation maintenance is balanced with the minimization goal, so (λ 0 , Φ i ) represents the additional error introduced from the choices of nominal weight λ 0 and the acceptable deviation bound Φ i . The specific choices of these parameters, and their impacts, are examined further in this section. Each agent's eventual nearest neighbours in the hexagon are its neighbour set N .
For the function to be minimised f k , we use convex quadratic function in two dimensions, The randomly generated quadratic in the following examples is, Q = 3.89 0.45 0.45 5.86 , to 2 decimal places. To simulate a moving source, the linear term ζ(t) is used to translate the quadratic along a path in the plane at a constant speed. The nominal weighting between the formation gradient and minimization gradient was λ 0 = 1, i.e. fully weighted on the minimization. This ensures that as long as the formation is "good enough", the agents will be attaining the best gradient estimate. The class K function σ(φ(x)) used in (16)   Immediately after the random initialization, the agents are not in formation. Therefore, their individual formation potential values φ i (x k ) far exceed the prescribed upper bound Φ i , and they coalesce into formation. Once in formation, or "close enough" as determined by φ i (x k ) ≤ Φ i , the formation begins converging to the neighbourhood of the minimisers of f k . The minimization error f k (x k ) − f k (x * k ) and neighbourhood bound from (27) are shown for individual agents in Fig. 3. The agents quickly converge to formation around the minimiser, and remains within the neighbourhood, oscillating beneath the bound as the source function f k changes. If the upper bound Φ i was decreased, representing a more stringent requirement on the formation control, the formation would converge to the minimisers X * k more slowly. The choice of upper bounded is also clearly tied with the choice of the potential function φ(x). If there is a critical safety distance, between UAVs for example, then the upper bound Φ i must be chosen such that the individual safety distance corresponds to a potential function value greater than the upper bound Φ i + c.
To demonstrate the benefits of the formation control, Fig. 4 shows the trajectories of the same agents without any formation control. The configuration of the agents is significantly looser, and though minimization error in Fig. 5 converges more quickly, it is approximately an order of magnitude larger than in Fig. 3. The error of the gradient estimate is high in this case largely due to the distance between the agents being significantly more than in the formation controlled case Fig. 2.

VI. CONCLUSION
In this paper we consider a formation of agents tracking the optimum of a time varying function f k with no gradient information. At each iteration, the agents take measurements, compute an approximate descent direction, and converge to a neighbourhood of the optimum. We derive bounds on the neighbourhood of convergence, as a function of the error in the gradient estimate, using minimal assumptions on the time-varying function. As the gradient approximation is constructed in a decentralised way, formation control is  used to encourage the agents to formations which improve the gradient estimates, while not overwhelming the task of minimizing the source function f k . We show that the formation control remains within a bounded distance of the optimal formation, and the implications to the convergence of the network to the minima of f k . In the future, a more flexible formation control approach with convergence guarantees as well as hardware experiments will be investigated.

A. Proof of Theorem 2
We first show that the trajectory of the global potential function φ(x) is a sum of local information for each agent i. Then we prove that if the potential function has violated the upper bound, i.e. φ i (x k ) ≥ Φ i , then φ i (x k+1 ) − φ i (x k ) ≤ 0 along trajectories. Then we show that if φ i (x k ) ≤ Φ i , then φ i (x k+1 )−φ i (x k ) is bounded, which finally gives that φ(x k ) is bounded for all k. Let A k be a diagonal matrix such that A i = α i k . Begin with the definition of a Lipschitz continuous gradient for the global potential function φ(x) which is equivalent to the sum of the Lipschitz conditions for each of the agents i individually, although using the Lipschitz constant L φ of the global function. Given Assumption 5, each agent has all the information required to compute their local Lipschitz bound. We proceed with the analysis of the Lipschitz bound of a single agent,