A heuristic for improving clustering in biomass supply chains

Abstract

Clustering is commonly used in various fields such as statistics, geospatial analysis, and machine learning. In supply chain modelling, clustering is applied when the number of potential origins and/or destinations exceeds the solvable problem size. Related methods allow the reduction of the models’ dimensionality, hence facilitating their solution in acceptable timeframes for business applications. The weighted minimum sum-of-square distances clustering problem (Weighted MSSC) is a typical problem encountered in many biomass supply chain management applications, where large numbers of fields exist. This task is usually approached with the weighted K-means heuristic algorithm. This study proposes a novel, more efficient algorithm for solving the occurring weighted sum-of-squared distances minimisation problem in a 2-dimensional Euclidean surface. The problem is formulated as a set-partitioning problem, and a column-generation inspired approach is applied, finding better solutions than the ones obtained from the weighted version of the K-means heuristic. Results from both benchmark datasets and a biomass supply chain case show that even for large values of K, the proposed approach consistently finds better solutions than the best solutions found by other heuristic algorithms. Ultimately, this study can contribute to more efficient clustering, which can lead to more realistic outcomes in supply chain optimisation.

KEYWORDS:

• Clustering
• weighted MSSC
• supply chain
• biomass
• heuristics

1. Introduction

Biomass is established as an important alternative to fossil fuels in the quest for renewable energy sources. Biomass plays a versatile role in energy systems, contributing not only to power and heat generation (Tziolas et al., 2017) but also serving as a source for biofuels with widespread applications, particularly in the transportation sector. Currently, biomass accounts for over 70% of the renewable energy production globally, while it contributes to a 4% on the fuel demand in road transport. According to the European Union's (EU) long-term strategy, the objective is to achieve a climate-neutral EU by 2050 (European Commission, 2020). To this end, sustainable biofuels have the potential to play a decisive role in decarbonising the transport sector (road, aviation, maritime, etc.), while ensuring a high level of affordable and accessible transport connectivity (European Parliament, 2021). The emphasis is especially on second-generation biofuels (known as advanced biofuels), which originate from dispersed sources, such as agricultural residues, biomass fractions from mixed municipal and industrial wastes, animal manure, wastes and residues from forestry, etc. According to the EU Renewable Energy Directive 2018/2001, advanced biofuels and biogas production are expected to increase at least to 1% of the final energy consumption in the transport sector in 2025 and at least to 3.5% in 2030 (IEA, 2017).

However, the economic viability of using biomass as a renewable energy source is affected by various factors, like uncertainty in feedstock availability, seasonal fluctuations in biomass production, and widespread geographical dispersion (Hassan et al., 2019; Psathas et al., 2022). Considering also the volatile and competitive global energy market, the continuous technological developments and the more stringent standards and commitments on promoting Sustainable Development Goals, the design of robust biomass supply chains becomes a highly demanding procedure needing a comprehensive system analysis. The overall challenge is to develop a decision-making framework able to amalgamate a multiplicity of factors, manage uncertainty and provide valuable key performance indicators (Abdussalam et al., 2021).

A commonly applied approach to support decision-making in such complex supply chain schemes is optimisation. Optimisation is tightly coupled with studying, modelling and solving a multiplicity of industrial and energy problems at any functional level: short-, medium- and long-term planning (Grossmann, 2012). Supply chain optimisation models help the decision makers and interested stakeholders make the right strategic and tactical decisions across the upstream, midstream and downstream stages, in order to generate a considerable profit and preserve overall efficiency and sustainability (Abdussalam et al., 2021). Mathematical Programming provides the framework to construct such optimisation models able to include a variety of feedstock types and sources, intermediate and final products, energy forms and product streams, technologies, and processes, while incorporating economic, environmental, and social aspects related to a typical biomass supply chain. These models can also consider various cost categories (harvesting, processing, transportation, etc.) and special constraints while they can express the necessary conditions such as biomass availability and geographical dispersion. Furthermore, these models can have interoperability with Geographical Information Systems (GIS) which provide rich spatial data (e.g. physical locations of fields, warehouses, conversion facilities, roads, etc.) and are equipped with analytical capabilities (e.g. K-means clustering) (Lozano-García et al., 2020).

The study and analysis of real-life biomass supply chains often includes large numbers of biomass fields and entails the expression of many series and types of variables and mathematical relationships, leading to the construction of detailed and large-scale models. When computational complexity increases, the commonly used exact solution algorithms, can become insufficient to obtain solutions in reasonable timeframes (Martín et al., 2023). In such cases, a reduction of the size of the optimisation problem, achieved by the reduction of the number of decision variables is required. This can be achieved by employing clustering and aggregation techniques, to formulate manageable and solvable biomass supply chain models (Potrč et al., 2020).

Clustering is a technique commonly used in various fields beyond the biomass supply chain field, such as statistics, geospatial analysis, and machine learning, being tightly linked with the optimisation concept. Specifically in supply chain modelling, clustering is naturally required when the number of potential origins and/or destinations is quite large and the native routing assignment significantly enhances the overall combinatorial optimisation load and the solvability of the cost/distance minimisation problem (Derya et al., 2023). Transforming the initial depiction of supply chain nodes to an equivalent partitioned representation can expand the applicability and suitability spectrum in supply chain management gaining benefits in model generation and solution times.

The motivation for dealing with the clustering problem in biomass supply chains originated from the research and modelling objectives of the Horizon 2020 CERESiS project, (CERESiS, 2023), which aims to produce sustainable liquid biofuels through the cultivation of suitable energy crops for the remediation and management of contaminated land. The optimisation of the respective biomass supply chain is significantly affected by the inherent geographical dispersion of the biomass sources and the fluctuation of the feedstock supply. Common clustering approaches in the literature either employ the standard K-means algorithm or spatially divide the map into a grid, arbitrarily placing the coordinates of the aggregated field in the centre. Both approaches overlook the total hectarage of each field, potentially leading to significantly inaccurate transportation cost estimates, hence affecting the accuracy of the whole supply chain modelling. This happens because these algorithms assume equal-sized clusters, while in real-world applications, biomass fields often have varying hectarage. It should be noted that similar issues exist in many other supply chain optimisation contexts where a large number of upstream supply nodes/origins exist, such as recycling, waste collection, agri-food, etc.

To address this issue, the present study introduces a novel weighted K-means algorithm accounting for both field coordinates as well as field size, resulting in a more accurate representation of transportation costs. Incorporating this innovative approach, the aim is to render large biomass supply chain models computationally feasible while also mitigating any inaccuracy linked to traditional and arbitrary clustering. The overall modelling improvement facilitates the end goal of drawing robust and economically viable biofuel or bioenergy production plans based on energy crops cultivated in geographically dispersed land parcels.

The methodological novelty lies in the development of a new heuristic approach to solve the weighted sum-of-squared distances minimisation problem, that considers field size in the clustering process, thus leading to more accurate clustering results in respect to the real system compared to the weighted version of the well-known K-means heuristic, when the size of each field/origin is important. The application novelty lies in applying a heuristics approach for clustering that can consider the field size specifically in the biomass supply chain context. The proposed approach can be equally applied in other supply chain contexts, where the weight of each origin could denote quantities of materials available or other values, instead of the field size of the biomass supply chain problem.

This paper is structured as follows: Section 2 presents the literature review. Section 3 introduces the model formulation. Section 4 describes the column generation approach. Section 5 analyses the results of the computation experiments. Finally, section 6 concludes the study and discusses future research perspectives.

2. Literature review

Numerous studies deal with the optimisation of biomass supply chains considering various logistics and processing stages involved, from the energy crops cultivation or feedstock supply up to biofuel or bioenergy production. For example, Fattahi and Govindan (2018) presented a multi-stage stochastic programme that minimises the cost of a biofuel supply chain to satisfy the demand, accounting for seasonal and stochastic biomass availability, and incorporating environmental and social impact. In a study by Momenitabar et al. (2023) the biofuel supply chain was optimised with a focus on sustainability criteria to determine the optimal siting of the conversion facilities. Sosa et al. (2015) developed a linear programming tool taking into account moisture in biomass, showing its impact on harvest area timeline and storage, as well as transportation elements within the supply chain, including strategies to optimise truckloads within legal weight and volume constraints. In the abovementioned studies as well as in much of the relevant literature, a detailed record of biomass/feedstock field data is often lacking or unavailable due to the inherent difficulty of gathering input from dispersed biomass fields. This is the reason that detailed data from biomass supply hubs’ location (e.g. fields) is usually consolidated or assumed to be available in a single geographical point, hence limiting the accuracy of the modelling process outcomes. However, in instances where this abundance of information is available and is used as the starting point for the biomass supply chain modelling, a method of field aggregation may become necessary for grouping the fields, to reduce the overall model complexity and allow its solution in reasonable times using mathematical optimisation.

For this purpose, clustering algorithms are important and inextricably linked with biomass supply chain design, management and optimisation. A common approach in the biomass supply chain literature lies in treating the substantial volume of data from supply hubs, which renders the supply chain optimisation problem complex and in certain cases unsolvable. For example, O’Neill et al. (2022) developed a stochastic mixed-integer linear programming (MILP) biofuel supply chain model, aiming simultaneously at the maximisation of economic performance and the minimisation of the environmental impact. Their dataset contained over 40,000 fields, which rendered the problem unsolvable with exact optimisation methods. To reduce computational cost, clustering was employed at resolutions of 25, 100, 225, 400, 625 km² and county level. Uen and Rodríguez (2023) presented a MILP model maximising the net present value of a food waste supply chain for renewable energy generation. The model involved a large dataset with 102 food waste sources and 621 candidate conversion facilities sites (anaerobic digestion and co-digestion), resulting in high computational cost. To address this issue, they used the K-means algorithm to create more representative centroids based on spatial proximity; the final cluster selection consisted of 9 clusters for 291 aerobic digestion candidate sites and 12 clusters for the 330 co-digestion candidates. Psathas et al. (2022) developed a biomass-to-biofuel supply chain that integrated large-scale centralised, smaller-scale decentralised as well as mobile processing facilities. The complex MILP model, applied to a use case involving over 3,000 biomass fields, showcased the inability to achieve a solution for such a high-dimension optimisation problem, prompting the need for clustering the fields in a more manageable group of 16 clusters using the K-means algorithm.

The application of clustering algorithms in the context of biomass and bioenergy is also apparent in Vehicle Routing Problems (VRPs). Ayoub et al. (2007) developed a bioenergy decision support system which combined optimisation and simulation for the design of bioenergy production. In their system, they employed fuzzy C-means clustering to define the location and optimal size of storage and bioenergy conversion facilities based on biomass collection points, with the end goal being the minimisation of transportation costs, CO₂ emissions and number of workers. Zamar et al. (2017) applied near-neighborhood search and a modified K-means clustering technique for the route optimisation in the bale collection problem, aiming to minimise travel time and fuel consumption. The constraint of this new K-means algorithm limited the temporary storage sites to valid locations. Zamar et al. (2017) presented a stochastic VRP model for the biomass residue collection from a set of sawmills to identify the optimal routing schedule. Due to variations in biomass quality, the aim was to maximise the ratio of energy returned on energy invested. This combination of uncertainty added into the simulation for the VRP required large computational effort. To make the problem more computationally manageable, the 25 sawmills were spatially clustered based on the travel distance matrix using the K-medoids technique into 4 distinct groups and the itinerary was then rerouted to the centre of each cluster.

Based on the above, clustering algorithms play a vital role in addressing challenges associated with upstream, midstream and downstream supply chain stages in biomass supply chain optimisation, where typically a large number of distinct locations are involved in real-life systems. They can identify potential processing hubs and deal with large datasets in MILP models with high complexity. Clustering algorithms, such as GIS proximity-based clustering, K-means, K-medoids, C-means prove to be essential in creating manageable data subsets, enabling the reduction of computational burden and ensuring reasonable solution times, allowing the construction and use of decision support systems suitable for business applications.

Very limited applications of the Weighted-MSSC problem can be found; one of them presenting the weighted version of the classical K-means algorithm for a use-case regarding locating distribution centres for COVID-19 vaccines throughout the US (https://towardsdatascience.com/determining-optimal-distribution-centers-locations-using-weighted-k-means-1dd099726307 [Online] (accessed April 2, 2024)). We use the related dataset in our computational results section. Similarly, the query ‘weighted MSSC clustering’ returns a number of papers on clustering ensembles where different clusters have different weights, but individual data points have no weights assigned to them as in our case. Overall, even though there have been a number of papers advancing the state-of-the-art in MSSC, with the exception of the article [https://towardsdatascience.com/determining-optimal-distribution-centers-locations-using-weighted-k-means-1dd099726307 [Online] (accessed April 2, 2024)] none have studied the weighted version of the problem. In the next three paragraphs, therefore, we revisit a number of papers on the MSSC problem instead.

In Arthur and Vassilvitskii (2007), the authors present the well-known ‘K-Means++’ heuristic for initialising the centres of the K-means algorithm, by modifying another popular heuristic, the so-called ‘farthest centers’ heuristic, which works by selecting at random one point from the dataset as first centre, and then repeatedly chooses the point that is ‘the furthest distance’ away from the current set of centres as the next centre. Because this heuristic was prone to choosing outliers in the dataset as the initial centres, the K-Means++ algorithm works by selecting at random but with probability proportional to the shortest distance of any point from the currently selected centres the next point to be selected as the next centre. This randomisation procedure guarantees theoretical good properties of the initial clustering for the K-means algorithm. Then, in Bahmani et al. (2012), the authors remove the serial bottleneck of the above-described procedure by fully parallelising it to allow utilising as many CPU cores as there are available in a cluster, and therefore allow much larger datasets than was possible before to be clustered efficiently; this K-Means variation is denoted as the ‘K-Means ||’.

Aloise et al. (2012) presented an improved column generation method for the MSSC that was able to solve to optimality problems with up to 2400 points in 2-D Euclidean space and K taking values up to 400. Very recently, the authors of (Piccialli et al., 2022) presented a Semi-Definite Programming approach (SDP) to the MSSC, that allows them to solve to provable optimality a number of problems that were not solvable before. In (Burgard et al., 2023), the authors use advanced mixed integer programming techniques such as cutting planes (including cardinality cuts, outer approximation cuts and bary-centre propagation) to enhance the performance of the Open-Source solver SCIP.

Peng and Xia (2005) present a ‘convexity cut method’ to solve to global optimality the MSSC problem. They use CPLEX for solving their corresponding cutting plane problem; they report results for problems up to 4600 data points, but they only work with very small values of K. Finally, in (Usman, 2023), the author presents a number of clustering algorithms for forest harvesting, mentioning the classical K-means algorithm as one possible method for solving their (unweighted) problem.

According to the literature review, clustering is not fully exploited under the concept of comprehensive biomass supply chain management. In the upstream level, where typically biomass fields of varying size exist, the existing clustering approaches are employed only to treat the location dimensionality while the hectarage (area) factor is ignored. The absence of incorporating hectarage in the cluster calculation may result in forming groups that don't adequately represent the initial set of fields, diminishing accuracy and affecting many aspects of the design and optimisation of the whole biomass supply chain such as the transportation costs, the location and capacity of midstream processing facilities, etc. The present study aims to fill this methodological gap and develop an effective clustering algorithm that allows a more precise representation of clustered sites considering both their size and location. Enhancing and improving the clustering process facilitates the management of more complex problems without increasing solution times, while preserving accuracy and generating robust results which support a more informed decision-making.

3. Model formulation

The clustering problem dealt in this study is the Weighted MSSC problem, which generalises the well-known MSSC problem (Hansen & Mladenović, 2001; Theodoridis & Koutroumbas, 2006) as follows: consider a finite set of points $S=\left\{s_1,\dots ,s_n\right\}$ in $d$-dimensional plane, with associated positive weights $w_1,\dots w_n$. The problem is to find $K>1$ ‘centers’ in the plane that will minimise the total weighted sum of square distances of each of the $n$ points to their closest centre multiplied by the point's corresponding weight. It is worth noting that when all the points have equal weights, the problem reduces to the well-known minimisation of sum of square errors problem that appears in many problems in signal processing (code-book design), unsupervised learning (K-Means clustering), supply chain design, and so on.

The problem can be stated then as follows:

We seek to find a partition $C$ of the set $S$ into $K$ disjoint sets $C_i,i=1,\dots K:C_i\cap C_j=\mathrm{\varnothing \forall i}\ne j,i,j=1\dots K,\underset{i=1}{\overset{K}{\cup }}{C}_i=S$ that minimises the function $f\left(C\right)={\sum }_{i=1}^K{\sum }_{s_m\in C_i}{w}_m\parallel s_m-\overline{s_i}{\parallel }^2$ where the point $\overline{s_i}$ is such that it minimises the sum of weighted distances of each point in the cluster $C_i$ to it.

Lemma 1:

Given the points $s_m$ in a partition cluster $C_i$, with accompanying weights $w_m$, the optimal point in the plane that minimises the sum of weighted distances of the points in the partition cluster to it, is the (weighted) barycenter $\overline{s_i}={\sum }_{s_m\in C_i}{w}_m{s}_m/{\sum }_{m:s_m\in C_i}{w}_m$.

Proof:

It is enough to show that for the function (1) $f\left(C_i,x\right)=\sum _{s_m\in C_i}{w}_m\parallel s_m-x{\parallel }^2$(1) the gradient ${\mathrm{\nabla }}_{x}f\left(C_i,{\overline{s}}_i\right)=0$. The gradient of the function at any point $x$ is (2) ${\mathrm{\nabla }}_{x}f\left(C_i,x\right)=-2\sum _{s_m\in C_i}{w}_m\left(s_m-x\right)$(2) Setting the above expression to zero, we obtain as the only solution the vector (3) $x^{\ast }=\overline{s_i}=\frac{{\sum }_{s_m\in C_i}{w}_m{s}_m}{{\sum }_{m:s_m\in C_i}{w}_m}$(3)

Therefore, mathematically, we seek to optimise the following problem (P): (P) $\underset{C}{min}\sum _{i=1}^K\sum _{s_m\in C_i}{w}_m\parallel s_m-\overline{s_i}{\parallel }^2$(P) subject to: $\{\begin{array}{ll}\overline{s_i}=\frac{{\sum }_{s_m\in C_i}{w}_m{s}_m}{{\sum }_{m:s_m\in C_i}{w}_m},& i=1\dots K\\ C_i\cap C_j=\mathrm{\varnothing },& \mathrm{\forall i}\ne j,i,j=1\dots K\\ \underset{i=1}{\overset{K}{\cup }}{C}_i=S\end{array}$

4. A column generation approach to solving the weighted MSSC problem

The minimum sum of square errors clustering problem (MSSC) has been the subject of intense study for many decades (Theodoridis & Koutroumbas, 2006). It is a combinatorial problem by nature: given a finite set $S=\left\{s_1,s_2,\dots s_n\right\}\in {\mathbb{R}}^d$of $n$ data points in some vector space, the problem is to find a partition of the data points in exactly $K>1$ disjoint partition clusters so that the sum of the square distances of each point to its closest centre (cluster mean point) is minimised. The problem is known to be NP-Hard (Christou, 2012). The MSSC problem is clearly a special case of the weighted MSSC problem where the weight of each data point is set to 1.

The most known and frequently used algorithm for solving the standard MSSC problem is the K-Means algorithm (Aloise et al., 2012; Burgard et al., 2023; Peng & Xia, 2005; Usman, 2023), an expectation-maximisation (EM) type algorithm that works by iterating two major phases: an expectation step that assigns points to their closest current centroid (mean), and a maximisation step that re-computes the centroids as the mean vector of every data point in a cluster, optimising the objective function under the current assignments. The algorithm iterates until some convergence criteria is met, or -for very large problems- for a fixed number of iterations. For most problems, the initial choice of cluster centres is a major factor that can determine the quality of the final solution as well as the speed of convergence (Arthur & Vassilvitskii, 2007; Bahmani et al., 2012).

K-means is a heuristic algorithm that does not guarantee the optimal solution of an MSSC problem, even though it often finds high quality solutions, at least for modest values of $K$. As described above, the K-means algorithm does not solve the weighted MSSC problem, however, it is not difficult to state a variant that attempts to. Essentially, the only modification to the K-means algorithm that is needed, is in the M-step (step 2) where we compute the new cluster centres. Given Lemma 1 in section 2, we see that the modification needed to be done to the original K-means algorithm to result in a weighted K-means algorithm is to compute the centre of each cluster in step 2 according to equation (3).

However, whereas the K-means algorithm, and the modified weighted K-means algorithm usually provide acceptable results for problems with a small number of partitions required (say, up to a few tens), the same is not usually the case for a larger number of partitions (Laszlo & Mukherjee, 2006). For such problem instances, we propose a column-generation type approach to a set covering formulation of the problem, in a spirit similar to the approach described in (Christou, 2011) for the standard MSSC problem.

In particular, notice first that the weighted MSSC problem, being also a combinatorial optimisation problem, can be written as a set partitioning problem. Consider all possible subsets of the set $S$ and denote by $A$ the $n\times 2^n$ matrix containing ones and zeros, whose columns represent all those possible subsets: every column of the matrix has exactly $n$ components, and the $i$-th component indicates whether the data point $s_i$ is contained or not in the subset. Given this (rather big) matrix, the weighted MSSC problem can be written down as the set-partitioning problem. (SPP) $\begin{array}{rl}& \underset{x}{min}{c}^{T}x\end{array}$(SPP) $\begin{array}{rl}& \text{subject to}:\{\begin{array}{l}\mathrm{Ax}=e\\ e^{T}x=K\\ x\in {\left\{0,1\right\}}^n\end{array}\end{array}$where the column vector $e\in {\mathbb{R}}^n\equiv [1,1,\dots 1{]}^T$ and the cost-vector $c=[c_1,c_2,\dots c_{2^n}{]}^T\in {\mathbb{R}}^{2^n}$ contains the (weighted) costs of all subsets of the set $S$. The set-partition problem formulation (SPP) asks to pick exactly $K$ subsets of the set $S$ (constraint $e^{T}x=K$ and $x_i\in \left\{0,1\right\}\text{for all}i=1,\dots 2^n$) such that they completely cover the set $S$ without any overlap (constraints $\mathrm{Ax}=e$) and minimise the total weighted cost of each selected subset.

Solving the problem (SPP) is not any easier than solving the problem (P), and both problems are NP-Hard. Even further, just writing down the initial model formulation of the problem (SPP) requires explicitly enumerating and storing all $2^n$ possible subsets of $S$ which is clearly infeasible for even small values of $n$. We, therefore, resort to a column-generation approach (originally described in (Christou, 2011) for the standard MSSC problem) whereby we begin with a base matrix $B$ that contains only a few of the columns of the full matrix $A$, and afterwards we solve the following variant of the initial problem (SPP), called (SCP_R): $\begin{array}{rl}& \underset{x}{min}{c}^T{x}_R\left({\mathrm{SCP}}_R\right)\end{array}$ $\begin{array}{rl}& \text{subject to}:\{\begin{array}{l}\mathrm{Bx}\ge e\\ e^{T}x=K\\ x\in {\left\{0,1\right\}}^n\end{array}\end{array}$

The components of the column-vector $c$ again correspond to the weighted clustering cost of the corresponding subsets described in the columns of the reduced matrix $B$, only now we seek to solve the set covering version of the problem (constraints $\mathrm{Bx}\ge e$) which means we don't care if some data points appear in more than one subsets selected in the final solution, as long as each point appears in at least one selected subset.

In case there exists at least one data point $s_{\mathrm{dup}}$ that appears in more than one selected subset in the solution of (SCP_R), we remove it from every selected subset in which it appears except the subset from which if we removed $s_{\mathrm{dup}}$ the weighted cluster cost increase would be the smallest among all other selected subsets that contain it; ties are arbitrarily broken. It is easy to see that such a procedure results in strictly lowering the objective function value of the weighted MSSC problem.

Given the new clusters that represent a valid clustering solution to the weighted MSSC, though it may or may not be optimal, we generate more columns to add to our base matrix $B$ by adding the new subsets that we created in the conversion process described above (if any), and then by running the weighted K-means algorithm on the solution found by the previous step, and adding all new clusters created during the run of the weighted K-means algorithm to the new base matrix $B$. We then solve the new problem (SCP_R) with the new base matrix $B$ again, and repeat this process until in one full iteration we no longer improve on the final objective value of the problem.

The entire algorithm, called W-EXAMCE, is a variant of the EXAMCE algorithm specified in (Christou, 2011) adapted to work with weighted clustering, is specified in Figure 1 and is available from https://www.github.com/ioannischristou/weighted_clustering.

Figure 1. W-EXAMCE algorithm.

The procedure Expand(C,τ) expands the set of available solutions to consider in step 5 by producing a (possibly empty) set of clusters that are perturbations of the input cluster C. We use a simple heuristic whereby we choose the τ closest neighbours to the centre of C that are not in C, as well as the τ members of C that are farthest from the centre of C, with τ being a small user-defined parameter (set after experimentation equal to 5). Distances here are weighted, meaning that the distance between a cluster centre $c$ and a data point $s_i$ with associated weight $w_i$ is the quantity $\surd w_{i}c-s_i$. We create and return a sequence of 2τ new clusters $N_1^{+}\subset N_2^{+}\dots \subset N_{\tau }^{+},N_1^{-}\supset N_2^{-}\dots \supset N_{\tau }^{-}$ where each of the $N_i^{+}$ clusters contains all points in C plus up to the $i$-th closest non-member neighbour of C, and each of the $N_i^{-}$ clusters contains all points in C except up to the $i$-th farthest member of C.

It is easy to see that in every iteration in steps 3–10 of the W-EXAMCE algorithm, the solution C’’ that is found is at least as good as the corresponding solution from the previous iteration:

Lemma 2.

The solutions C’’ found in step 7 of the W-EXAMCE algorithm are non-increasing.

Proof

Notice that the columns of the matrix B in problem (SCP_R) to be solved at each iteration always form a super-set of the columns of the same matrix in the previous iteration, even after that previous matrix is augmented by the solution found in the previous execution of step 7. Therefore, the solution of the problem in step 5 is always at least as good as the solution of step 7 found in the previous iteration. And since both the Rm_Dup(.) procedure in step 6 by design can only further reduce the objective value of the problem, and the same is true for the application of the W-K-Means algorithm in step 7, the solution C’’ found in step 7 at each iteration is always monotonically decreasing.

Because of the above monotone convergence of the algorithm (convergence is guaranteed since the solution sequence is monotonically decreasing, and is bounded from below by zero), we have the following:

Corollary 3.

The W-EXAMCE algorithm converges in a finite number of iterations.

Proof

The convergence of the algorithm is already established. The finiteness of the number of steps follows from the fact that there are only finitely many different partitions to partition the data, and the algorithm stops as soon as no progress is made in two consecutive iterations.

A small worked example of the W-EXAMCE algorithm is presented in Appendix A.

5. Computational results

In all experiments below, we use the GUROBI optimiser (v. 11) to solve the SCP problem in step 5 of the W-EXAMCE algorithm in Figure 1. We use a PC with an Intel Core i9-10920X CPU running at 4.5 GHz, equipped with 64GB RAM, running Windows 11. Both the standard K-means and the scalable K-Means|| algorithm implementations are fully parallelised using Java threads and utilise all 24 virtual CPU cores that the CPU offers. All algorithms are written in Java.

In Table 1, the meta-data for the major datasets used for our experiments are presented. All datasets represent points in the Euclidean 2-dimensional plane. Due to the lack of publicly available biomass-related land plot datasets, several more datasets from well-known libraries and other sector applications were used to demonstrate the performance of the proposed method. In Table 2 we present the results obtained from 4 datasets in the TSPLIB that are commonly used to benchmark algorithms for the standard MSSC problem (see Christou, 2011; Hansen & Mladenović, 2001; Laszlo & Mukherjee, 2006), with attached random positive weights for each data point. We compare results from applying standard K-means with 100 restarts, with the scalable K-means|| algorithm with 100 restarts, and the W-EXAMCE algorithm that uses the clusters obtained from 10 restarts of the standard K-means algorithm as base clusters. The column labelled ‘Soln Improvement%’ under the W-EXAMCE results columns lists the percentage improvement of W-EXAMCE over the baseline K-means algorithm. As can be seen from the results, W-EXAMCE is always the clear winner of the three algorithms even though it takes a little more time to complete.

Table 1. Used dataset characteristics.

Dataset	Library	Studies using the dataset	Number of points	Number of dimensions	Weights
u1060	TSPLIB	(Christou, 2011; Laszlo & Mukherjee, 2006)	1,060	2	Random
pcb3038	TSPLIB	(Christou, 2011; Laszlo & Mukherjee, 2006)	3,038	2	Random
pr2392	TSPLIB	(Aloise et al., 2012)	2,392	2	Random
rl11849	TSPLIB	(Aloise et al., 2012; Christou, 2011)	11,849	2	Random
Scotland Vacant and Derelict Land Register	Land Register	Psathas et al. (2022)	3,398	2	Actual Plot Areas
US COVID-19 Cases	N/A	https://towardsdatascience.com/	4,478	2	Actual number of cases per district
		determining-optimal-
		distribution-centers-locations-
		using-weighted-k-
		means-1dd099726307
		[Online] (accessed April 2, 2024)

Table 2. Comparison of results of W-EXAMCE, K-Means and scalable K-Means|| w/ 100 restarts on TSPLIB datasets.

		Best K-means / 100 restarts		Best K-means|| / 100 restarts		W-EXAMCE
Dataset	k	Solution	Time (secs)	Solution	Time (secs)	Solution	Solution Improvement (%)	Time (secs)
rl11849	100	1.74E+09	7.78	1.73E+09	10.2	1.72E+09	1.5	155.0
	200	8.36E+08	8.9	8.24E+08	13.6	8.13E+08	2.8	242.0
	300	5.40E+08	9.9	5.30E+08	16.7	5.14E+08	5.2	133.0
	400	3.95E+08	10.24	3.82E+08	19.5	3.68E+08	7.4	177.0
	500	3.11E+08	10.3	2.94E+08	24	2.81E+08	10.4	707.0
pr2392	50	4.74E+08	3.8	4.81E+08	3.98	4.70E+08	0.8	21.0
	100	2.23E+08	3.28	2.18E+08	4.06	2.02E+08	10.2	26.5
	200	9.67E+07	1.77	9.57E+07	5.45	8.43E+07	14.7	36.6
	300	6.02E+07	1.94	5.64E+07	5.41	4.73E+07	27.4	73.1
	400	4.09E+07	1.98	3.75E+07	6.19	3.01E+07	36.2	40.7
	500	2.99E+07	1.99	2.64E+07	7.94	2.08E+07	44.1	47.9
u1060	50	2.97E+08	2.1	2.90E+08	2.62	2.72E+08	9.1	30
	100	1.32E+08	3	1.31E+08	3.5	1.05E+08	25.6	11.5
	150	8.48E+07	3	7.96E+07	4.5	6.45E+07	31.4	20.4
	200	6.14E+07	4.5	6.15E+07	4.6	4.21E+07	46.0	16.1
	250	4.42E+07	4.1	4.13E+07	4.8	3.02E+07	46.3	20.4
	300	3.41E+07	4.9	3.26E+07	5.5	2.30E+07	48.2	21.4
pcb3038	50	1.05E+08	5.1	1.04E+08	5.5	4.94E+07	113.4	15.3
	100	5.29E+07	6	5.09E+07	6.1	2.39E+07	121.4	36
	200	2.59E+07	12	2.08E+07	12	1.06E+07	144.9	16
	300	1.65E+07	22.1	1.02E+07	23.1	6.21E+06	166.3	25.4
	400	1.25E+07	20.1	1.00E+07	22.9	4.12E+06	203.6	30.8
	500	9.84E+06	25	7.53E+06	26.2	2.87E+06	242.9	34

In Table 3, we present results from one large-scale use-case for phyto-remediation-based supply chain management, with real land plot data from the Scottish Vacant and Derelict Land Register (SVDLS, 2022): the dataset contains the coordinates of 3,398 previously-developed brownfield land plots, in the British National Grid Coordinate system, a plane coordinate system that is based upon a Transverse Mercator projection. The dataset also contains the area that each plot covers that serves as the weight of each data point for our clustering purposes. This dataset refers to one of the phytoremediation use cases of the project CERESiS (CERESiS, 2023), where the objective was to optimise the respective supply chain of energy crops grown in potentially contaminated land to produce biofuels, while at the same time decontaminating the soil. The dataset is visualised in Figure 2, however the size of each land plot is not drawn to scale.

Figure 2. Locations of Vacant and Derelict Land Register Sites in Scotland in the British Coordinate Grid. The size of each plot is not drawn to scale.

Table 3. Comparison of results of W-EXAMCE, K-Means w/ 100 Restarts and Scalable K-Means|| w/ 100 Restarts on Scotland Vacant and Derelict Land Register dataset.

		Best K-Means / 100 Restarts		Best K-Means|| / 100 Restarts		W-EXAMCE
Dataset	k	Solution	Time (secs)	Solution	Time (secs)	Solution	Solution Improvement (%)	Time (secs)
Scotland Vacant and Derelict Land Register	100	2.50E+11	4.5	2.11E+11	7.8	1.68E+11	48.8	37.8
	200	1.49E+11	19.8	1.01E+11	25.3	6.05E+10	146.3	80.9
	300	8.89E+10	16.9	9.91E+10	30.1	2.37E+10	275.1	77.9
	400	6.95E+10	21.8	4.99E+10	25.9	1.44E+10	382.6	91.8
	500	4.94E+10	24	2.51E+10	28.5	8.64E+09	471.8	126.2

Finally, to assess the algorithm on another real-life large dataset, in Table 4 we present results using the public (weighted) dataset created for determining the best locations for placing distribution and vaccination centres for the COVID-19 pandemic in the U.S.A. This dataset contains 4,478 data points in 2D. We have converted the original standard lat-lon coordinates of the dataset into a flat grid using Mercator Projection (UTM), and normalised the weights of the original dataset to sum up to 1 (otherwise, the sum of the weighted square distances would lead to arithmetic overflows due to resulting large numbers.)

Table 4. Comparison of results of W-EXAMCE, K-Means w/ 100 Restarts and Scalable K-Means|| w/ 100 Restarts on the US COVID-19 Cases dataset.

		Best K-Means / 100 Restarts		Best K-Means|| / 100 Restarts		W-EXAMCE
Dataset	k	Solution	Time (secs)	Solution	Time (secs)	Solution	Solution Improvement (%)	Time (secs)
US COVID-19 Cases	50	1.16E+10	8	1.06E+10	37	9.63E+09	20.1	30.6
	100	4.92E+09	7.5	4.35E+09	34.9	3.70E+09	33.0	35.7
	200	1.94E+09	9.8	2.16E+09	30.9	1.31E+09	47.7	51.3
	300	1.16E+09	9.1	1.31E+09	37	6.66E+08	74.5	71.5
	400	7.86E+08	12.7	9.00E+08	34.7	4.03E+08	94.8	99.7
	500	5.55E+08	11.6	6.31E+08	37.5	2.73E+08	103.1	103

In Figure 3 we show the percentage gap in solution quality between the proposed W-EXAMCE and the weighted K-means with 100 restarts on the Scotland Vacant and Derelict Land Register data; as can be seen clearly both from Figure 3 and from Tables 2–4, the solution found by W-EXAMCE is always superior to that found by 100 restarts of the K-means or the scalable K-Means|| algorithms. It can also be observed that the superiority of the W-EXAMCE algorithm is greatly enhanced when the number of partitions increases, while it requires more time to be solved; however, the time needed is reasonable for practical applications. By employing the proposed algorithm, a more efficient clustering is achieved, in order to support the supply chain optimisation process that uses the clusters as input.

Figure 3. Comparison chart of the solution gap between W-EXAMCE and weighted K-Means with 100 Restarts on the Scotland Land Register dataset.

6. Conclusions and future directions

This study proposed a novel, more effective clustering method for use within biomass supply chain optimisation modelling, as well as for any other case where clustering is required within supply chain analysis and optimisation context. In the context of biomass supply chain analysis, clustering is not often utilised; instead, clustering primarily finds applications in VRP problems (such as the ‘cluster-first-route-second’ heuristic by Fisher and Jaikumar (1981)), where weights are not usually taken into account.

The proposed algorithm (W-EXAMCE) for the weighted MSSC problem proved to be very effective in identifying more representative clusters of biomass fields than the commonly used technique of the weighted K-means and the scalable K-Means II. The solution quality gap between the proposed algorithm against the weighted K-means increases significantly as the number of partitions increases. Therefore, when more refined supply chain analysis is required, the proposed algorithm can offer more accurate clusters, leading to more accurate transportation cost estimates, hence more accurate supply chain optimisation outcomes.

The applicability of the proposed algorithm for the weighted MSSC problem (W-EXAMCE) in biomass supply chain modelling, where exact locations of fields or biomass hubs are known, has been demonstrated for the Scottish Vacant and Derelict Land Register dataset (used as a use case in the CERESiS project (CERESiS, 2023), the US COVID-19 cases dataset [https://towardsdatascience.com/determining-optimal-distribution-centers-locations-using-weighted-k-means-1dd099726307 [Online] (accessed April 2, 2024)], and several TSPLIB datasets (with random weights). The proposed algorithm offers a viable option for increasing the optimisation results’ accuracy. Given the spatial dispersion, typically low availability of biomass feedstock, and diverse land hectarage of the fields under analysis, a refined approach is necessary, since simple aggregation and clustering techniques may introduce inaccuracies in the results. Here, the weighted K-Means algorithm provides robust solutions for typical biomass to biofuels or bioenergy supply chains, offering a detailed grouped view of the fields based on hectarage that is essential for effective decision-making.

The proposed W-EXAMCE algorithm could have wider applicability in all supply chain modelling contexts where many upstream, midstream or downstream nodes are involved, provided that each node of the physical system carries a different weight. In this case, using the proposed algorithm would lead to more accurate clustering compared to the typical K-means algorithm, ultimately leading to supply chain optimisation IEA, 2017 outputs that are more relevant and based on more accurate inputs.

In the future, the W-EXAMCE algorithm could be applied in more diverse datasets of biomass supply chains, to understand its performance gap and expected benefits in a wider set of conditions. Its applicability could also be tested in other supply chain contexts to demonstrate its potential for generalisation.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data in this study can be provided upon request from the authors by emailing them. The entire W-EXAMCE algorithm is available from https://www.github.com/ioannischristou/weighted_clustering.

Additional information

Funding

This study is part of the CERESiS project which has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement No 101006717.

Appendix A

For a better understanding of the workings of the W-EXAMCE Algorithm described in Figure 1, we demonstrate its application with a small example.

In Figure A1 below, we show a dataset of 10 points in 2 dimensions. Each data-point has an associated weight according to Table A5 below, and the ‘X’ that is used to draw each point in Figure A1 is drawn approximately according to the scale shown in Table A5.

Figure A1. Toy dataset with 10 Points in 2D to trace the steps of algorithm W-EXAMCE.

Table A5. Toy dataset weights.

Point index	Weight
1	0.01
2	0.01
3	0.01
4	0.45
5	0.45
6	0.02
7	0.01
8	0.01
9	0.02
10	0.01

We are going to trace the execution of the algorithm W-EXAMCE with k = 4. For simplicity, we also set τ = 0 for our procedure $\mathrm{Expand}\left(C,\tau \right)$ so that in essence the procedure does not do anything. We produce the initial matrix $A_B$ by executing the clusterers in the set Ba = (W-KMeans(4), W-KMeans(4), W-KMeans(4), W-KMeans(3), W-KMeans(3), W-KMeans(5), W-KMeans(5)); this means that we apply 3 times the W-KMeans algorithm with k = 4 (starting from a different random initial solution each time), and then we apply 2 times the W-KMeans algorithm with k = 3, and with k = 5 (step 1). This results in a total of 28 clusters whose members and total weighted costs are shown in Table A6 below. The clustering solutions are colour-coded in Table A6 below, so that the first clustering solution comprises of the clusters 1–4 (inclusive), the second solution comprises of clusters 5–8 and so on.

Table A6. Results of applying clusterers in the set Ba to the dataset.

Cluster index	Cluster members	Weighted cluster cost
1	2, 4, 6, 10	5.73E-04
2	1, 5, 7	9.79E-05
3	9	0
4	3, 8	5.00E-05
5	3, 5, 7	1.00E-04
6	8, 9	1.33E-04
7	2, 4, 6, 10	5.73E-04
8	1	0
9	1, 2, 4	9.79E-05
10	6, 10	1.33E-04
11	3, 5, 7	1.00E-04
12	8, 9	1.33E-04
13	1, 3, 5, 7	1.49E-04
14	8, 9	1.33E-04
15	2, 4, 6, 10	5.73E-04
16	1, 2, 3, 4, 5, 6, 7, 8	5.30E-03
17	10	0
18	9	0
19	2, 4, 6	1.45E-04
20	10	0
21	5, 7	4.89E-05
22	8, 9	1.33E-04
23	1, 3	5.00E-05
24	10	0
25	3	0
26	1, 5, 7	9.79E-05
27	2, 4, 6	1.45E-04
28	8, 9	1.33E-04

The best W-KMeans(4) solution found in this step has value 4.65E-04, calculated as the sum of the values in the column ‘Weighted Cluster Cost’ for rows 9–12.

Given the data in the 2nd column of Table A6, the constraints of the problem in step 5 are then as follows:

Figure A2. Matrix representation of the set covering problem formulation constraints.

Where we have only written the matrix $A_B$ as the LHS and the values of the RHS vector of the constraints.

When the MIP solver solves the SCP problem with the above data, the best solution it finds picks the columns 6, 13, 17 and 19 which represent the following clusters:

{8, 9}

{1, 3, 5, 7}

{10}

{2, 4, 6}

This solution has total weighted cost equal to 4.27E-04 ( = 1.33E-04 + 1.49E-04 + 0 + 1.45E-04) which happens to be less than any cost found by the application of the base clusterers for $k=4$. Because the solution found is also a solution to the set partitioning problem (SPP), the step 6 is a no-op (does not do anything). Similarly, step 7 (applying W-KMeans(4) from the solution just found) does not create any new clusters because the new solution found is already a local optimum for the W-MSSC problem. Therefore (given that τ = 0) steps 8–10 are also no-ops, and no new columns are added to the matrix $A_B$. Therefore, in the next iteration of steps 3–10, no improvement will occur, and the algorithm stops after two iterations of its main loop. The best solution found by applying the W-EXAMCE algorithm is shown in Figure A3 below.

Figure A3. Best clustering solution found by W-EXAMCE on toy dataset.