24.2.1 Step 1: Replicate the Apparent Best Player

Leapfrogging (LF) is a multiplayer, direct‐search approach. Trial solutions (players) are randomly initialized through decision variable (DV) space. One player will have the worst objective function (OF) value, while another, the best. The worst player leaps over the best and “lands” at a random spot in the reflected DV space on the other side of the best. If the new leap‐over position is infeasible, then the same player leaps from the infeasible spot back over the best player and into a new random location within the new reflected DV space. Normally, for deterministic functions, whether the leap‐to location is better or worse, when the new location is feasible, the move is complete. Normally, an iteration is defined as the number of leap‐overs the same as the DV dimension. This means that if leap‐to positions are feasible, there are N function evaluations per iteration.

However, that definition of an iteration is arbitrary. An iteration in the CHD method has one function evaluation per dimension, but HJ has 1 + 1.5N function evaluations per iteration, ISD (with central difference) has 1 + 2N, and CSLS and SQ and Newton types have many more.

If a leap‐to position is up on a cliff, or otherwise worse than the leap from position, then that player will be the next to leap. As will be seen, if it jumps on and off of a cliff, this will make the local variation seem greater than the stochastic properties at the local optimum and confound the steady‐state identification. Accordingly, for stochastic functions if the leap‐to position is either worse or infeasible, then immediately leap it back. If it jumped up on a cliff, then the leap back will be an improvement. Call an iteration when there are N leap‐overs with improvement.

In LF, the player with the worst OF value leaps over the player with the best, but the best player might have had a fortuitous realization. Rather than truly being best, it might be a phantom appearance. One would not want that to lead the optimization to move more players into an undesired DV region. To accommodate the stochastic OF attribute, here, the best player is reevaluated at each leap‐over, and the worst of replicate values is taken as its representative OF value.

How many replicates? If you want to be c confident that the number of replicates will have found one of the worst fraction, f, of OF realization values, then

(24.4)

For example, if you desire to be 95% confident that there are enough replicates to find 1 of the 75 percentile worst values (1 in the 25% worst fraction), then

(24.5)

This would indicate the current OF value and 9 additional samples. However, if the best remains the best, 9 additional samples at the next leap‐over would represent a total of 19, and 9 more on the next leap‐over would represent 28. The worst of 28 is roughly the 90 percentile worst. I find that 1–3 replicates per leap‐over is adequate. If the player is a phantom best, then after a few leap‐overs, it will have enough replicates to acquire a value representing a local worst outcome.

The greater the number of replicates, the lower the chance that a converged DV* value will represent a bad solution masquerading as a good one. Choose c and f values to represent the situation. If it is a normal business economic analysis, then could be an appropriate number; however, where there is less tolerance for risk, perhaps is appropriate.

If reevaluation of what appears to have been the best player reveals that it is not the best, then the second best becomes the lead player and basis for the leap‐over, redirecting the search. This prevents end‐game convergence of the cluster in DV space and motivates the novel approach demonstrated here—stop when there is no improvement in the stochastic response with respect to iteration. Stop at steady state.

24.2.2 Step 2: Steady‐State Detection

Normally, in optimization of deterministic functions, the best trial solution will progressively improve with each iteration, asymptotically approaching its optimum. However, with replicating stochastic functions, at the optimum, the value of the best player will generate what appears to be a noisy value at a steady‐state average. As the optimizer moves toward the optimum on a stochastic surface, this noisy value will approach its noisy steady state (SS).

If one player is initialized in the vicinity of the global, then it may remain in the best spot, while others are gathering, and observing the OF response of the best player might identify SS prior to all players converging. However, if the worst OF is observed, then when it settles to SS, all players will be in the vicinity of the optimum. Accordingly the method is to observe the OF of the worst player and claim SS when it settles to a noisy SS w.r.t. iteration.

Appendix D provides a method to detect steady state in a signal. There are many methods to determine steady state (SS) in a signal. The method demonstrated here is based on a probable steady‐state detection algorithm, chosen for simplicity and effectiveness. The method first calculates a filtered value, X_f, of the process measurements, and then the variance in the data is measured by two methods. Deviation d₂ is the difference between measurement and the filtered trend, while deviation d₁ is the difference between sequential data measurements.

If the process is at SS, X_f is almost the middle of the data. Then a process variance, σ², estimated by d₂ will ideally be equal to σ² estimated by d₁, and the ratio of the variances, , will be approximately equal to unity, . Alternately, if the process is in a transient state, then X_f lags behind the data trend, the variance as measured by d₂ will be much larger than the variance as estimated by d₁, and ratio will be much greater than unity, . Being a ratio of variances, the statistic is scaled by the inherent noise level in the data and is dimensionless. Critical values for the ratio to confidently trigger the process being at steady state with independent and identically distributed variation (white noise) were developed. Filter coefficient values of 0.1 and an r‐critical value of 0.85 provide agreement with human perception to claim SS.

The method presumes no autocorrelation in the “time series” (iteration sequence) of the data. This is ensured by observing the worst player position after each iteration, not the best. The best player makes sequential moves toward the optimum or remains stationary for several iterations, creating autocorrelation in the sequential OF values. By contrast, the worst player is from an iteration‐to‐iteration, new randomized surface sampling. As the players converge on the optimum, the worst player moves from the high OF values to the low ones. With iteration, the worst OF value progressively drops until, at the minimum, its stationary value represents independent stochastic OF values from the proximity of the optimum.

24.4.1 GMC Control of Hot and Cold Mixing

Function 36 in the 2‐D examples code presents this application. Hot and cold fluids are mixed in line, and the objective is to determine the controller outputs (ideally these are the control valve positions), o₁ and o₂, to produce the desired mixed temperature and total flow rate.

With no uncertainty on model values, the contour of the 2‐D search for o₁ and o₂ appears as Figure 24.1a. However, with a relative uncertainty on each model coefficient (Gaussian with a 5% standard deviation), in Figure 24.1b the contours of one realization are obviously irregular, and each realization will generate a different irregular set of contours. The corresponding 3‐D plot of OF w.r.t. DVs (Figure 24.1c) reveals one realization of the irregular surface, and each realization will distribute the spikes and wrinkles differently. In general, Figure 24.1c shows that the minimum is somewhere along a curved path from the front right to the back left.

24.4.2 MBC of Hot and Cold Mixing

This is similar to the aforementioned physical process; however, it is a single controlled variable application. The objective is to determine the optimum sequence of future control actions to make the temperature best match a linear reference trajectory to the set point. Control is subject to rate of change constraints.

24.4.3 Batch Reaction Management

This simulates an algae farm and offered as Function 63 in the 2‐D Optimization Examples program. Algae is seeded in a wastewater pond, permitted to grow in response to sunlight and nutrients, and after a time the water is processed to harvest the algae. Profit depends on the mass of algae harvested. The longer the growth period, the more algae in the pond; but the mass asymptotically approaches a maximum as limits of nutrients make the algal death rate approach growth rate. Waiting a long time to harvest 100% of what is possible reduces the number of batches that can be generated per year. It might be better to wait until the batch is only 80% developed and to be able to harvest twice as many.

Additionally, the deeper the pond, the more volume there is to grow algae and, consequently, mass that can be harvested. But as the pond gets deeper, less sunlight and CO₂ will diffuse to the lower levels. Too deep a pond leads to the cost of processing excess water with little extra product to be harvested.

The objective is to determine the batch growth time and pond depth to maximize annual profit. Figure 24.2 illustrates one realization of the surface. Contrasting Figure 24.1c, in which the smooth area is in the vicinity of the optimum, in Figure 24.2 the fluctuations are most severe in the region of the optimum.

24.4.4 Reservoir and Stochastic Boot Print

These two simulators represent a water reservoir design problem and are Functions 18 and 20 in the 2‐D Optimization Examples program. The objectives of a reservoir are to trap excessive rain water to prevent downstream floods, to release water downstream to compensate for upstream droughts, and to provide water for human recreational and security needs. The design objectives also include minimizing the cost of the dam and land. The taller it is, the more it costs to build; but with its greater capacity, the probability of flood or drought impact will be lower, and the recreational and security features will be better. The fuller it is kept, the less it can absorb floods, but the drought, security, or recreational performance will be better. On the other hand, if kept nearly empty, it can mitigate any flood, but cannot provide recreation or drought protection. The OF represents the reservoir cost plus the economic risk (probability of an event times the cost of the event). The lower left axis in Figure 24.3 (one realization of Function 20) represents the dam height, while the lower right axis is the set point portion of full.

The OF value in Figure 24.3 depends on a probability of the flood or drought event, and each realization of the contour will yield a slightly different appearance. Note that starting in the middle of the DV space on the planar portion, a downhill optimizer will progressively move toward the far side of the illustration to a less costly reservoir, but into the region with the spikes. In that region for a small reservoir, or one that is kept too full or too empty, there is a probability of encountering a costly flood or drought event that the reservoir cannot mitigate. Moving in the downhill direction the optimizer may, or may not, encounter a costly event in the Monte Carlo simulation. If it does not, it continues to place trial solutions into a region with a high probability of a disastrous event and continues into the bad region as the vagaries of probability generate fortuitous appearing OF values.

Stochastic boot print, Function 18, is a surrogate for the reservoir, a simplified set of equations that computes faster. It’s response is similar.

24.4.5 Optimization Results

Figure 24.4a shows the results of a Levenberg–Marquardt search on the GMC control problem with convergence based on sequential changes in the DV, a classical measure. Trial solutions were initialized at random locations within the DV space, 100 times, and converged at 100 wrong places. The pattern is independent of the threshold on the change in DV or whether convergence is based on a change in OF value. The results are similar when other single‐point algorithms are used, even direct‐search optimizers such as Hooke–Jeeves. These algorithms get stuck on a local fortuitous good spot and cannot find a way out of the apparent, phantom minimum. In contrast, Figure 24.4b reveals the results of a particle swarm optimizer search starting with 20 particles, with an end‐game perching attribute to converge. The broad surface exploration and particle attraction toward the best leads the swarm to the vicinity of the global minimum. But, in that vicinity, the swarm converges on a fortuitous best, resulting in a range of phantom solutions.

The range of converged solutions is relatively insensitive to several orders of magnitude on the convergence criterion precision. In all figures the threshold value is 0.01 as an rms change on the DVs; however, a four‐order‐of‐magnitude range of 0.001–1 was explored. By contrast, the number of function evaluations to converge is strongly dependent on the convergence threshold. This is an issue for a user, who must specify the threshold values for the convergence criterion, but cannot have rational, a priori criteria for relating criterion thresholds to results. In the 100 trials, PSO required 2180 ANOFE (average number of function evaluations per trial) to converge.

Figure 24.4c reveals similar results for LF. For the same convergence criterion, the ending range of solutions is nearly identical to PSO, but LF averaged 1150 ANOFE, therefore the computational work is significantly lower.

In contrast, Figure 24.5 illustrates the impact of the method presented here. It is the same problem as that in Figure 24.4c; however, the OF value of the best player is based on the worst of the prior and one replicate value, and convergence is based on the automated, scale‐independent recognition of steady state. The replicates reduce the impact of fortuitous good spots from continuing to act as a draw, which significantly reduces the range of solutions nearer to the global precision. One might expect that replicates double the optimizer effort, but because of the SS convergence criterion, evaluations stop when continued work is substantially ineffective. The ANOFE for this set of 100 trials is 1090 (about 10% fewer), and the converged solutions are visibly more precise.

These illustrations motivate a choice of dual criteria for evaluating the method: (i) the range of converged solutions indicating precision and (ii) ANOFE indicating data generation cost.

Table 24.1 summarizes findings for each of the stochastic functions described previously, comparing LF with a conventional convergence criterion of DV rms less than 0.01 and LF with a replicate evaluation and the SS detection on the worst player value. The table entries are the precision of the cluster of converged solutions (as indicated by rms value of the scaled DV range) and the work required to obtain the solution (as measured by the average number of function evaluations to converge). The values are based on 100 trials, rounded to reflect run‐to‐run variability. For similar converged cluster size (DV rms), the proposed convergence criteria finds the optimum with nearly half the ANOFE. Or for similar ANOFE the proposed method provides greater precision on the DV optimum.