Nondisturbing Extremum Seeking Control for Multiagent Industrial Systems

Industrial applications of extremum seeking control (ESC) can be a hit and miss affair. Although a gain in performance can be achieved, the dither applied to excite the system causes unwanted fluctuations in the performance of the system. The fluctuations in systems with a single extremum seeking loop are generally small. However, for systems with many extremum seeking loops, the fluctuations in each loop may add up to an intolerable amount of fluctuation in the total performance. In this article, we propose a method to cancel the dither-induced fluctuations in the overall system performance to a large extent by smartly constructing the dither signals in each extremum seeking loop using a centralized coordinator. The novelty of our method lies in the direct calculation of the dither signals that avoids the heavy computations required by other methods. Moreover, we provide a solvability analysis for the problem of cancelling dither-induced fluctuations in the total performance of the system. Furthermore, a complete stability analysis of the overall ESC scheme with dither coordination is given.


I. INTRODUCTION
I NDUSTRIAL applications of extremum seeking control (ESC) can be found in sectors that encounter a large level of uncertainty in daily operation. These uncertainties are often caused by inconsistencies in the composition of raw materials, such as in the process industry, or by changing environmental conditions, commonly encountered in the energy industry. Typical examples of industrial applications of ESC are power optimization of wind turbines [8], [12], [18] and photovoltaic arrays [7], [22], [23], oil well optimization [26], [29], and the operation of smelting furnaces [5]. For many of these applications, the system to be controlled can be divided into several subsystems, where each subsystem is locally controlled by ESC. Although controlled by individual local controllers, the subsystems may share input resources, output infrastructure, or system-wide constraints that must be honored while steering the control inputs toward their optimal values. Examples of distributed ESC for multiagent systems are presented in [9], [13], [27], [30], [39], and [40].
Although ESC comes in many forms (see, e.g., [3], [35], and [41] and references therein), if the uncertainties acting on the system are large and slowly time varying, dither is commonly used to track changes in the unknown optimal operating conditions. Too large or too fast dither signals in the input may be harmful to equipment, violate constraints, or cause unacceptably large fluctuations in the system's output. For individual systems controlled by ESC, dither is often chosen to be sufficiently small to not pose this problem, although a positive lower bound on the amount of dither commonly exists in practice due to the presence of measurement noise [14]. For large systems with many ESC loops, the combined effects and consequences of the introduced dither signals may be substantial, even if all dither amplitudes are small. For instance, for a multiagent system comprised of many parallel subsystems, the amplitude of the dither-induced fluctuation in the overall output of the system may be as large as the absolute sum of the fluctuation amplitudes of all individual subsystems [31]. Therefore, it can be challenging to choose suitable dither signals that are sufficiently large to ensure a proper tracking of changes in the optimum of each individual subsystem, and at the same time, are small enough to avoid losses and damages due to overly large fluctuations in the overall system output.
Previous works addressing dither signals concern mostly their shape, frequency, and amplitude. For instance, the use case in [36] presents results that show a relation between the size of dither amplitudes and the achieved convergence speed. Besides, the amplitudes should be large enough to satisfy the persistence of excitation requirement, and thus, ensure convergence to the optimum; see, e.g., [1] and [2]. When it comes to the dither frequencies, an important criterion is to select them such that time-scale separation is obtained between the dynamics of the system, the dither, and all other time scales of the controller; This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see http://creativecommons.org/licenses/by/4.0/ see, e.g., [14], [20], and [33]. Performance improvement of multiagent systems can also be achieved by reusing agent frequencies that are not significantly coupled with regards to their steady-state performance; see, e.g., [21].
Preliminary ideas on cancelling dither-induced fluctuations in the total system input were reported in [27] for a multiagent resource allocation problem. In that article, the cancellation is achieved by coordinating phases of individual dither signals. Some form of minimizing dither signal amplitudes while ensuring an appropriate persistence of excitation condition was presented in [2]. However, that result does not guarantee cancellation of dither-induced fluctuations. The challenge of designing dither signals to cancel fluctuations both in the total input and total output of a multiagent system, while at the same time ensuring that each subsystem is sufficiently excited, was first treated in [31]. The solution presented in that article calculates dither signals by solving a nonconvex, yet computationally feasible optimization problem. The proposed algorithm implicitly coordinates both amplitudes and phases of individual dither signals to achieve minimization of dither-induced fluctuations both in the total input and output. This approach was then successfully applied to the case of a network of subsystems sharing a common constraint [32], with practical benefits demonstrated for the case of constrained ESC optimization for an oil production system. This article focuses on minimizing dither-induced fluctuations in the total output of a multiagent system. It proposes an algorithm for computing the dither-signals for each subsystem's extremum seeking controller such that the fluctuations in the total output (defined as the sum of the outputs of individual subsystems) are as small as possible. Considerations to avoid rapid changes in dither signals, and to maintain a constant excitation level in each subsystem's input, make a complete cancellation in total output infeasible. In contrast to the method proposed in [31] and [32], we propose an algorithm for direct calculation of the dither signals, which is much lighter than the computationally heavy method proposed in these articles. Furthermore, the dither signal optimizer of [31] and [32] considers only sinusoidal dither signals, while the new proposed method holds for a larger class of dither signals. Finally, this article presents a solvability analysis for the problem of cancelling dither-induced fluctuations in the total output and a complete stability proof for each extremum seeking loop consisting of a subsystem, an extremum seeking controller, and the dither signal generated by the proposed algorithm.
This article is organized as follows. Preliminaries are given in Section II. The ESC problem is presented in Section III. Our method to solve the dither coordination problem is given in Section IV. In Section V, an example of a stabilizing extremum seeking controller is given, which is subsequently used in Section VI to show stability and convergence of the resulting extremum seeking loops. Section VII presents two use cases to illustrate our method. Finally, Section VIII concludes this article.

II. PRELIMINARIES
The sets of real numbers, the set of nonnegative real numbers, and the set of positive real numbers are denoted by R, R ≥0 , and R >0 , respectively. The sets of natural numbers (nonnegative integers) and positive integers are denoted by N and N >0 . A scalar function α : R ≥0 → R ≥0 is said to belong to class K if it is continuous, strictly increasing and α(0) = 0. A continuous function β : R ≥0 × R ≥0 → R ≥0 is said to belong to the class KL if, for each fixed s, the mapping β(r, s) belongs to the class K with respect to r, and for each fixed r, the mapping β(r, s) is decreasing with respect to s, and β(r, s) → 0 as s → ∞ [19]. Let x be a vector in R n , where n is a positive integer. The transpose of the vector x is denoted by x T . The Euclidean norm of x is denoted by x .

III. EXTREMUM-SEEKING PROBLEM FORMULATION
Consider a system consisting of N parallel subsystems. We assume that there is no interaction between the subsystems, or at least that the interaction between the subsystems is negligibly small. The dynamics of each subsystem i ∈ {1, 2, . . . , N} are given by a function f i : is a controllable input, and t ∈ R ≥0 is the time. The set of possible input values U i is assumed to be a closed interval defined by The measured output of each subsystem is given by Due to limited knowledge about the system, we consider the state x i (t) of each subsystem, the measurement noise d i (t), and the functions f i in (1) and h i in (3) to be unknown. Let us denote the true (noiseless) output of the subsystem by z i (t) = h i (x i (t)).
(4) The true output z i (t) is measure for the performance cost of the subsystem. The total performance cost of the system is defined as the sum of the true outputs of all subsystemš To achieve optimal operating conditions under steady-state conditions, we are interested in minimizing the total steadystate performance cost. To clarify what we mean with the total steady-state performance cost, let us consider fixed input values u i (t) = u i . We assume the following.
Assumption 1: For each i ∈ {1, 2, . . . , N} and each fixed In addition, we assume that the constant solution given by the mapping X i : U i → R n i is locally attractive, uniformly on its domain as well as in time.
Assumption 2: The following holds for each i ∈ {1, 2, . . . , N}. Let us define the closed ball centered around X i (u i ) with radius r i ∈ R >0 . There exists a radius ρ 0 i ∈ R >0 and a function β xi ∈ KL such that for all x i (0) ∈ B i (u i ; ρ 0 i ), u i ∈ U i , and t ≥ 0. Because all solutions of the subsystem that start sufficiently close to the constant solution X i (u i ) eventually converge to it, we say that X i (u i ) is the steady state of the subsystem. In correspondence to (4), the true steady-state output of the subsystem is given by Therefore, the total steady-state performance cost may be expressed asŽ with u = [u 1 , u 2 , . . . , u N ] T . In order to minimize the total steady-state performance cost, we have to take into account all constraints on the system. As each controllable input u i must remain within its feasible region U i in (2), we may formulate the optimization problem as follows: min with However, because the function h i and the state x i (as well as the steady-state X i (u i )) of each subsystem are unknown, the cost functionŽ in (10) is also unknown. To optimize the system's performance despite the limited knowledge about the system, we note that the optimizer of the optimization problem in (11) corresponds to the optimizers of the optimization problems min for each subsystem i ∈ {1, 2, . . . , N}; see (10). Moreover, the output y i (t) in (3) is a reasonably accurate approximation of the steady-state performance cost Z i (u i (t)) in (9) for sufficiently low noise levels and near-steady-state conditions. Therefore, we may minimize the total steady-state performance cost by using N extremum seeking controllers to minimize the steady-state performance cost of each subsystem. A depiction of the system with N parallel control loops is given in Fig. 1. Let the controllers that optimize the steady-state performance of each subsystem i ∈ {1, 2, . . . , N} be perturbation based. Although details may vary, the general idea behind a perturbationbased controller is to excite the subsystem by adding a dither signal to its input in order to retrieve sufficient information to approximate the first-and sometimes, second-order derivatives of its steady-state performance cost function Z i . Subsequently, these derivative estimates are utilized to steer the nominal input of the subsystem toward its optimal value using a gradient-based [3], [35] or Newton-based [11], [25] optimizer. For simplicity, we limit ourselves to controllers that rely on gradient-based optimization in this work. A schematic representation of a  gradient-based controller is given in Fig. 2, whereû i (t) denotes the nominal input of the subsystem and p i (t) denotes the dither signal. A detailed example of such a controller is given in Section V.

A. Dither-Induced Fluctuations
Although dither is often necessary to optimize the steady-state performance of the system, it introduces a fluctuation in the output of each subsystem. For subsystem i ∈ {1, 2, . . . , N}, let the input u i (t) of the subsystem be the sum of its nominal valuê u i (t) and its dither signal p i (t) as Assuming that the subsystem is close to steady state and that the nominal input is much slower than the dither signal, we obtain from a first-order Taylor series approximation of Z i in (9) that the subsystem's true output in (4) can be approximated by Here, Z i (û i (t)) is the nominal value of the performance cost and dZ i du i (û i (t))p i (t) is a dither-induced fluctuation. Because the total performance cost in (5) is the sum of all subsystem outputs, we obtain thať where the first summation in the right-hand side of (16) is the nominal value of the total performance cost of the system, and the second summation is the corresponding fluctuation. Because the fluctuation of the total performance cost is the sum of the fluctuations of the outputs of all subsystems, we may have that the fluctuation of the total performance cost is large even if all fluctuations of the outputs of the subsystems are small. This situation is illustrated in Fig. 3. A large dither-induced fluctuation of the performance cost is generally undesirable as it may lead to damaged equipment, violated constraints, and large variations in production rates. A large performance fluctuation may rule out ESC as a feasible optimization method for industrial applications, especially for systems for which a steady production rate is critical, such as energy production systems.
Similarly to [31] and [32], we aim to minimize the fluctuation in the total performance cost by cleverly coordinating the dither signals of the subsystems such that the fluctuations of the outputs of the subsystems cancel each other out instead of amplifying each other. Our solution in the next section is simpler to compute than the computationally heavy solutions in [31] and [32]. Moreover, we consider a much broader class of dither signals than the single-frequency sinusoids used therein.

IV. DITHER SIGNAL COORDINATION
For any subsystem i ∈ {1, 2, . . . , N}, we define the dither signal p i (t) in (14) as the linear sum of M time-varying basis functions, denoted by b j (ωt), with j ∈ {1, 2, . . . , M} as where a i,j (t) ∈ R represent the corresponding amplitudes for i ∈ {1, 2, . . . , N} and j ∈ {1, 2, . . . , M}. By scaling the input of the basis functions by the positive parameter ω ∈ R >0 , we can select the frequency of the dither signal based on the value of ω. For simplicity, the basis functions are chosen to be periodic with common period T ∈ R >0 (for ω = 1). In addition, they are designed such that for all j, k ∈ {1, 2, . . . , M} and all t ≥ 0. Hence, the basis functions have zero mean and are orthogonal. Moreover, the basis functions and their first-order time derivatives are uniformly bounded. That means that there exist constants and for all t ≥ 0 and all j ∈ {1, 2, . . . , M}. As an example, the basis functions can be chosen as b 1 Note that a wide class of dither signals can be captured by (17) given that M is sufficiently large.
Although the amplitudes a i,j (t) may change over time, we assume that this change is much slower than the variations in values of the basis functions; in the time scale of the basis functions, the amplitudes can be regarded as quasi-constant. As a measure for the magnitude of the fluctuation of the total performance cost [see (16)], we introduce where we average over one period of the basis functions to obtain a measure that represents the average magnitude of the fluctuation instead of its instant value. Note that if the amplitudes a i,j (t) and the gradients dZ i du i (û i (t)) are slowly time varying with respect to the basis functions of the dither signals. Using the approximation in (23), γ(t) is zero if the fluctuation of the total performance cost is zero, and positive otherwise. Utilizing the orthogonality property of the basis functions in (19), we obtain from (22) that For any gradients dZ i du i (û i (t)), we aim to keep the fluctuation in the total performance cost at a minimum by selecting suitable amplitudes a ij (t) that result in a low value of the measure γ(t). However, the selection of amplitudes should not have a profound effect on the optimization of the steady-state performance of each subsystem i ∈ {1, 2, . . . , N}. Motivated by the stability analysis in Section VI, the following three important requirements must be met when selecting the dither signals.
1) For a robust optimization process, there must exist a uniform positive lower bound on the level of excitation provided by the dither signals in order to estimate the gradient dZ i du i (û i (t)) that is required for the optimization of the steady-state performance of each subsystem; see Section III. 2) All dither signals are required to be small; large dither signals result in inputs with large oscillations around their performance-optimal values, which may lead to a substantial performance loss on average. 3) The rate of change of all dither signals must be small for each subsystem to remain close to steady-state operation. Despite these three requirements, there is a large freedom in selecting the amplitudes a ij (t) in (17). If the amplitudes of the dither signals are slowly time varying with respect to their basis functions, the level of excitation provided by the dither signals is given by for all t ≥ 0 and all i ∈ {1, 2, . . . , N}; see (17) and (19). To satisfy the first requirement, we fix the level of excitation of each dither signal as for all t ≥ 0 and i ∈ {1, 2, . . . , N}, where ν ∈ R >0 is a constant. Because (26) implies that |a ij (t)| ≤ ν for all t ≥ 0, i ∈ {1, 2, . . . , N}, and j ∈ {1, 2, . . . , M}, and because all basis functions are uniformly bounded [see (20)], the second requirement is also satisfied by choosing a sufficiently small value of ν. From (24), it follows that γ(t) is a function of the dither amplitudes a ij (t) and the gradients dZ i du i (û i (t)). We aim to find mappings from the gradients dZ i du i (û i (t)) to the amplitudes a ij (t) that lead to a small value of γ(t) for each possible set of gradient values. For each i ∈ {1, 2, . . . , N} and each j ∈ {1, 2, . . . , M}, let this mapping be denoted by We assume that each mapping is defined such that the fixedexcitation condition in (26) is satisfied for all possible gradient values. To satisfy the third requirement, we additionally assume that each mapping is globally Lipschitz so that for some Lipschitz constant L A ∈ R >0 . In this case, the rate of change of all dither signals can be made arbitrarily small by tuning the extremum seeking controller of each subsystem such that the rate of change of all gradients is small, and by choosing a small value of ω such that the basis functions b j (ωt) of the dither signals are slowly time varying; see (21). We note that the are not known. However, as mentioned in Section III, for each i ∈ {1, 2, . . . , N}, the extremum seeking controller of the subsystem i contains a gradient estimator that produces an estimate of dZ i du i (û i (t)) to optimize the performance of the subsystem. This estimate is also used for computing the dither amplitudes.

A. Theoretical Lower Fluctuation Bound for a Fixed Excitation Level
As a benchmark for the reduction of the fluctuation of the total performance cost that can be achieved with the described approach, we determine the theoretical minimal value of γ(t) for any (not necessarily Lipschitz) mappings in (27) that satisfy the fixed-excitation condition in (26). This theoretical lower bound on the value of γ(t) is given in the following lemma.
Under the fixed-excitation condition in (26), the minimal value of γ(t) that can be achieved for any (not necessarily Lipschitz) mappings of the dither amplitudes is given by (30) Proof: See the Appendix A. It follows from Lemma 3 that a zero total fluctuation can only be achieved whenever the magnitude of the largest gradient is smaller than or equal to the sum of magnitudes of all other gradients. Only in that case, it is possible to choose the dither amplitudes such that the fluctuations of the performance cost of individual subsystems cancel each other out. If the magnitude of the largest gradient is larger than the sum of all others, complete cancellation is not possible under the fixed-excitation condition in (26). The smallest fluctuation of the total performance cost can then be achieved by cancelling out the fluctuation of the performance cost of the subsystem with the largest gradient as much as possible by choosing the same perturbation with an effective opposite sign for all other subsystems. This approach results in the minimal value of γ(t) in (30).
The mappings of the dither amplitudes that correspond to the minimal value of γ(t) in Lemma 3 are discontinuous, and therefore, not globally Lipschitz. These discontinuities occur at the zero crossings of the gradient values. Alterations to obtain globally Lipschitz mappings are described in the next section. Note that these alternations inevitably lead to a larger value of γ(t).

B. Dither Coordination Method
In this section, we present our solution to the dither coordination problem. Our method requires that the number of the basis function M is equal to the number of subsystems N . Let I(i) denote the numbering of the subsystems from the largest to the smallest gradient value, as in Lemma 3. For each i, j ∈ {1, 2, . . . , N}, let us define the mappings where the nonnegative function ϕ ki (t) is given by for all i ∈ {1, 2, . . . , N} and some shape parameter λ 1 ∈ R >0 , and where the function s ij (t) is given by Due to the difference in magnitude of the gradients, the dither amplitudes of the various subsystems may have different effects on the fluctuation of the total performance cost. The function ϕ ki (t) in (32) is introduced to counteract these differences. To ensure that the fluctuations cancel each other out instead of amplifying each other, the function s ij (t) in (34) is introduced to give each amplitude the appropriate sign. In addition, to guarantee that the resulting dither signals satisfy the excitation condition in (26), we define the scaling parameters ξ k (t) ∈ R ≥0 for all k ∈ {1, 2, . . . , N}, such that  (33), the mappings in (31) are not globally Lipschitz (similar to the optimal mappings in Lemma 3). With χ i (t), it can be verified that the mappings of the dither amplitudes in (31) are globally Lipschitz for any positive value of the shape parameter λ 1 . By substituting a ij (t) = A ij (g Z (t)) into (24), we obtain that for all t ≥ 0 and any set of gradient values. This indicates that, for large values of λ 1 , the mappings A ij in (31) result in a small total fluctuation that is reasonably close to the minimal total fluctuation computed in Lemma 3, while selecting a small value of λ 1 may lead to a significant increase in total fluctuation of the performance cost. However, the resulting global Lipschitz constant is large for small values of λ 1 . As a result, this may require that the extremum seeking controller is tuned very conservatively to make gradients sufficiently slow to satisfy the third requirement in Section IV. Hence, there is a tradeoff to be found between maximizing the optimization speed of the total performance cost and minimizing its fluctuation when selecting the parameter λ 1 . Although the mappings in (31) are not optimal, it should be noted that the worst-case fluctuation of the total performance cost for the mappings in (31) is equal to and often much smaller than the worst-case fluctuation for fixed dither amplitudes. This claim is easy to prove by noting that, for any dither amplitudes that satisfy the excitation condition in (26), there exist gradients such that and that the right-hand side of (38) is equal to or larger than the equality (but not necessarily the inequality) in (37); see the definitions of χ i (t) and ξ k (t).
As mentioned in Section IV, estimates of the gradients need to be used because their true values are unknown. Let the vector of gradient estimates be denoted bŷ Although the true gradient values may be slowly time varying, their estimates may not be. To guarantee that the estimates are sufficiently slow, we propose the following rate limiters: for i ∈ {1, 2, . . . , N}, with where λ 2 and λ 3 ∈ R >0 are tuning parameters. Note that |ġ Z,i (t)| ≤ λ 2 for all t ≥ 0. Hence, λ 2 determines the upper limit on magnitude of the time derivative ofḡ Z,i (t). The value of the parameter λ 3 should be large so thatḡ Z,i (t) converges fast to the gradient estimateĝ Z, The dither amplitudes that are used to attenuate the fluctuation of the total performance cost are thus given by withḡ Z (t) = [ḡ Z,1 (t),ḡ Z,2 (t), . . . ,ḡ Z,N (t)] T . Note that, for all i, j ∈ {1, 2, . . . , N}, the following bound on the time derivatives of the dither amplitudes can be given: where the parameter λ 2 is defined in (39) and L A is the global Lipschitz constant of the mappings in (31).

V. EXTREMUM SEEKING CONTROLLER
The proposed dither coordination method can be applied to new and existing ESC schemes; see Section III. We introduce the following controllers as an illustration. For each subsystem i ∈ {1, 2, . . . , N}, let us consider the following gradient estimator, similar to the one in [15]: Here,ĝ Z,i (t) ∈ R is an estimate of the gradient dZ i du i (û i (t)). The signals m i (t), q 1,i (t) ∈ R and q 2,i (t) ∈ R >0 are internal variables of the gradient estimator. The tuning parameter η ∈ R >0 is chosen sufficiently small, such that the variable q 1,i (t) remains close to zero and the time variations in the (positive) variable q 2,i (t) are small. Let the optimizer for each subsystem be given byu where κ 1 and κ 2 ∈ R >0 are tuning parameters. The function proj i is the projection operator that confinesû i (t) to a set U i ⊂ U i . From (17), (20), and (26), it follows that there exists a constant c p , ∈ R >0 such that for all i ∈ {1, 2, . . . , N} and all t ≥ 0. By defining the setÛ i aŝ we ensure that, ifû i (t) ∈Û i , then u i (t) ∈ U i . Here, we assume that ν is sufficiently small, such that the setÛ i is nonempty. The projection operator proj i can now be defined as for all i ∈ {1, 2, . . . , N}. In the next section, we prove asymptotic stability under suitable conditions for each closed-loop subsystem consisting of the subsystem in (1) and (3), the controller in (43) and (44), and the dither signal in (17) with the amplitudes in (41).

VI. STABILITY ANALYSIS
In order to proof convergence of each extremum seeking loop, we require that true steady-state output of each subsystem given by the mapping Z i in (9) has a unique minimum that corresponds to the optimal performance. Assumption 4: For each i ∈ {1, 2, . . . , N}, there exist a constant u * i ∈ U i and a function α Z,i ∈ K such that for all u i ∈ U i . Hence, u * i is a unique minimizer of the function Z i in (9) on the domain U i in (2). Assumption 4 implies that, from any initial condition u i (0), we will end up at the optimal value u * i as long as we follow a gradient-descent path with the gradient-descent optimizer in (44). In addition to Assumption 4, we require that all functions are sufficiently smooth on their relevant domains. In order to define the relevant domain of the state vector x i (t), let us introduce the set where B i (u i ; ρ i ) is the closed ball defined in (7) and ρ i = β xi (ρ 0 i , 0), with β xi and ρ 0 i defined in Assumption 2. We note that X i is compact because the set U i in (2) is a closed interval. We assume the following.
Assumption 5: The following holds for each i ∈ {1, 2, . . . , N}. The function With these assumptions, we are able to prove asymptotic convergence of each extremum seeking loop to a small region of the optimal steady-state conditions, and give a bound on the gradient estimation error.
We note that we can make the bound in the right-hand side of (51) arbitrarily small in the absence of the disturbance d i (t) by choosing small values of ν and ω. In the presence of the disturbance d i (t), the bound cannot be made arbitrarily small. However, boundedness of the solutions can still be guaranteed. This stability result is similar to the ones in [14,Corollary 13] and [16,Th. 14].

A. Proof of the Stability Result
We begin the proof of Theorem 6 by deriving a bound on the magnitude of time derivative of the input u i (t). We obtain from (14) and (17) thaṫ where we have used that M = N for the dither coordination method in Section IV-B. From (20), (21), |a ij (t)| ≤ ν in (26), (44), and the bounds on the parameters in Theorem 6, it subsequently follows that |u i (t)| ≤ ωνcu (53) for all t ≥ 0, all η ≤ ω 4 , all κ 1 ≤ ην 5 and all λ 2 ≤ ην 7 (where L A in (42) is a global Lipschitz constant of the mappings in (31) for all λ 1 ≤ 6 , for any given value of 6 ). Now, let us focus our attention to the state of the subsystem in (1). Assumption 2 ensures convergence of the state x i (t) to the steady state X i (u i ) for fixed inputs u i . In the following lemma, we prove convergence to a neighborhood of the steady state for time-varying inputs u i (t). Note that the size of this neighborhood depends on how fast the input u i (t) changes over time. Therefore, to ensure convergence to a small neighborhood, we must make sure that the time derivative of u i (t) is sufficiently small.
Lemma 7: Let us definẽ Under the assumptions of Theorem 6, there exist class-K functions α x1 and α x2 , such that the solutions of (1) satisfy for all initial conditions where the values of 1 , 2 , 3 > 0 are sufficiently small to guarantee that x i (t) ∈ X i for all t ≥ 0, all ν ≤ 2 and all ω ≤ 3 by using the bound oṅ u i (t) in (53).
Proof: The proof of the lemma follows similar lines as the proofs of [37,Prop. 2] and [14,Lem. 8]. In short, first a converse Lyapunov theorem (e.g., [24] or [19,Th. 4.16]) is used to show that there exists a Lyapunov function for fixed inputs u i (t) = u i . Subsequently, this Lyapunov function is used to prove the input-to-state stability of the subsystem with respect to the time derivativeu i (t) of the input; see (55) and (56). For the construction of the converse Lyapunov function, we require that the functions f i and X i are Lipschitz on their domains. This is satisfied due to the differentiability of the functions (see Assumption 5) and the compactness of the sets X i and U i . Further details of the proof are omitted for brevity.
The state errorx i (t) in (54) influences how close the measured output y i (t) in (3) of the subsystem is to the true steady-state output Z i (u i (t)) in (9). Let us define the measurement error Because x i (t) ∈ X i and u i (t) ∈ U i for all t ≥ 0 (see Lemma 7 and Section V), it follows (3), (9), and Assumption 5 that there exists a constant cỹ ∈ R >0 , such that for all t ≥ 0. Ifx i (t) and d i (t) are small, then the measurement errorỹ i (t) is also small and the measurement output y i (t) is approximately equal to the true steady-state output Z i (u i (t)).
In turn, the signalĝ Z i (t) produced by the estimator in (43) is an accurate approximation of the gradient dZ i du i (u i (t)), as shown next.
Proof: See Appendix B.
Using an accurate gradient estimate, the gradient-descent optimizer in (44) steers the nominal inputû i (t) to a small neighborhood of the valuê see Assumption 4 and the definition of the setÛ i in (46). We obtain the following result. It follows from (61) that |u * i −û * i | ≤ νc p . Therefore, we obtain using (14) and (45) that for all t ≥ 0. Note that all solutions of the extremum seeking scheme are bounded under the conditions of Theorem 6. The

A. Power Production of Solar Arrays
Consider three photovoltaic arrays with n s cells in series affected by different sun light and temperature conditions. According to the model presented in [38], the light generated current i s and the reverse saturation current i o in a photovoltaic cell with temperature T and solar irradiance S can be described as where the parameters of the equations are given in Table I. The corresponding output current i of a photovoltaic array with n s connected cells in series is where v is the output voltage. A dc-dc buck converter is used to connect each array to a dc load. The converter dynamics are represented by the model in [23] as where i(v; T, S) is the nonlinear mapping from the duty cycle u to the output current. This mapping can be computed from (64) and (65). Three photovoltaic arrays are simulated for the temperatures T and irradiances S given in Table II. A subscript i ∈ {1, 2, 3} is used to distinguish between the three arrays. The steady-state relation between the duty cycle and the power output of each array is depicted in Fig. 4. To maximize the total produced power, we apply the extremum seeking controllers in Section V to each of the photovoltaic arrays. Here, we use the negative of the power as a measure for the performance cost (67) We apply the dither coordination method in Section IV-B to minimize the fluctuation of the total produced power. The base functions of the dither signals and the controller parameters are presented in Table II, where tri(·) is the triangle wave function defined in [36]. We constrain the duty cycle of each array to be between 0.45 and 1. Therefore, two out of the three arrays are not able to produce at their peak values in Fig. 4. As a reference, we compare the results for the varying coordinated dither signals with those of fixed dither signals (where the amplitudes of the dither signals do not change over time). The fixed dither signals have the same basis functions and effective amplitude as the varying dither signals.
We observe from the simulation results in Fig. 5 that each extremum seeking controller steers the duty cycle u to its constrained optimal value. As a result, the power output P of each photovoltaic array goes to its constrained maximum. All optimal values are indicated by the dashed lines in Fig. 5. The figure shows the results for the varying dither signals. Due to their large similarity, the corresponding results for the fixed dither  signals are omitted. The main difference between the varying dither signals and the fixed dither signals is showcased in Fig. 6; the fluctuation in the total power P total of the three arrays is visually much smaller for the varying dither signals with dither coordination than for the fixed dither signals without dither coordination. Due to the smaller fluctuation, we see in Fig. 6 that the value of the fluctuation measure γ in (24) is much smaller overall for the varying dither signals than for the fixed dither signals. Moreover, after the transient response has died out in the first two seconds, the value of γ for the varying dither signals remains close to the theoretical lower bound derived in Lemma 3. This lower bound is indicated by the dashed line in Fig. 6. Fig. 7   presents the dither amplitude a ij that are used to compute the varying dither signals. Note that there are no rapid changes in the values of a ij , which is one of the three requirements in Section IV. Thus, the proposed dither coordination method in Section IV-B is able to drastically reduce the fluctuation in the total produced power, even in the presence of dynamics and local input constraints.

B. Power Optimization in Electric Submersible Pump (ESP)-Diluent Lifted Wells
In this example we consider oil wells equipped with ESPs [34] installed several hundred meters down in the well. High viscosity of the reservoir fluid (consisting of viscous oil and water) increased by emulsion formation at the ESP, significantly reduces ESP efficiency and increases its power consumption. To reduce the ESP power consumption, diluent (lighter oil) is injected upstream the ESP [10], [26]. A diagram of this production system for the case of several wells is shown in Fig. 8.
In the initial phase of oil production, when the percentage of water in the reservoir fluid is low, increasing the diluent injection rate initially reduces the ESP power consumption; after a certain optimal rate, it starts to increase the ESP power consumption. The objective of diluent optimization is to find the optimal diluent injection rate that provides a minimal power consumption for the given ESP intake pressure setpoint.
Due to high uncertainties (emulsions can increase fluid viscosity up to six times [28]), extremum seeking control is a proper optimization tool to address this problem. However, dither-induced fluctuations in ESP power consumption from several ESP-lifted wells can add up to formidable power fluctuations negatively affecting power generators on the oil platform. The dither coordination method presented in this article can be applied to reduce these power fluctuations.
For simplicity, we present a static model of an ESP lifted well, which is derived from the dynamic model in [4] where p in is the ESP intake pressure; p wh is the well-head pressure; ρgh is the hydrostatic pressure at the ESP intake for a given fluid density ρ and a vertical distance h; Δp f is the frictional pressure loss (calculated based on the Darcy-Weisbach and Colebrook equations [17]); and Δp p is the pressure increase yielded by the pump. The production rate of the reservoir fluid q r is calculated as a product of the productivity index P I and the difference between the reservoir pressure p r and the pressure p bh at the bottom of the well. The total flow rate through the pump q p consists of the reservoir inflow rate q r and the diluent injection rate q d . The pressure increase provided by the ESP Δp p and the ESP break horse power P depend on the pump speed f , flow rate through the pump q p , the fluid density ρ, and the viscosity μ. Δp p and P are calculated from the pump characteristics H 0 (q) and P 0 (q) corresponding to the pump performance at the pump speed f 0 , the flow rate q, the fluid density ρ 0 , and the viscosity μ = 1 cP. The terms C H (μ), C Q (μ), and C P (μ) are viscosity correction factors for using the model with fluids of higher viscosities; see [4] for more details on the model.
Diluent injection enters the model through the diluent injection rate q d , the fluid viscosity μ, and the density ρ, which is calculated by mass balance equations from the densities and flow rates of the reservoir fluid and the diluent. The viscosity of the diluted fluid μ is calculated using formulas for oil-water emulsions from [6]. The effect of diluent depends on the water cut WC-The percentage of water in the reservoir fluid.
The ESP speed f is controlled by a proportional-integral controller that keeps ESP intake pressure p in at a set-point p sp in , and thus, provides a constant inflow rate of the reservoir fluid consisting of oil and water.  Based on the presented model, we calculate the curves describing the effect of diluent injection on the ESP power consumption for three identical wells with different water cuts: WC=32%, 40%, and 35%, respectively, as shown in Fig. 9. Numerical values for these curves are taken from [4] with the additional parameters being: the viscosities and densities μ d = 1 cP and ρ d = 800 kg/m 3 of the diluent, μ o = 125 cP and ρ o = 970 kg/m 3 of the oil, and μ w = 1 cP and ρ w = 1000 kg/m 3 of the water. In reality, these curves, apart from being convex, are highly uncertain due to inaccurate pump models, uncertainties in reservoir fluid properties and a varying diluent mixing efficiency at the ESP. Extremum seeking control is, thus, a suitable method to find the optimal diluent injection rate [26], although one needs to avoid or reduce the dither-induced fluctuation in the total ESP power consumption.
This can be achieved with the method presented in this article. The basis functions and parameters for the extremum seeking controller are presented in Table III. The basis functions have periods of about 15-20 min to ensure quasi steady-state operation of the well under these perturbations. The diluent injection is constrained to the interval of [100,400] m 3 /d. Fig. 10 shows the performance of the controller for the three ESP-lifted wells with varying dithers signals. As seen from the figure, the diluent injection rates u i = q d i for each well  i are steered toward the optimal solutions by the extremum seeking controllers so that the power P i consumed by each ESP reaches its minimum value. Fig. 11 presents the total power consumed by all the ESPs and the fluctuation measure γ in (24) for fixed and varying dither signals. As seen in the figure, the extremum seeking controllers with coordinated varying dither signals noticeably mitigate the fluctuation of the total power consumption compared to the standard method with fixed dither signals. The fluctuation measure γ is considerably lower for the coordinated dither signals, being much closer to the theoretical lower bound indicated by the dashed line in the figure. Fig. 12 depicts the dither amplitudes a ij calculated for the coordinated dither signals shown in Fig. 10. The absence of spikes or abrupt changes in the dither amplitudes is essential for the smooth and safe operation of this algorithm in the presented oil production system.

VIII. CONCLUSION
In this article, we have introduced a novel method to coordinate the dither signals utilized in the extremum seeking loops of individual subsystems in order to minimize the dither-induced fluctuation in the total output of a multiagent industrial system. Therefore, we do not only optimize the quantity of the total output by employing the optimizing capabilities of ESC but also its quality by minimizing the undesired dither-induced fluctuation. The presented method is computationally much lighter than current alternatives. This greatly benefits its practical deployment. The stability analysis in this article tells us that convergence to a small region of the optimum can be guaranteed under conditions that are similar to those of other ESC methods. Two case studies show that the presented dither-coordination method can drastically decrease the magnitude of the fluctuation in the total output, while having a negligible effect on the performance of the extremum seeking loop of each individual subsystem.
APPENDIX A PROOF OF LEMMA 3 For brevity of notation, we drop the time index t. The proof of the lemma consists of two parts. In the first part, we prove that, if (1) )|, it is always possible to choose the dither amplitudes such that γ = 0 and M j=1 a 2 ij = ν 2 for all i ∈ {1, 2, . . . , N}. In the second part, we prove that the minimal value of γ is (| du I(i) (û I(i) )|. Part 1: For γ to be zero, we require that Due to the numbering specified in the lemma, we have that dZ I (1) du I(1) (û I(1) ) = 0 implies that dZ I(i) du I(i) (û I(i) ) = 0 for all i ∈ {1, 2, . . . , N}. If dZ I (1) du I(1) (û I(1) ) = 0, then M(t) = 0. Therefore, any vector a col j ∈ R N is in the nullspace of M. As a consequence, any dither amplitudes that satisfy M j=1 a 2 ij = ν 2 for all i ∈ {1, 2, . . . , N} will result in the theoretically minimal value γ = 0. We can always find dither amplitudes that satisfy these conditions. Thus, the minimal value of γ is zero if dZ I (1) du I(1) (û I(1) ) = 0. Now, let dZ I (1) du I(1) (û I(1) ) = 0. The nullspace of the matrix M is given by where for each k ∈ {2, 3, . . . , N},v k is a vector with two nonzero elements: its first element is (1) ), 0, . . . , 0] T , etc.). Therefore, any vector a col j ∈ ker(M) may be written as a col j = N k=2 q jkvk for some coefficients q jk ∈ R. It follows from this and the definitions of a col j andv k that, if there exist dither amplitudes that satisfy γ = 0 and M j=1 a 2 ij = ν 2 for all i ∈ {1, 2, . . . , N}, then there must exist coefficients q jk such that M j=1 if N is even.
(76) Note that the right-hand side of (76) is smaller than or equal to ( dZ I (1) du I(1) (û I(1) )) 2 using the specified numbering of the gradients. Moreover, noting that r col k for k ∈ {2, 3, . . . , N} are N − 1 vectors on an M -dimensional unit sphere, M ≥ N − 1 implies that there exist a continuous path on the M -dimensional unit sphere between the vector values that correspond to (75) and the vector values that correspond to (76). Therefore, following this path, the value of the left-hand side of (74) will change from the value ( N k=2 | dZ I(k) du I(k) (û I(k) )|) 2 to a value that is smaller than or equal to ( dZ I (1) du I(1) (û I(1) )) 2 in a continuous manner. Thus, somewhere on this path, we will always encounter values of r col k for k ∈ {2, 3, . . . , N} for which (74) holds if and only if N k=2 | dZ I(k) du I(k) (û I(k) )| ≥ | dZ I (1) du I(1) (û I(1) )|. Hence, we obtain that there exist dither amplitudes for which γ = 0 and Using the sign function, we may rewrite this equation as follows: withg Z,i (t) in (59). For brevity, we omit the time argument in the following expressions. Using a first-order Taylor expansion of Z i (u i ), it follows that see (17). By combining (43), (59), (95), and (97), we get that the time derivative of the Lyapunov-function candidate in (96) can be written aṡ (99) Subsequently, applying Young's inequality yieldṡ Because the function Z i in (9) is twice continuously differentiable (see Assumption 5), and because its domain U i is compact, we have that d 2 Z i du 2 i is uniformly bounded. Using this and (45), we obtain that w i in (98) can be bounded by |w i | ≤ ν 2 c w (101) for some constant c w ∈ R >0 . By applying the comparison lemma [19,Lem. 3.4], and by using |u i (t)| ≤ κ 1 [see (44)] and the bounds in (82), (83), and (101), we obtain from (100) that V g (m i (t),g Z,i (t), q 2,i (t)) ≤ max ν 4 , sup τ ∈[0,t) |ỹ(τ )| 2 e − η 2 t V g (m i (0),g Z,i (0), q 2,i (0)) c V for all t ≥ 0, η ≤ ω 4 , κ 1 ≤ ην 5 , and some constant c V ∈ R >0 . From (53), (55), and (58), it follows thatỹ i (t) is uniformly bounded under the conditions of Theorem 6. Therefore, we obtain from (83), (96), and (102) thatm i (t) in (95) andg Z,i (t) in (59) are uniformly bounded. Moreover, the bound on the gradient estimation error in (60) follows from the same arguments.