I have been meaning to play around with this algorithm for a very long time, and with a few days off work it seemed like the ideal time to get stuck in. Climate change is a hot topic, with a lot of research into the health and economic impacts of climate change. Physiological effects of climate change are modeled using heat stress indicators, such as the **wet bulb globe temperature** (WBGT), and there is a lot of interest in the development of accurate predictive models for heat stress indicators. Lemke and Kjellstrom (2012) published a comparative study on the most widely used models for WBGT, finding that Bernard’s model was the most suitable for their use in subsequent research and publications. (1) It is the computational algorithm behind Bernard’s model that we explore in this article.

The vast majority of computational algorithms in science are iterative. While there may be numerical solutions, there is not always a numerical solution and often numerical solutions do not scale well (for example, the normal equation for linear regression). Iterative models provide a way to model increasingly complex problems over very large data sets. Euler’s method is a widely used iterative algorithm which continues to update a hypothesis, h(x), while some boundary condition holds. The general form of Euler’s method is shown below:

repeat { h(x) := h(x) + alpha.h(x) } where alpha << 1

Given the general form above, where alpha << 1, then the *“future hypothesis”* is* “mostly what it is now, plus a little bit”*. Generally, iteration will continue until some boundary condition is met. For example in linear regression, the hypothesis is updated until the mean squared error is minimised. Often the choice of ** alpha **is the single most critical element to the efficiency of the algorithm, where if alpha is too large the hypothesis may oscillate wildly and never converge, or if alpha is too small the rate of convergence will be slow.

**The Relaxation of Tw: The Existing Algorithm**

The computational modeling of WBGT is a practical example of Euler’s method in action. The critical step is the *relaxation of Tw*, an iterative process strongly reliant on the input values. Currently, there is no accurate numerical solution to calculate Tw and the physical measurement of atmospheric Tw is inherently inaccurate.The algorithm, using R, is given below:

esat <- function (t) { return (0.6106 * exp(17.27 * t / (237.3 + t))) } mcphersons <- function (ta, td, tw) { return ( (esat(td) - esat(tw)) * (1556 - 1.484 * tw) + 101 * (ta - tw) ) } relax <- function (ta, td) { # two additional variables to record intermediate state for visualisation twX <- c(td) mcpX <- c(10000) # initial conditions mcp_i <- 10000 mcp_j <- 10000 tw <- td # Increment Tw while these boundary conditions hold while (mcp_i > 0 & mcp_j > 0 & tw < ta) { mcp_j <- mcp_i mcp_i <- mcphersons(ta, td, tw) tw <- tw + 0.3 twX <- c(twX, tw) mcpX <- c(mcpX, mcp_i) } return (list(tw=twX, m=mcpX)) }

In the algorithm above, Tw is slowly “relaxed” until McPhersons formula reaches zero. With each iteration, Tw is advanced by 0.3 degrees and McPhersons formula is re-evaluated. Figure 1 shows the values of McPhersons Formula at each iteration:

The plot above was recorded with inputs (ta, td) = (20, -8). While it appears the the relaxation is linear within this range, as the number of iterations increase the relaxation is noticeably non-linear. Relaxation continues until the value of McPhersons formula crosses y = 0.

The number of steps required to fully relax Tw is strongly dependent on the input values Ta & Td and can take several thousand iterations with real-world data. When applied to large data sets (for example the CRU data sets (2) for which there is currently more than 30 years worth of measured temperature observations) the relaxation can take many hours to complete.

In previous experiments we have successfully been able to short-circuit this iteration by projecting a linear trend and then stepping backwards to the target value. However, there are a number of issues with this approach. Firstly, it adds complexity to the code to be able to determine the linear model, extrapolate out to y = 0 and then step backwards. Secondly, while it is effective at reducing the total number of iterations and improving the overall performance, it departs from the classic iterative model which makes it harder to understand and share with other researchers. Below, we explore a modified iterative scheme using ** gradient descent** which results in a faster rate-of-convergence than the existing method and has the advantage that gradient descent is a well understood and widely-used algorithm.

**Exploring the Relaxation of Tw**

We know, from observation of the relaxation profile in Figure 1, that the general trend in relaxation is downwards. Looking at the code, we can see that the relaxation of Tw depends on the value of McPhersons formula and that the boundary condition is approximately where McPhersons formula dips below the x-axis. Given this boundary condition, the target value of Tw is approximate to within +/- 0.3 degrees. However, if we plot the square of McPhersons formula we observe a more accurate limit:

Figure 2 provides a different perspective of the relaxation of Tw. The plot of the square of Mcphersons formula is a classic convex function with a global minimum that corresponds to our target value of Tw. Looking at Figure 2, it seems that we could employ *gradient descent* to relax Tw. An algorithm for gradient descent is given below:

descend <- function (ta, td, alpha=0.001, eps=0.05) { # initial conditions twX <- c(td) tw_i <- td tw_j <- 100000 # this is just an impractical temperature # increment Tw, while delta(tw) > epsilon while (abs(tw_i - tw_j) > eps & tw_i < ta) { tw_j <- tw_i tw_i <- tw_i + alpha*mcphersons(ta, td, tw_i) twX <- c(twX, tw_i) } return (twX) }

And if we plot Tw, we can see that the *descend()* function asymptotes at the target value of Tw, ~ 8.2.

Figure 3 shows the relaxation of Tw via the gradient descent algorithm, which is able to accurately and efficiently approximate Tw. By using gradient descent, the accuracy of the approximation is improved (accurate to within the threshold *epsilon*) and the total number of iterations has been greatly reduced. The model can be further refined by adjusting the learning rate, *alpha*, and the threshold, *epsilon.*

**Comparing the two models**

To assess the suitability of the gradient descent algorithm, we compared the resulting values of Tw and WBGT to the original algorithm on a randomly generated test set with 10,000 observations. The experimental setup is shown below:

n = 10000 data <- data.frame( ta = sample(5:40, n, replace=TRUE) , td = sample(-20:20, n, replace=TRUE) , tw_relax = rep(NA, n) , tw_descent = rep(NA, n) ) for (i in 1:n) { original <- relax(data$ta[i], data$td[i])$tw descent <- descend(data$ta[i], data$td[i], eps=0.05) data$tw_relax[i] <- max(original) data$tw_descent[i] <- max(descent) } data$delta_tw <- data$tw_relax - data$tw_descent

The agreement between these two models was determined as the difference in Tw by the original *relax()* *algorithm* and the new *descend() algorithm* and is shown in Figure 4:

Figure 4 shows a reasonable agreement between the original algorithm and the gradient descent algorithm. It seems that the original algorithm has a tendency to over-shoot the true minimum of McPhersons formula and typically produces values of Tw that are 0.5 degrees greater than the gradient descent algorithm.

WBGT is calculated as:

We compared the values of WBGT using Tw determined by the original algorithm with the same for gradient descent:

Figure 5 shows good agreement of WBGT calculated by the original algorithm and the gradient descent algorithm. A 95 % confidence interval was determined for the agreement of WBGT(relax) and WBGT(descend):

Based on the 95 % confidence interval above, the gradient descent algorithm is entirely suitable for the approximation of WBGT and could serve as an efficient and reasonable substitute for the existing algorithm.

**Conclusion**

Accurate and efficient prediction of heat stress indicators, such as WBGT, is important for the modeling of climate change. There are no accurate numerical solutions for the calculation of WBGT, where Tw is unknown, which has led to the development of iterative algorithms. However, iterative algorithms are expensive when applied to large databases of real-world climate data. To improve the accuracy and efficiency of calculating WBGT we propose a method for the relaxation of Tw that uses the well-known and widely used gradient descent algorithm. We have shown excellent agreement between the original algorithm and the gradient descent algorithm to within 0.18 +/- 0.0015 degrees Celcius with a 95 % confidence interval. Overall, it appears that the gradient descent algorithm is able to accurately and efficiently approximate WBGT. Also, as gradient descent is widely used in learning algorithms, such as linear regression and logistic regression, the algorithm has the advantage that is it relatively easy to implement and share across research groups.

**References**

(1) Lemke, B., & Kjellstrom, T. (2012). Calculating workplace WBGT from meteorological data: a tool for climate change assessment. *Industrial health*,*50*(4), 267-278.

(2) The University of East Anglia, Climatic Research Unit. http://www.cru.uea.ac.uk/data

Pingback: Visual Learning & Discovery | :: NickBurns

Pingback: Visual Discovery | adatastory