Lightning optimizes: a threshold mechanism ensures minimum-path flow

A well-known property of linear resistive electrical networks is that the current distribution minimizes the total dissipated energy. When the circuit includes resistors with nonlinear monotonic characteristic, the current distribution minimizes in general a different functional. We show that, if the nonlinear characteristic is a threshold-like function and the energy generator is concentrated in a single point, as in the case of lightning or dielectric discharge, then the current flow is concentrated along a single path, which is a minimum path to the ground with respect to the threshold. We also propose a dynamic model that explains and qualitatively reproduces the lightning transient behaviour: initial generation of several plasma branches and subsequent dismissal of all branches but the one reaching the ground first, which is the optimal one.


Introduction
In lightning or gas electrical discharge, the current flow is essentially concentrated along a single path. Under very slow motion it can be seen that lightning starts by generating several branches and then develops by dismissing all of them but a single one, along which the energy is discharged [1]. This phenomenon has been deeply investigated. Several types of lightning are known, which are carefully described, e.g., in [2,3]. As far as the numerical modeling of the phenomenon is concerned, computational models for lightning simulation have been proposed in [4,5,6,7], while the fractal nature of lightning discharge has been investigated in [8,9,10,11]. Detailed surveys on the subject are also available; see, e.g., [12] and [13].
Here, we do not investigate the whole phenomenon in its complexity. We rather focus on a specific question about path formation in lightning discharge: we are interested in the initial phase of the process, when the lightning path is formed. Also, we consider the ideal case in which the lightning source is a single point and the final destination is a zero-potential ground.
This type of lightning, classified as Category 1 Lightning, "is the most common cloud-to-ground lightning. It accounts for over 90% of the worldwide cloud-to-ground flashes" [3]. Cloud-toground lightning begins with an initial breakdown and the consequent creation of a ionized channel, the stepped leader, which generates several branches. Once the stepped leader is close to the ground, it may be approached by channels originating from the ground, the connecting leaders. When the stepped leader finally connects the ground to the cloud, the return stroke is triggered, which is a ground-potential upward wave [2,14]. After the return stroke reaches back the cloud, the main branch reaching the ground is crossed by a long-duration discharging current: the continuing currents. In the meanwhile, the secondary branches originally established by the stepped leader are depleted: the continuing currents flow along the main path only [2].
Why is the lightning current eventually concentrated along a single path? Does this path enjoy any optimality property?
To mathematically address these questions, we consider an idealized model based on the assumption that lightning is mainly due to a dielectric breakdown of the air (gas in the case of discharges). The current-voltage diagram of a gas is characterized by two regions: for all the voltage values belonging to a symmetric interval around the origin, the current is very low (high resistivity); for voltage values outside this interval, the current becomes very large (low resistivity). The voltage value corresponding to the ends of the interval is called breakdown threshold. Then, we consider an ideal characteristic with conductivity approaching infinity when the electric field is larger than a threshold [15,5].
Lightning path can be interpreted as the solution of an optimization problem over a network. To formulate the problem, we consider a graph describing an electrical network, where capacitances and resistors with possibly nonlinear characteristic are associated with the links.
The grid model we use is akin to that proposed in [5,Eq. 5], which is the discretized version of nonlinear field equations [15,5]. We show that the steady-state solution minimizes a convex functional that, in the special case of linear resistors, turns out to be the total dissipated energy.
Conversely, if the resistor characteristic is a threshold-like function, the steady-state solution becomes the minimum path, where each current link is weighted by its local threshold voltage; hence, the optimal path is the "minimum-threshold path". Our main result is supported by a theorem and reinforced by simulations of randomly generated graphs with random threshold values, showing that the transient behavior of the model can faithfully reproduce the qualitative lightning evolution (see Materials and Methods).

Results
So, how does lightning optimize its path? To build a model, we consider the electrical grid network in Fig. 1, with a capacitive and a possibly nonlinear resistive effect between adjacent ground v=0 d i j u k Figure 1: The electrical grid model with capacitances and possibly nonlinear resistors connecting adjacent terminals. In the graph representation, each terminal is associated with a node and each electric component with a link.
terminals. Ground terminals are connected among them with zero resistance (ideal conductor ground) and the ground voltage is v = 0. The network is associated with a graph where the n+1 nodes represent the terminals and the links represent the electric impedances. In particular, node n corresponds to the zero-potential ground, while at the source node 0 a (current or voltage) generator is applied, with its other terminal grounded, inducing an input current d that enters the network.
As shown in Fig. 1, each link is assumed to be the parallel connection of a capacitor and a possibly nonlinear resistor. Injecting a current d in a node of the network leads, after a transient, to a steady state in which the currents flow only through the resistors. If these are linear, the steady-state solution corresponds to the current distribution that minimizes the total dissipated energy [16]: where R k is the resistance value associated with link k andū k is the steady-state current flowing through it. In this minimum energy configuration, steady-state currents are scattered all over the network (as in Fig. 2, left). In phenomena such as lightning, the situation is completely different [2,3]: after a transient, lightning "chooses" a single path (as in Fig. 2 Figure 2: The electrical current distribution at steady state: scattered in the case of linear resistances (left); concentrated along a single path in the case of threshold characteristics (right), when the current flows along the path that minimizes the sum of the dielectric rigidities of its links. Yellow (resp. blue) means presence (resp. absence) of flowing current. Assume that, denoting by v i and v j the potentials at the extreme nodes of link k, the resistance obeys the nonlinear law where the value of the threshold V k is the local dielectric rigidity of link k. Figure 3: Nonlinear current-voltage and voltage-current characteristics, where v denotes voltage and u current. A generic threshold-like current-voltage characteristic φ, with threshold V (red, left) and its inverse φ −1 (red, right). The ideal, sharp threshold characteristics (blue) can be seen as the limit of a sequence of sharper and sharper threshold-like functions φ and φ −1 .
If the nonlinear law approaches the ideal limit characteristic, depicted in Fig. 3 (left) along with its inverse characteristics (right), then we have the following results (derived in Materials and Methods).
• The steady-state current distribution minimizes the energy function J th = k V k |u k |.
• Considering the family P of all possible paths connecting the source node 0 to the ground node n, the whole injected current d flows along the path P h * ∈ P that minimizes the associated total energy: where k ∈ P h denotes that path P h crosses link k.
• In the transient, the injected current starts flowing along several "tentative" branches.
When one of these branches -corresponding to the minimum-threshold path -first connects to the ground (in general with the aid of a connecting leader), all the others branches are depleted, as shown in Fig. 4, obtained through simulations that faithfully reproduce this behavior.

Discussion
Therefore, the analysis of a grid circuit with capacitors and resistors having nonlinear characteristics unravels why flow phenomena such as lightning tend to concentrate the whole current flow along a single path, despite the availability of several admissible routes: this phenomenon is due to the threshold mechanism associated with the dielectric rigidity. In fact, for a nonlinear resistive network model, the solution of the flow equations minimizes a convex functional.
In the special case of linear resistance, this functional is the dissipated energy. In the case of ); the optimal path is crossed by the long-duration current (steady state). The color map goes from blue (no current) and light blue (low intensity current) to yellow (high intensity current).
threshold-like nonlinear characteristics, we have proven that the minimized functional is the sum of the currents along the links, weighted by the link threshold; hence, all the current eventually flows through the global minimum path if the links connecting the grid nodes are weighed by their local dielectric rigidity. In real situations, the dielectric rigidity can vary randomly and drastically in space, being a function of the local humidity, temperature, pressure and pollution.
This explains the seemingly random path of lightning: such a randomness is due to the gas current condition, because the lightning actually searches the right path.
The threshold model, including capacitive effects among nodes, faithfully describes also the transient and our simulations reproduce the behavior described in [2] and analyzed in [15,5].
Our model does not take into account inductive effects, considered by some authors in the return stroke analysis [2, pp. 169-170], but this does not invalidate our results because: (i) the minimum path analysis is carried out at steady state,v = 0, when the inductances are equivalent to shortcuts; (ii) the return stroke starts when the stepped leader has reached the ground, hence the path has been already "decided".
Our analysis reveals that lightning is one of the many phenomena in nature where a spontaneous optimization appears to take place [17,18,19], leading to the most efficient path choice [20]. In our lightning discharge model, the resulting steady-state current flow is globally optimal, even though the current flow is locally determined on the basis of the impedance characteristic of each single link.

Materials and Methods
The network is modeled as a grid graph with n+1 nodes and m links (more details on our model and our assumptions are in the SI Section 1.1). Node 0 represents the node where a current d is injected. The kth electric component is associated with link k connecting nodes i and j. Its impedance is given by the parallel connection of a capacitance and a possibly nonlinear resistor (cf. Fig. 1 in the main paper), so that the current u k flowing through the component can be written as where v i and v j are the voltages at the terminals i and j, while φ k (·) is the resistor currentvoltage characteristic function and C k is the capacitance.
Function φ k (·) is any symmetric increasing locally Lipschitz, or twice differentiable, function (the case in which φ k (·) is non-decreasing only is considered in the SI Section 1.1). In the case of a linear resistor, As mentioned in the main paper, we are interested in threshold-like characteristic functions whose value is close to zero in an interval and becomes very large if the voltage crosses the threshold value V k . In Fig. 3, the generic function φ is depicted (red curve, left) along with its inverse function φ −1 (red curve, right).
Following the description of Category 1 Lightning in [2,3], we distinguish two phases.
• First, the stepped leader "seeks the path to the ground": the current is relatively low and the capacitance effect dominates, leading to a fast variation of the voltages at the nodes, in the transient evolution of the model.
• Then, once the stepped leader has connected the cloud to the ground, dielectric breakdown is fully developed and the long-duration discharging current is triggered. We analyze this phase assuming steady-state conditions.
The initial branching phase is by far shorter than the second phase. Indeed, only by means of special very fast video equipments the first stage can be observed, while the second one can be captured by the human eye as we commonly experience. During the second phase, the current actually varies with time but its variation rate is extremely low with respect to the first phase.
In both phases, we show that the presence of a threshold mechanism is crucial to enable the observed behavior: it explains both the transient evolution of the phenomenon and the achievement of a minimum-path steady-state configuration.
We analyze by simulations the initial transient (first phase), during which the path is "decided" and the currents converge to a steady-state distribution.
First of all, however, we show that in the second phase, after an initial transient, the system converges to a steady-state, i.e., terminal voltages satisfy the conditionv i = 0. In this state, the non-null currents define a single flow along the minimum path in terms of dielectric rigidity.

Steady-state analysis: the chosen path is the minimum path
We start by showing that a threshold mechanism yields steady-state minimum-path flow.
We consider the functional where the index k refers to the links (see SI Section 1.1 for details). As a first result we have the following proposition, proven in the SI Section 1.3.
Proposition 1 Given the injected current d, at steady state (v i = 0) the current distribution in the network minimizes the functional given by (3). If the φ k are strictly increasing, the optimal current distribution is unique.
Note that dimensionally J(u k ) is a power, since φ −1 k (I) is a tension integrated with respect to a current I. Consistently, in the special case of linear resistances R k , namely when Here we are interested in the case in which the resistor characteristic is threshold-like. The ideal threshold function corresponding to a dielectric rigidity value V k is (see Fig. 3, blue) This curve is represented in blue in Fig. 3 (left). This ideal characteristics is physically unfeasible and will not be used for our simulations. However, the corresponding optimization problem is still well defined. Indeed, the inverse function of φ th k (the blue curve in Fig. 3, right) is and for this choice the functional in (3) becomes The following proposition holds and is proven in the SI Section 1.4.
Proposition 2 Given the injected current d, the admissible (compatible with Kirchhoff's laws) distribution of the steady-state link currentsū th k , k = 0, . . . , m − 1, which minimizes functional (6), corresponds to all the current flowing from the source node to the ground along a minimumthreshold path, namely a path P h * ∈ P (where P is the family of paths from the source node to the ground) that minimizes the cost which is the sum of all dielectric rigidities of the links along the path.
Functional (6) is not strictly convex, hence uniqueness is not ensured (see the SI Sections 1.1 and 1.4 for further details). However, the uniqueness assumption is generically satisfied; in fact, if the dielectric rigidities are randomly generated, the probability of finding two or more minimum paths with the same rigidity is zero, hence we can assume that the minimum path is unique.
The next step is to show that the closer a characteristic function (red) is to the ideal threshold φ th k (v) (blue) (Fig. 3), the closer the current distribution is to the minimum-path distribution.
Given a sequence of characteristics φ r k (v), r = 1, 2, . . . which "converge to the ideal one" and are physically feasible, so that the corresponding steady-state solutions are uniquely defined, these steady-state solutions converge to the minimum-path distribution. This property is formalized in the following theorem.
Theorem 1 Consider a sequence of characteristics φ r k , r = 1, 2, . . . , such that, for all r, the corresponding steady-state current distributionū r k in the links are uniquely defined. Assume that the threshold-like characteristics converge to the ideal one Moreover, assume that the minimum path in terms of sum of dielectric rigidities is unique. Then the link current distributionsū r k converge toū th k , u r k →ū th k as r → ∞ , namely, the current distributions converge to the one with the whole current d flowing along the minimum path.
Function φ th k is an idealized version of the gas dielectric characteristics in which the admittance is virtually zero for small voltage values and very large if the voltage is larger than the threshold value known as dielectric rigidity. In practice, true characteristics can be reasonably approximated [15,5] by a continuous curve that drastically increases after the threshold. The meaning of the theorem is that, if these characteristics are sharp and close to φ th k , then the current tends to flow along the minimum path. The result does not rely on any specific characteristic model: only the property of the characteristics becoming "close to the ideal" is essential.
In the model, we consider cell-to-cell capacities, but other capacities, such as capacities with respect to the ground, can be considered and the analysis remains valid, since at steady state the current through the capacities is zero. It is also fundamental to remark that the result is topology-independent: we could consider any network topology, not necessarily a grid. Also, we could consider conductive elements on the ground; in this case the lightning may find the minimum cost path as the one that connects the source to the grounded object. Some examples are in Fig. 5.

Transient analysis: seeking the minimum path
We show here that a threshold mechanism also explains the lighting transient behavior, which can be described by the dynamic model This system asymptotically converges to the steady-state conditionv(t) = 0, which leads to the condition Bφ B v(t) = d corresponding to the constraint of the optimization problem considered in the steady-state analysis. A detailed stability analysis is in the SI Section 1.2. The transient analysis has been carried out via simulation, using a standard ODE solver. In particular, to numerically demonstrate the dynamic behavior of the system, we have performed many simulations for different values of the dielectric rigidity. Videos are available to display some particularly significant cases (see the Supplementary Material Movies S1-S6 for details). 1 The characteristic function φ k (v k ) can be any locally Lipschitz or continuously differentiable function that is non-decreasing and has a very high slope after the threshold. For simulations purposes, we have adopted the piecewise-linear threshold-like functions (shown in Fig. 6, left), all with the same plasma conductivity r = 800, whereas V k has been randomly chosen for each link in the interval V k ∈ [0.5 − δ/2, 0.5 + δ/2] with uniform distribution. The variability of the dielectric rigidity is described by δ, while is a very small number representing the negligible conductivity under the threshold V k (we have set = 10 −5 ).
Other "sharp" characteristic functions would produce the same behavior. We also simulated the system with the polynomial φ k (v k ) = (v k /V k ) 2r+1 (shown in Fig. 6, right), which yields comparable results for large enough r, as expected. Yet, the non-Lipschitz nature of the polynomial function is numerically challenging and requires large computational times and specialized integration routines for stiff systems.
In all our numerical experiments, the capacitances have been taken all equal. We have set C k = 1 for all k, without restriction, since changing the capacitance value is equivalent to scaling time, hence the steady-state value is unaffected. Fig. 4 reports four instants of the simulation with δ = 0.7 (see Supplementary Material, Movie S5). It can be seen that, for larger variability of the dielectric rigidity, i.e. larger δ, the initial branching activity is more intense.
However, the asymptotic behavior is qualitatively the same regardless of the value of δ, with no exception: a single branch survives, which is numerically verified to be the minimum path in terms of total dielectric rigidity, as expected.

Decentralization, topology-independence and limitations
Remarkably, as we have stressed in the main paper, the steady-state current flow in our lightning discharge model is globally optimal, even though the current flow is locally determined on the basis of the impedance characteristic of each single link. In the context of distributed flow control in networks [21], this kind of mechanism is called network-decentralized [22,23,24], and localized strategies have been shown to lead to a globally optimal behavior [25,26,27,28,29,30]. In the considered network-decentralized control strategy, each links locally decides how much current flows through it. This approach is completely different from Dijkstra's decentralized minimum-path algorithm, which is based on decentralized dynamic programming techniques [31,32] and in which the routing decision is made at the nodes: each node locally decides to which of the outgoing links an incoming unit of flow must be redirected. In our setup, a "link decision" is made: each admittance can be interpreted as an agent that locally decides how much current is allowed to flow.
Our results are independent of the topology of the network, which does not necessarily need to be a grid graph with square cells: other topologies would lead to a minimum path solution.
However, we stress that our model is far from capturing all the complex aspects of lightning.
Its validity is limited to the beginning of the phenomenon, until the return stroke is triggered, because in this initial stage the path is chosen. After the stepped leader has reached the ground, the current follows the "chosen minimum path" until the end, as it is experimentally well documented (and confirmed by our simulations), because this path becomes a ionized channel with low resistance. So our model and simulations are not expected to be a faithful quantitative reproduction of the whole lightning process (including discharge endurance, power dissipation and so on), but their significance is limited to the first part. Yet, the qualitative behavior, with the discharge following the minimum path, has been always confirmed with no exception. Moreover, our model does not consider other aspects such as the ground currents, which are not relevant to the path choice. We have considered the so-called Category 1 Lightning, cloud-to-ground, which is the most common type of lightning; however, the model can be adapted to any type of lightning of gas discharge. n corresponds to the zero-potential ground, while at node 0 a (current or voltage) generator is applied, with its other terminal grounded, inducing an input current d that enters the network.
Consider the k-th electric component of the network, associated with the graph link k ∈ L, k = (i, j), which connects node i to node j: its impedance is given by the parallel connection of a capacitance and a possibly nonlinear resistor, so that the current u k flowing through the component can be written as corresponding to the admittance equation (2)), where v i and v j are the voltages at the terminals i and j, while φ k (·) is the resistor current-voltage characteristic function and C k is the capacitance.
We consider the generalized node-link incidence matrix of the graph G, which is the matrix B ∈ {−1, 0, 1} n×m obtained by assigning an arbitrary direction to each link k = (i, j) of G and setting a 1 entry in position i (source node) and a −1 entry in position j (destination node), and zero elsewhere, in the corresponding k-th column of B, and then removing the row corresponding to node n. In particular, the m columns of B are associated with the links representing the electric components and its n rows are associated with the nodes representing terminals. Links coming from the external environment (associated for instance with an injected current) have a single nonzero entry, equal to −1, corresponding to their destination node and links going to the external environment have a single nonzero entry, equal to 1, corresponding to their source node: in our model, the connections to the external environment are the (zerovoltage) ground and the source of supplied power. The following assumption holds, because we have considered the ground zero-potential node as an external node.
Assumption 1 The network graph G is connected internally and connected to the external environment. As a consequence, matrix B has full row rank.
Links are associated with the currents u 0 , . . . , u m−1 flowing through the individual electrical components, which we group in the vector u ∈ R m ; nodes are associated with terminal voltages v 0 , . . . , v n−1 , grouped in the vector v ∈ R n . Then, the admittance equation (2) can be written as while the current balance at the node i is where B i is the i-th row of B and d i is the i-th element of the vector of externally supplied current.
Merging equations (9) and (10) yields the dynamics of the overall circuit G, in terms of voltages and currents, which is described by the discretized space model [15,5] with equations where C = diag{C 0 , C 1 , . . . , C m−1 } is a diagonal matrix whose diagonal elements are the capacities C k , u = [u 0 , u 1 , . . . , u m−1 ] is the vector whose components are the currents along the links of G,d = [d, 0, . . . , 0] ∈ R n is the input current vector, v = [v 0 , v 1 , . . . , v n−1 ] is the vector whose components are the voltages at the nodes of G,v is the time derivative of vector v and φ(·) = [φ 0 (·), φ 1 (·), . . . , φ m−1 (·)] is the vector of the characteristic functions.
The following general assumption is considered.
Assumption 2 Each characteristic function φ k : R → R is a possibly nonlinear odd monotonically increasing function and locally Lipschitz.
This assumption can be weakened by requiring φ k to be monotonically non-decreasing only.
The shape of a generic characteristic function satisfying Assumption 2 is shown in red in Fig. 3. The assumption implies that, for each link k ∈ L, φ k is invertible. Then, for each k, the function where g k . = φ −1 k is the (monotonically increasing) inverse function of φ k , is well defined in (−∞, +∞).
Functions f k are continuously differentiable. In addition, they are strictly convex, since their derivative g k is an increasing function almost everywhere (f k = g k > 0 almost everywhere; in fact, g k may be not defined in some isolated points, e.g. in u = 0 for g k (u) = u 1/3 ).
If we assume that φ k is non-decreasing only, then we have convexity but not strict convexity of f k .

The dynamic model
We report here the stability analysis of the complete model (11), which we can rewrite in the equivalent formv The stability of this type of systems has been studied in the literature [33,28]. Consider the steady-state vectorv, such that and denote by x the shifted variable defined as x(t) = v(t) −v, whose time variation iṡ Since φ is a vector of strictly increasing functions, we can write where ∆(v) is a diagonal matrix of strictly positive continuous functions [26,28,29]. Hencė Consider the positive definite Lyapunov function candidate V (x) = 1 2 x BCB x, which is the energy stored in the capacitors. Its derivative is negative definite: This ensures asymptotic stability of the steady-state solution.

Proof of Proposition 1
We have to prove that the steady-state current distribution in the network induced by a constant current injectiond, achieved whenv = 0, namely when is indeed the current that solves the optimization problem where (15) is the flow constraint imposed by Kirchhoff's current law [16].
In view of the assumptions on φ k , the function f k defined in (12) is continuously differentiable with increasing derivative, hence strictly convex. Therefore, the optimization problem (14)-(15) is strictly convex and has a unique solution, achieved by applying the first order Karush-Kuhn-Tucker conditions to the Lagrangian function where λ ∈ R n is the vector of Lagrangian multipliers. The derivative with respect to u must be zero, hence we get Now, the first derivative of the elements of f is f k = g k , which is invertible with inverse φ k . As The solution of the optimization problem is therefore the unique solution of the system (13), Interestingly, the Lagrange multiplier vector is the steady-state voltage, In the case of non-decreasing functions φ k , we still have convexity but not strict convexity: the result holds, but the minimizing distribution may be non-unique.
Finally note that, in the special case of linear resistances R k , namely when Hence, in the nonlinear case, the minimized functional is not the dissipated energy as in the case of linear resistances.

Proof of Proposition 2
We have to show that the limit optimization problem associated with the limit characteristic φ th k (the limit of φ r k as r → ∞), admits as its optimal solution the current distribution with all current d channeled along the shortest path.
Assume that the injected current is positive: d > 0 (the case d < 0 is identical). Let u * denote the current distribution solving the optimization problem (18)- (19). To keep the proof simple we assume that all links in the network have been oriented in such a way that u k ≥ 0.
This is not a restriction since link orientation is arbitrary.

Now consider the modified problem
where the absolute value has been removed and a positivity constraint has been added.
The solution u * of the previous problem is a feasible solution of the new problem, because its elements are nonnegative by construction. It is also optimal for the new problem. Indeed, if another solutionũ were found with a lower cost, this would be a feasible solution also for the original problem (18)- (19) and would have a cost smaller than that of u * . Now note that, to solve (20)-(22), we can just take d = 1 and then scale the solution (by the true value d > 0). The proof is concluded by noticing that (20)- (22) with d = 1 gives the minimum cost path [16], with optimal cost d m−1 k=0 V k and the whole flow through the minimum cost path. An energetic interpretation is that the solution minimizes the overall power, measured as the product between the current flowing in a link and its dielectric rigidity.

Proof of Theorem 1
We have to prove that the steady-state solutions u r * associated with the characteristic functions φ r k , which have been shown to be the minimizers of converge to the solution of (18)- (19) if this is unique (equivalently, the minimum path is unique).
Denote by J r and J th the cost functionals of the considered optimization problems, Since g r k is strictly increasing, f r k is strictly convex and, in turn, also J r is strictly convex. Hence, the minimizer vector u r * of (14)-(15) is unique.
Let J th * = J th (u * ) be the optimal cost of the limit problem. For any y, f r k (y) −→ r→∞ V k |y|.
Then, J r (u) −→ r→∞ J th (u) and, in particular, J r (u * ) −→ r→∞ J th (u * ) = J th * . This means that the sequence of optimal costs {J r (u r * )} r∈N is upper bounded by a sequence {J r (u * )} r∈N that converges to J th * as J r (u r * ) ≤ J r (u * ) −→ r→∞ J th * .
Functionals J r , for all r, as well as J th , are radially unbounded because they are the sum of non-negative radially unbounded functions f r k . Then, their optimal solutions u * , respectively u r * are finite. In view of (25) there existsJ > 0 for which these solutions u r * are inside the compact set S r * = {u ∈ R m : J r * (u) ≤J} .
By construction J r (u) converges to J th (u) = V k |u k |, which is radially unbounded. Then, the sequence of sets S r * is uniformly bounded in a compact set S, hence all optimal solutions are uniformly bounded: {u r * } r∈N ∈ S.
We prove the convergence u r * −→ r→∞ u * by contradiction. We assume that u r * −→ r→∞ u * .
Negating convergence to u * implies that there exist an open neighborhood, U ⊂ S, of u * and a sub-sequence of u r * that is in the complement of U in S, namely in the compact set S \ U. In turn, this sub-sequence confined in the compact set admits a sub-sub-sequence that converges to some point u • ∈ S \ U. Hence there exist a sub-sequence {ur * }r ∈N of the original {u r * } r∈N that converges to some vector u • = u * , being N = {N 1 , N 2 , . . . } an infinite ordered set of increasing integers. All vectors ur * satisfy the constraint Bur * =d as they are solutions to problem (14)- (15).
Hence, also the limit vector u • does satisfy Bu • =d. Then, the proof can be concluded by showing that which is a contradiction, because it would imply that either J th * is not the optimal as assumed, or (if equality holds) that the optimization problem (18) has two minimum points, u * and u • , against the uniqueness assumption.
To prove (26), the first step is to note that there exists a finite value g such that, for all k and for all sufficiently large r, g r k (u k ) ≤ g for all u k such that u ∈ S, since g r k (u k ) −→ r→∞ V k and S is a compact set, and hence is bounded. As a consequence, for a sufficiently large r, J r has a uniformly bounded gradient since, for all y, |∂J r /∂u k | (y) = g r k (y) ≤ max u k :u∈S {g r k (u k )} ≤ g .
Hence, as J r is convex for all r, there exists a constant L such that, for sufficiently large r and for all p, q ∈ S, |J r (p) − J r (q)| ≤ L p − q .
Then, forr ∈ N sufficiently large, the following inequalities hold: where the last inequality and the limit come from (25) and the fact that {ur * }r ∈N is a subsequence of {u r * } r∈N with limit u o . Then we have shown the contradiction (26), which concludes the proof.