A Petri Net-based Algorithm for Solving the One-Dimensional Cutting Stock Problem

1. Introduction

Cutting stock problems are combinatorial problems involving the creation of small objects from large objects according to a cutting plan made of cutting patterns. The larger the number of small objects, the greater the number of possible cutting patterns, and hence these problems are classified as NP-hard optimization problems [1]. In this study, we focus on the single stock-size cutting problem according to Wascher’s classification [2] for the one-dimensional case, or one-dimensional cutting stock problem (1D-CSP), as it is best known in most of the literature. The 1D-CSP arises in industries manufacturing goods from paper, metal, glass, and wood, among others. The goal is the reduction of material usage and, therefore, the minimization of the related costs. Apart from cost reduction, industries such as the building industry are concerned with reducing gas emissions due to metal and concrete usage, as pointed out in [3]. Some other studies also address the problem by considering mathematical models that include the reuse of leftovers, as in [4].

Most approaches to solving the 1D-CSP are based on linear programming (LP) methods, with the column generation technique being one of the most outstanding ones. Reference [5] presents some LP formulations proposed in the early stages of addressing the 1D-CSP. Additionally, a comprehensive survey of LP methods used to solve the 1D-CSP up to 2015 can be found in [6]. Column generation and other techniques focus on generating and reducing the number of cutting patterns. In the 1D-CSP, a cutting pattern is an ordered sequence of items cut from a large object or stock. One-cut models and arc-flow models are other paradigms used to address the 1D-CSP, as solutions to the 1D-CSP can be constructed as sequences of items. It is well-known that LP methods can yield optimal solutions to the 1D-CSP, provided that the complexity of the instances is small, that is, with a small number of small objects or items. This limitation has led to the use of heuristic and metaheuristic approaches. Although these methods usually result in near-optimal solutions, the computational time is typically higher compared to LP methods. Additionally, a drawback of using metaheuristics is the complexity of solution representation and the numerous parameters that must be set and tuned. Since a solution for the 1D-CSP can be sequentially constructed, we propose using ordinary Petri nets (PN) to model the solution construction. We take advantage of the reachability tree to find the best solution through the implementation of a filtered beam search algorithm (FBS). Beam search algorithms are a type of any-time search algorithm, suitable for optimization problems where solutions can be constructed sequentially and represented as a tree structure. Unlike other any-time search methods, beam search algorithms do not perform backtracking. Only the most promising nodes at each level are retained to continue searching for the best path from the root node to a desired final node, while non-promising nodes are pruned permanently. The FBS is an improvement to the classical beam search algorithm. At each level, it first filters the nodes through a local evaluation and then performs a global evaluation among the filtered nodes to retain the best ones for further investigation.

In the literature, most mathematical models for the 1D-CSP consider single-objective functions. However, [7] demonstrated the effectiveness of using a bi-objective function. In this study, we consider the bi-objective function of the mathematical model (1)-(3) based on the one used in [8] for the 1D-CSP without contiguity. The objective function to be minimized accounts for both the minimization of waste and the minimization of the number of stocks with waste.

$$\begin{array}{cc}\hfill \min & W={\displaystyle \frac{1}{k+1}}\left(\sum _{j=1}^k\sqrt{{\displaystyle \frac{w_j}{L_j}}}+\sum _{j=1}^k{\displaystyle \frac{V_j}{k}}\right)\hfill \end{array}$$

(1)

$$\begin{array}{cc}\mathrm{subject}\mathrm{to}\hfill & \\ \hfill & \hfill \sum _{j=1}^k{x}_{ij}=d_i,i=1,2,\dots ,m\end{array}$$

(2)

$$\begin{array}{ccc}\hfill & \hfill & \hfill x_{ij}\ge 0\mathrm{and}\mathrm{integer}.\end{array}$$

(3)

In model (1)-(3), $L_j$ is the length of the stock, m is the number of different items, k is the total number of stocks used. $d_i$ represents the demand or total number of orders of item i, while $l_i$ is the length of item i, $x_{ij}$ is the number of orders of item i in stock j. $w_j$ and $V_j$ are defined in Equations (4) and (5), respectively, where $w_j$ represents the wastage of stock j.

$$\begin{array}{cc}\hfill w_j& =L_j-\sum _{i=1}^m{x}_{ij}{l}_i,j=1,2,\dots ,k\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill V_j& =\left\{\begin{array}{cccc}1\hfill & \mathrm{if}& w_j>0,\hfill & j=1,2,\dots ,k\\ 0\hfill & & \mathrm{otherwise}.\hfill \end{array}\right.\hfill \end{array}$$

(5)

Even though the main drawback of PN is state explosion, we do not generate the complete reachability tree. Instead, we generate promising states from the initial state that lead to the desired final state with minimum waste, according to the FBS, using the bi-objective function (1) to evaluate the nodes. The initial state or marking of the PN is where no items have been taken to form a solution, and the final state or marking is where all items have been taken to form a solution. At each level of the tree, partial solutions are evaluated locally, and the best ones are then evaluated globally. Nodes with the best partial solutions are kept for further investigation. The advantages of using PN are twofold: they are useful for validating solutions when PN is deadlock-free, and they can be combined with graph search methods to solve optimization problems. Some of the problems addressed under the PN paradigm include job and flow shop scheduling problems. In these cases, PN and graph search methods like ${\mathrm{A}}^{*}$ algorithms [9,10,11,12] and beam search algorithms [13,14,15,16] have been used to obtain the best scheduling solutions. However, to the best of our knowledge, PN has not been used to solve cutting problems. Our experimental results show that PN is suitable for modeling pattern generation and solving the 1D-CSP, with acceptable performance in both solution quality and computational time.

2. Literature Review

It is widely considered that Kantorovich first addressed the 1D-CSP in [17] providing an LP formulation of the problem and using the multipliers method to solve it. Later, [18,19] used the column generation technique to solve the integer linear programming formulation by addressing an associated knapsack problem. Following the introduction of the column generation technique, several procedures were developed to improve its computational performance using branch-and-price algorithms [20,21,22,23]. Other algorithms have focused on constructing and reducing the number of patterns [24,25,26,27,28]. Alternative approaches to modeling and solving the 1D-CSP include the arc-flow and one-cut formulations. The arc-flow formulation, a graph-based approach, was introduced in [29]. Recently [30] proposed and compared arc-flow formulations based on other algorithms. One-cut formulations have received little attention since their use in [31]. Moreover, [32] proved the equivalence among pattern-based formulations, such as that of Gilmore & Gomory [18,19], and the arc-flow and one-cut formulations. [33] used dynamic programming to solve the 1D-CSP, and [34] showed the relationship between arc-flow formulations and dynamic programming in solving optimization problems.

LP methods are restricted to instances of a small number of items, and to overcome this limitation, many researchers have proposed metaheuristic methods like genetic algorithms [35,36,37]. However, the use of genetic algorithms faces the difficulty of finding the best way to encode the solutions or cutting patterns, as pointed out in [38,39,40,41]. Evolutionary programming was used in [8,42] to address the 1D-CSP with and without contiguity. Additionally, [43] concluded that evolutionary programming performs better than genetic algorithms for solving the 1D-CSP. Swarm intelligence algorithms, such as particle swarm optimization algorithms [44,45,46,47] and ant colony optimization algorithms [48,49,50], have also been used to solve the 1D-CSP. The main issues with these two algorithms are parameter tuning and the need for hybridizations to avoid premature convergence and control the search space. Additionally, discretization procedures are required since these algorithms are better suited for continuous problems. Other metaheuristic algorithms for solving the 1D-CSP include tabu search, simulated annealing [51], and iterative local search [52,53]. More recently, [54] introduced the least-lost algorithm to simultaneously minimize waste and stocks with waste in the 1D-CSP, showing better results than the evolutionary programming algorithm of [8]. [55] compares a Generate & Solve algorithm to the column generation technique. The Generate & Solve algorithm uses a reduction technique based on the one presented in [56] to deal with large instances of the 1D-CSP. Although the algorithm is not better than the column generation technique, it can solve larger instances in acceptable computational time with quasi-optimal solutions. [57] introduces a deep reinforcement learning algorithm to minimize the number of stocks and maximize the length of leftovers. The algorithm performs well in terms of computational time and handling large instances.

3. Basic Concepts

In this section, we present the basic concepts of ordinary PN and the beam search algorithm necessary for introducing our algorithm. However, PN theory is extensive, so readers are referred to [58] for more information on the properties and variants of PN.

3.1. Petri Nets

PN were introduced by Carl A. Petri in 1960 as part of his doctoral dissertation. Since then, many researchers have studied PN theory and its applications. Additionally, other types of PN, such as timed PN, colored PN, and stochastic PN, have been developed to enhance their modeling power. In this paper, we will refer to the ordinary PN.

An ordinary PN, or simply PN, is a bipartite graph consisting of two sets of nodes: $P=\{p_1,p_2,\dots ,p_r\}$, and transitions $T=\{t_1,t_2,\dots ,t_s\}$, such that $P\cap T=\varnothing$. These nodes are connected by a set of directed arcs $F\subseteq \left\{\right(P\times T)\cup (T\times P\left)\right\}$. Places may represent activities or conditions in systems modeling, whereas transitions represent the occurrence of a condition or the execution of an activity [58]. Figure 1 shows a typical example of a PN model, where circles represent places and rectangles represent transitions, respectively. Each arc carries a number or arc weight, defined by the function $W:F\to N^{+}$ , which helps simulate the system’s dynamics. Moreover, a particular state of the system is represented by a distribution of tokens (black dots) within the places, resulting in the marking of the PN. Formally, a marking is defined by the function $M:P\to N$. A marked place means that the corresponding activity or condition is being executed.

The dynamics or evolution of PN markings involves removing and adding tokens when a firing occurs. Let ^•$t$ and $t^{•}$ denote the sets of input and output places of transition t. If $M\left(p\right)\ge W(p,t)$ for each $p{\in }^{•}t$ at some marking M, then the transition t is enabled and can be fired. When transition t is fired, $W(p,t)$ tokens are removed for each $p{\in }^{•}t$, and $W(t,p)$ tokens are added to each output place $p\in t^{•}$. The dynamics of a PN can also be described by matrix operations using Equation (6), known as the state equation where ${\mathbf{M}}_k$ is an $r\times 1$ vector representing the current state or marking of the PN, and ${\mathbf{M}}_0$ is the initial marking or state of the PN. $\mathbf{I}$ is the incidence matrix of size $r\times s$, given by Equation (7) where ${\mathbf{I}}^{+}$ and ${\mathbf{I}}^{-}$ are the input and output incidence matrices, respectivily.

$${\mathbf{M}}_{k+1}={\mathbf{M}}_k+\mathbf{Iu}$$

(6)

$$\mathbf{I}={\mathbf{I}}^{+}-{\mathbf{I}}^{-}$$

(7)

Each entry of ${\mathbf{I}}^{+}$ corresponds to the weights $W(t,p)$, while in ${\mathbf{I}}^{-}$ the entries correspond to the weights $W(p,t)$, as illustrated in Figure 2. Finally, $\mathbf{u}$ is the firing vector of size $s\times 1$, where one entry is 1 (indicating the transition to be fired at marking ${\mathbf{M}}_k$), and the remaining entries are zero. Typically, from a set of enabled transitions at a marking $\mathbf{M}$, one transition is fired at a time to explore the possible states that can be reached from that marking. This is represented by the notation $\mathbf{M}\stackrel{t}{\to }{\mathbf{M}}^{\prime }$, where t is the fired transition and ${\mathbf{M}}^{\prime }$ is the marking reached from $\mathbf{M}$.

Given a PN with an initial marking ${\mathbf{M}}_0$, a firing path $\sigma$ is a sequence of transitions fired from ${\mathbf{M}}_0$ to marking ${\mathbf{M}}^{\prime }$, denoted as $\mathbf{M}\stackrel{\sigma }{\to }{\mathbf{M}}^{\prime }$. When this occurs, ${\mathbf{M}}^{\prime }$ is said to be reachable from ${\mathbf{M}}_0$. The set of all possible reachable markings from ${\mathbf{M}}_0$ can be represented by a reachability tree, provided that this set is finite. An example of a reachability tree, corresponding to the PN in Figure 1, is shown in Figure 3. Taking the marking ${\mathbf{M}}_0$ as the root of the reachability tree, a new branch is created for each transition fired from ${\mathbf{M}}_0$. After exhausting all possible firings from ${\mathbf{M}}_0$, a new level of markings is created, such as ${\mathbf{M}}_1$ in Figure 3. At the new level, for each marking, new branches are created using the same procedure. When all the possible markings have been obtained, a new level is added to the tree. The process is repeated for each of the markings at every new level until all the reachable markings from ${\mathbf{M}}_0$ have been obtained. However, when no transition can be fired from a specified marking, it is considered a dead marking. If the evolution of markings is stopped by dead markings without reaching a desired final marking, the PN is blocked.

It should be noted that although the set of reachable markings from ${\mathbf{M}}_0$ can be finite, the number of these markings could be so large that it becomes impossible to enumerate or represent them using a reachability tree. This issue is known as the state explosion of PN.

3.2. Filtered Beam Search Algorithm

The beam search algorithm was first used as an artificial intelligence technique for speech recognition. It has since been applied to combinatorial optimization problems, such as scheduling problems [59], the travelling salesman problem [60,61,62], and cutting and packing problems [63,64,65]. This method is efficient because solutions to these problems can be constructed sequentially by appending the next operation, city, or item to the partial solution based on the cost it adds. Thus, the construction of the solution can be viewed as a decision tree. At each level, the corresponding nodes represent partial solutions with given costs, while at the last level, the nodes represent complete solutions. One of these complete solutions is then selected as the best or optimal solution. The beam search algorithm performs a breadth-first search followed by a depth-first search. At each level of the tree, the nodes are branched, and those that are not promising for further branching based on a cost function evaluation are pruned. The number of nodes $\beta$ that are branched at each level are the beams. Once these nodes have been selected, the others are discarded permanently. Since the beam search has no backtracking, it is possible to discard good solutions, but it saves much computational effort. In [66], a modification to the beam search was introduced known as the filtered beam search algorithm. As the name suggests, at each level, a local evaluation is carried out over the nodes that share the same parent node. Only those nodes that satisfy a local evaluation condition are considered beam nodes, thereby filtering the nodes. Once all the filtered nodes have been obtained from the parent nodes, only those that satisfy a global evaluation condition become the new beam nodes or parent nodes for the next level.

The number of filtered nodes $\alpha$ (filter width) and the number of beam nodes $\beta$ (beam width) are user-defined. However, it is known that larger values require more computational effort to search the tree, while smaller values reduce the likelihood of finding the best or optimal solution. A representation of the filtered beam search algorithm is shown in Figure 4 for a solution with a sequence of n nodes. Sometimes, it is desired to maintain diversity in the search by setting $\alpha <\beta$, which guarantees that at least one beam node is selected from a different parent node [64].

4. PN-Based Algorithm for Solving the 1D-CSP

In this section, we first provide details about the PN model used to construct solutions for a given instance of the 1D-CSP. Next, we explain how the beam search algorithm, based on the reachability tree, is used to find the best solution.

4.1. PN Model for the 1D-CSP

As mentioned above, we use ordinary PN to model a given instance of the 1D-CSP. Figure 5 illustrates the general PN model for a given instance of the 1D-CSP, with m different items and stock length L. The construction of a solution for a 1D-CSP instance through the PN model follows the next-fit heuristic principle. This means that items are added sequentially to the solution as long as the total length of the items, or pattern length, does not exceed the length of the stock. Otherwise, the next item is put into a new stock, and the process is repeated until all items have been added to the solution. In this way, there will be patterns with or without waste, and the total waste of the solution will be the sum of the waste yielded from each pattern. The PN model allows the formation of patterns with and without waste, enabling the beam search algorithm to find the solution with the minimum waste and the minimum number of stocks with waste.

According to the general PN model, there will be m places for each item type labeled $p_2,p_3,\dots p_{m+1}$, as shown in Figure 5. We refer to these places as item places. The number of tokens inside the item places represents the number of orders of each item, denoted $d_1,d_2,\dots ,d_m$. The stock is represented by a single place labeled $p_1$, and the number of tokens inside this place is equal to the stock length, L. We refer to this place as the stock place. Item places, as well as the stock place, are input places for transitions $t_1,t_2,\dots ,t_m$, referred to as item transitions. Each item transition consumes a token from its corresponding item place, indicating that an order of that item is being put into a solution. Similarly, each item transition removes $l_i$ tokens from the stock place, which corresponds to the length of the item of the corresponding item place. As items are taken, a pattern is formed, and its length is accumulated in the place labeled $p_{m+2}$, which we name the pattern place. As shown in Figure 5, the pattern place is an output place of each item transition, where each one adds $l_i$ tokens to the pattern place according to the item length it is being added to the solution. Therefore, the number of tokens inside the pattern place varies from 0 to L. However, the minimum pattern length will be $\lambda =min\{l_1,l_2,\dots ,l_m\}$, representing a pattern formed by only one order of the item with the smallest length.

The transitions $t_{m+1},t_{m+2},t_{m+3},\dots ,t_{m+L}$ are the pattern transitions, and there will be L transitions of this type. A pattern transition removes tokens from the pattern place according to the corresponding pattern length that has been formed. For example, $t_{m+1}$ removes L tokens, $t_{m+2}$ removes $L-1$ tokens, and so on, such that $t_{m+L}$ removes $\lambda$ tokens. Since a pattern includes the lengths of the items as well as the unused material of the stock (if any), pattern transitions $t_{m+2},t_{m+3},\dots ,t_{m+L}$ are output transitions of the stock place to complete the pattern formation. Thus, $t_{m+2}$ removes one token, $t_{m+3}$ removes two tokens, and so on, such that $t_{m+L}$ removes $L-\lambda$ tokens from the stock place. Transition $t_{m+1}$ is not an output transition of the stock place since it corresponds to a pattern without waste. At the same time, the pattern transitions are input transitions of the stock place so that once a pattern has been formed and included in a solution, another stock must be available to continue the process of forming the solution. Each pattern transition adds L tokens to the stock place as shown in Figure 5. So far, the PN model can construct any possible solution by starting with any item, but we included three more places to enhance the model. These places are $p_{m+3}$, $p_{m+4}$ and $p_{m+5}$, where $p_{m+3}$ records the amount of waste, $p_{m+4}$ helps to count the number of stocks without waste, and $p_{m+5}$ counts the number of stocks with waste. Places $p_{m+3}$ and $p_{m+5}$ are output places of the pattern transitions. Each pattern transition adds one token to $p_{m+5}$ while $p_{m+3}$ receives $1,2,\dots ,L-\lambda$ tokens from transitions $t_{m+2},t_{m+3},\dots ,t_{m+L}$, respectively. Place $p_{m+4}$ is an output place of transition $t_{m+1}$, which adds one token each time there is a pattern without waste. Based on the above descriptions, Figure 6 shows the general form of the incidence matrices of the PN model for an instance of the 1D-CSP.

4.2. Marking Evolution and Solution Construction

To obtain a solution for a given instance of the 1D-CSP with the PN model described in the previous section, referred to as PN-CSP, we begin with the initial marking as shown in Figure 5 where there is an unused stock (place $p_1$) and m items with their corresponding number of orders (places $p_2$ to $p_{m+1}$). Hence, the general form of the initial marking of the PN-CSP is:

$${\mathbf{M}}_0={[L{d}_1{d}_2\cdots d_{m+1}0000]}^T.$$

From the initial marking, the enabled transitions are $t_1,t_2,\dots ,t_m$ so that any of the items could be the first in the solution. The solution is constructed in a linear sequence of items arranged from left to right. Once an enabled transition is fired, the pattern place accumulates tokens representing the length of the item placed first in the solution. Suppose transition $t_2$ is fired; then the resulting marking is:

$${\mathbf{M}}_1={[L-l_2{d}_1{d}_2-l_2\cdots d_{m+1}{l}_{2}000]}^T.$$

While the number of tokens in $p_{m+2}$ is not equal to or greater than any of the weights of the arcs $L,L-1,L-2,\dots \lambda$, transitions $t_1,t_2,\dots ,t_m$ can be fired one at a time because they will be the only ones enabled. This means the pattern formation continues as long as there is space in the current stock. However, if a determined pattern length less than the stock length is reached, such as $L-1,L-2$, or $\lambda$, the corresponding transition from $t_{m+2}$ to $t_{m+L}$ is enabled, considering there will be $1,2,$ or $L-\lambda$ tokens in place $p_1$, representing the waste of the stock. If one of the transitions $t_{m+2},\dots ,t_{m+L}$ is selected to be fired, the pattern formation stops to renew the stock for a new pattern. For example, consider the marking:

$${\mathbf{M}}_k={[2ab\cdots cL-2000]}^T$$

where $a,b,\dots ,c$ are the corresponding remaining number of items’ orders after $k-1$ items have been included in the first pattern of the solution, and $M_k\left(p_{m+2}\right)=L-2$. Then $t_{m+3}$ is enabled and can be fired, yielding the marking:

$${\mathbf{M}}_{k+1}={[Lab\cdots c0201]}^T.$$

Here, there is a new empty stock of length L, and one stock with waste ($M_{k+1}\left(p_{m+5}\right)=1$) has been generated as part of the solution with a waste of 2 ($M_{k+1}\left(p_{m+3}\right)=2$). With the new empty stock, the construction of the solution continues with the remaining items’ orders. When a pattern transition is fired, it means that a cutting pattern has been formed, adding a new stock with or without waste to the solution. These markings that enable and fire a pattern transition are called pattern markings. On the other hand, if an item transition is fired, it means that an item has been added to the current pattern in construction. These markings that enable and fire an item transition are called item markings. The marking evolution is performed as described, forming patterns and renewing the stocks when needed until all of the items have been included in a solution.

When the number of tokens in places $p_2,p_3,\dots ,p_{m+1}$ is zero, the solution has been completed, and the marking ${\mathbf{M}}_f={[L00\cdots 00\omega \nu \eta ]}^T$ is considered the final marking, where $\omega$ is the final total waste, $\nu$ is the final number of stocks without waste, and $\eta$ is the final number of stocks with waste. Therefore, the solution can be obtained with the firing sequence $\sigma$ obtained from ${\mathbf{M}}_0$ to ${\mathbf{M}}_f$, that is ${\mathbf{M}}_0\stackrel{\sigma }{\to }{\mathbf{M}}_f$. A final firing sequence would have the following form:

$$\sigma =t_2,t_m,\dots ,t_1,t_{m+3},t_1,t_3,\dots ,t_m,t_{m+1},\cdots ,t_{m+L},t_m,t_1,\dots ,t_3,t_{m+2}.$$

To decode $\sigma$ into a sequence of item lengths, the pattern transitions $t_{m+3}$, $t_{m+1}$, $t_{m+L}$ and $t_{m+2}$ are dropped, and a new sequence is made up with the remaining transitions. This new sequence is:

$$\varsigma =t_2,t_m,\dots ,t_1,t_1,t_3,\dots ,t_m,\cdots ,t_m,t_1,\dots ,t_3.$$

The solution is formed with either the weights $W(t_i,p_{m+2})$ or $W(p_1,t_i)$ for $1\le i\le m$ for each of the transitions in $\varsigma$, such that these weights correspond to the lengths of the items. Therefore, for the final firing sequence $\sigma$, the solution would have the following form:

$$s=l_2,l_m,\dots ,l_1,l_1,l_3,\dots ,l_m,\cdots ,l_m,l_1,\dots ,l_3.$$

5. Filtered Beam Search Algorithm Implementation

Consider the reachability tree construction of the PN-CSP with initial marking ${\mathbf{M}}_0={[L{d}_1{d}_2\cdots d_{m+1}0000]}^T$ for a given instance of the 1D-CSP. Starting from this initial marking, all potential solutions are derived by investigating firing sequences leading to the final marking ${\mathbf{M}}_f={[L00\cdots 00\omega \nu \eta ]}^T$. However, exhaustive exploration of all sequences is computationally intensive. To solve this issue, we use the FBS, which reduces explored markings as we progress in the reachability tree. FBS does not perform backtracking like traditional any-time search methods, optimizing computational efficiency. The effectiveness of FBS depends on parameters $\alpha$ and $\beta$, chosen to balance exploration and exploitation. Algorithm 1 outlines our FBS implementation, starting by generating all child nodes from firing transitions $t_1,t_2,\dots ,t_m$ at the initial marking or the root node (level 0 or $L0$), resulting in m markings at level 1 ($L1$), as shown in Figure 7. In our FBS, each marking at L1 is neither locally nor globally evaluated, treating each node as a potential beam node, contrary to [59]. We made this decision because, by the local evaluation, the algorithm will tend to select the solutions with the largest items placed first, minimizing waste in the initial stock. Preliminary experiments showed this approach yields better solutions compared to selecting only $\beta$ nodes at $L1$. Thus, starting from L2, we apply both local and global evaluations, setting $\alpha \le \beta$ to ensure solution diversity.

In most beam search algorithm implementations, node costs at any node are computed using a function that considers cumulative costs up to the node and an estimation of the remaining solution’s cost. The effectiveness of these algorithms heavily depends on choosing the right estimation cost function. Usually, in the FBS algorithm, local evaluation uses cumulative costs, while global evaluation considers both cumulative and estimated costs. Our FBS implementation employs the bi-objective function (1) as a fitness function to perform both local and global evaluations. This function calculates the complete solution, capturing the number of stocks used, the total waste generated, and the number of stocks with waste that have been generated in the solution.

As we progress in constructing the reachability tree, each node or marking represents only a partial solution, thus its corresponding cost is also partial. Despite not having a complete solution, partial solutions provide insights into the number of stocks used, those with or without waste, and the amount of unused material or potential waste on the current stock. For example, at the initial marking or root node, no stock has been used, resulting in no stocks with or without waste, and an unused material length of L. At level L1, item markings represent solutions where one order of each item has been taken, but no stock has been fully utilized with or without waste. Eventually, as the reachability tree progresses, pattern markings are reached at some level, indicating that a cutting pattern has been obtained. At the first pattern marking, a stock has been used, and the waste amount for the stock is known. Subsequent item markings’ costs must consider the used stock and its waste. It can be seen from the bi-objective function (1) that the first term sums the partial wastage of each used stock, while the second sums the number of stocks with waste. We compute two cumulative costs accordingly: one for waste using Equation (8) and another for the number of stocks with waste using Equation (9). Both are employed in Equation (10) for the global evaluation of nodes. These two cumulative costs are computed upon reaching a pattern marking, while Equation (11) is utilized for local evaluations, considering the unused material of the current stock in each solution.

$$W_w=W_w+\sqrt{{\displaystyle \frac{w_j}{L}}}$$

(8)

$$W_{sw}={\displaystyle \frac{M\left(p_{m+5}\right)}{M\left(p_{m+5}\right)+M\left(p_{m+4}\right)}}$$

(9)

$$F_{global}={\displaystyle \frac{W_w+W_{sw}}{M\left(p_{m+5}\right)+M\left(p_{m+4}\right)+1}}$$

(10)

$$\begin{array}{cc}\hfill F_{local}& =\left\{\begin{array}{cc}{\displaystyle \frac{1}{M\left(p_{m+5}\right)+M\left(p_{m+4}\right)+1}}\left({\displaystyle \frac{M\left(p_1\right)}{L}}+W_{sw}\right)\hfill & \mathrm{if}M\mathrm{is}\mathrm{an}\mathrm{item}\mathrm{marking}.\hfill \\ & \\ {\displaystyle \frac{1}{M\left(p_{m+5}\right)+M\left(p_{m+4}\right)+1}}\left({\displaystyle \frac{w_j}{L}}+W_{sw}\right)\hfill & \mathrm{if}M\mathrm{is}\mathrm{a}\mathrm{pattern}\mathrm{marking}.\hfill \end{array}\right.\hfill \end{array}$$

(11)

In Equations (8) and (11), $w_j$ is obtained from the incidence matrix of PN-CSP, where $w_j=\mathbf{I}(p_{m+3},t)$ for $t\in \{t_{m+2},t_{m+3},\dots ,t_{m+L}\}$. This is because when a pattern transition is fired, the marking of place $p_1$ is immediately updated to L, indicating that the amount of waste in the already used stock is lost. However, this amount of waste can be determined from the incidence matrix. The cumulative costs $W_w$ and $W_{sw}$ are set to zero for the item markings between the initial marking and the first pattern marking. These costs are computed and updated at each reached pattern marking. Consequently, the outcomes of $W_w$ and $W_{sw}$ are retained for the item markings between two subsequent pattern markings. Once the final marking is reached, the last $F_{global}$ outcome is set as the best solution $F_{best}$. To determine whether the final marking has been reached, at each level, we compare the markings obtained with the desired final marking ${\mathbf{M}}_f\left(p_1\right)$. However, we only consider the marking of places $p_1,p_2,\dots ,p_{m+2}$. Specifically, we compare the marking $\mathbf{M}=[M\left(p_1\right)M\left(p_2\right)M\left(p_3\right)\cdots M\left(p_{m+1}\right)M\left(p_{m+2}\right)***]$, obtained at some level, with the final marking ${\mathbf{M}}_f=[L00\cdots 00***]$. The ’*` means that the marking of the corresponding places is not considered for the comparison. This is because we do not know in advance how much waste would be generated by the best solution or the number of stocks with and without waste. However, we know that at the final marking, there are no more item orders to place, and the last stock is generated when the last item order is taken for the final cutting pattern.

Observe from Algorithm 1, lines 17 and 29, that when performing global and local evaluations, we do not necessarily select $\beta$ or $\alpha$ number of markings, respectively. Instead, we select the minimum value between $\beta$ (or $\alpha$) and the number of markings obtained globally (or locally). This is because, in a reachability tree, a marking could be generated by firing different transitions, resulting in fewer markings than the $\beta$ markings needed at a level or fewer than the $\alpha$ markings needed from a parent marking or beam node. In some beam search algorithm implementations, if the required $\beta$ nodes are not met for some level, then more parent nodes at the previous level are explored to meet the required $\beta$ nodes. However, we do not do this in our FBS implementation, as it is possible to encounter dead markings. This might make it impossible to complete the required number of nodes. Dead markings are related to the liveness property of a PN and this property is closely related to the presence of siphons [67]. It is known that a PN with siphons is likely to be non-live, as once the places of a siphon become unmarked at some marking, they remain unmarked in all subsequent markings. We performed a structural analysis for an instance with ten different items and a total of 20 orders, using an open-source PN modeling platform called PIPE2 version 4.3.0 [68] to determine the presence of siphons. The results are shown in Figure 8, indicating that place $p_2$ ($p_1$ according to the enumeration of the PIPE2 platform) is a minimal siphon. Extending this result for the general PN-CSP model with the initial marking in the general form, and considering the dynamics described in the above paragraphs, places $p_1,p_2,\dots ,p_m$ could be minimal siphons and could become unmarked at some marking. The PN-CSP would become deadlocked if all the places $p_1,p_2,\dots ,p_m$ are unmarked, although this situation indicates that the desired final marking has been reached. However, to prevent our algorithm from getting stuck due to incomplete $\beta$ nodes at some levels and potentially failing to reach the desired final marking, we decided to use the minimum number of nodes as outlined in Algorithm 1.

Algorithm 1 Pseudocode of algorithm FBS-PN-CSP

Input: Number of different item types (m), stock length (L), item lengths $(l_1,l_2,\dots ,l_m)$, and item orders $(d_1,d_2,\dots ,d_m)$.  Output: $F_{best}$ 1: With the instance data, generate the incidence matrices ${\mathbf{I}}^{+}$, ${\mathbf{I}}^{-}$ and obtain the incidence matrix $\mathbf{I}$. 2: Generate the initial marking as ${\mathbf{M}}_0=[L{d}_1{d}_2\cdots d_{m+1}0000]$. 3: Set the beamwidth $\beta$ and the filterwidth $\alpha$. 4: $F_{best}\leftarrow \infty$ 5: $W_w\leftarrow 0$ 6: $W_{sw}\leftarrow 0$ 7: Generate level $L0$ of the reachability tree with the marking ${\mathbf{M}}_0$ as the root node such that $Currentlevel=\left\{{\mathbf{M}}_0\right\}$. 8: $level\leftarrow 0$ 9: $Bestsolution\leftarrow False$ 10: while $BestsolutionisFalse$ do 11:    for each marking in $Currentlevel$ do 12:      Get $W_w$, $W_{sw}$ and $F_{global}$. 13:    end for 14:    if level > 1 then 15:      Kept the $\mu =\min (\beta ,|Currentlevel\left|\right)$ best markings in $Currentlevel$ with the lowest $F_{global}$ cost and form the set $\mathcal{B}$ with these markings. Eliminate the remaining markings from $Currentlevel$. 16:    else 17:      Kept the $\mu =\left|Currentlevel\right|$ markings in $Currentlevel$ and form the set $\mathcal{B}$ with these markings. 18:    end if 19:    $Newlevel=\varnothing$ 20:    for each marking $\mathbf{M}$ in $\mathcal{B}$ do 21:      Determine the reachable markings ${\mathbf{M}}^{\prime }$ from $\mathbf{M}$ with Equation (6) and form the set ${\mathcal{M}}_r$ with those reachable markings, where $r=1,2,\dots ,\left|\mathcal{B}\right|$. 22:      if ${\mathcal{M}}_r\ne \varnothing$ then 23:         Obtain $W_w$, $W_{sw}$, $F_{local}$ and $F_{global}$ for each ${\mathbf{M}}^{\prime }\in {\mathcal{M}}_r$. 24:         if level >0 then 25:           Check whether any marking ${\mathbf{M}}^{\prime }\in {\mathcal{M}}_r$ has already been generated previously. If it has, let ${\mathbf{M}}_{prev}^{\prime }\in {\mathcal{M}}_1\cup {\mathcal{M}}_2\cup \cdots \cup {\mathcal{M}}_{r-1}$ be the previously generated marking that is equal to ${\mathbf{M}}^{\prime }\in {\mathcal{M}}_r$, and then eliminate from the reachability tree the one with the maximum $F_{global}$ cost. 26:           if ${\mathcal{M}}_r\ne \varnothing$ then 27:              Kept the $\tau =\min (\alpha ,|{\mathcal{M}}_r\left|\right)$ markings with the lowest $F_{local}$ cost and eliminate the remaining ones from ${\mathcal{M}}_r$. 28:           end if 29:         else 30:           Keep the markings in ${\mathcal{M}}_r$. 31:         end if 32:         Make $Newlevel\cup {\mathcal{M}}_r$ 33:      end if 34:    end for 35:    Check whether the final marking ${\mathbf{M}}_f=[L00\cdots 00***]$ has been reached. 36:    if ${\mathbf{M}}_f$ has been reached then 37:      $F_{best}\leftarrow F_{global}$ 38:      $Bestsolution\leftarrow True$ 39:    else 40:      $level\leftarrow level+1$ 41:      $Currentlevel=Newlevel$ 42:    end if 43: end while

6. Experiments

Although we identified a wide variety of datasets in the literature, most studies do not use common datasets or performance measures to compare the algorithms. Additionally, implementing the published algorithms correctly can be challenging, and often, the code is not shared. Nevertheless, we identified two algorithms that share some similarities with our algorithm and are not based on populations. These algorithms were tested on five datasets that we found suitable for comparison with our algorithm. The first algorithm is the Least Lost Algorithm (LLA), introduced in [54], which uses a dataset of ten instances with a single stock size. We refer to this dataset as hinterding because, to the best of our knowledge, it was first used in [35]. The LLA is executed in two phases. In the first one, the item lengths are arranged in decreasing order while in the second, the second half of the arrangement is put before the first half. In both phases, the items are taken from the arrangement to construct cutting patterns by beginning with those that have no waste. If no additional cutting patterns without waste can be found, the search continues with patterns that have one unit of waste. This process continues by increasing the allowed waste by one unit until all the item orders have been included in a cutting pattern. If the number of stocks used is above the optimal theoretical value by more than 5%, the algorithm performs the second phase. The feature of increasing the allowed waste by one unit is represented in our PN-CSP model with the pattern transitions. The LLA was compared to the EP algorithm shown in [8] over the same dataset and the LLA yielded better results.

The other algorithm we compare to the FBS-PN-CSP is the Generate & Solve algorithm (G&S) presented in [55]. This algorithm is also executed in two phases. The first one reduces the problem size using a tree method that constructs all possible cutting patterns. Feasible cutting patterns are then selected and sent to the second phase, where an integer linear programming (ILP) method is implemented to obtain the best solution for the reduced set of solutions. In the first phase, a new reduced set of solutions is generated by considering an aging mechanism for previously used solutions, which can be re-selected as long as they have not reached a specified age limit. These new solutions are then sent to the second phase to determine the best solution among them, replacing the older solution if the new one is better. This process continues until a predefined time limit is reached, and the last solution found is returned as the best solution for the problem. The G&S algorithm was compared to the column generation technique. Both algorithms performed similarly for small instances, but the G&S algorithm was more effective for large instances. From the datasets used in [55], we selected the ones called hard28, hard10, and wae_gau1, available at [69]. Additionally, we included the dataset falkenauer_U from the same repository, which consists of four groups named binpack1, binpack2, binpack3, and binpack4. Table 1 provides details about the datasets. The third column lists the minimum and maximum number of item types among the instances of each dataset, while the fourth column shows the number of different item types, along with the standard deviation in rounded brackets. Both the hinterding dataset and the hard10 dataset have four different item types, but the variability in the former is greater than in the latter.

The comparisons among the three algorithms were made with respect to the number of stocks used and its optimal theoretical value, which can be obtained using Equation 12. The experiments were conducted on a Quad-Core Intel^(R) Core^(TM) i5 computer with a 3.4 GHz processor and 8 GB RAM. All three algorithms were coded in Python 3.11.1. For the second phase of the G&S algorithm, we used the Scipy package with the linprog library. We also set the same parameter values as in [55], that is, specifically: a time limit of 600 seconds for the G&S algorithm, a mutation rate of 10%, a maximum age of 2, and a maximum number of cutting patterns for the tree construction of 1,000,000. Additionally, we also executed the G&S algorithm 10 times for each instance. However, we did not set a time limit in the second phase of the G&S algorithm, as the solutions from the ILP solver were obtained in an acceptable time.

$$lb=⌈{\displaystyle \frac{{\sum }_i^m{l}_i{d}_i}{L}}⌉$$

(12)

One of the issues with beam search algorithms is the selection of the beam width $\beta$. The smaller the $\beta$ value, the faster the search, but this comes at the loss of optimality due to aggressive pruning. We conducted preliminary experiments to determine suitable $\beta$ and $\alpha$ values, considering the condition $\alpha \le \beta$. These experiments were carried out with the hinterding dataset. In all ten instances of this dataset, we observed that the best solutions were found with a beam width between 5 and 25; with higher values, the solution improvement was negligible. Thus, for all the experiments with the FBS-PN-CSP algorithm, we used the interval $5\le \beta \le 35$ for the beam width, at steps 5, 15, 25, and 35, with $2\le \alpha \le \beta$ for each value of $\beta$. However, for the hinterding dataset, we additionally implemented a $\beta$ value of 45, and for the hard10 dataset, we only could implement a beam width of $\beta =5$. With $\alpha =1$, no good solution was found in any instance.

7. Results

Table 2, Table 3, Table 4, Table 5, Table 6, Table 7, Table 8 and Table 9 show the best results obtained with each algorithm for the number of stocks used and the time taken. Additionally, we provide the $lb$ value for all the instances and the percentage above this value (%$>lb$) obtained by each algorithm. The FBS-PN-CSP, G&S, and LLA algorithms were compared only with the hinterding and falkenauer_U datasets since the LLA algorithm ran out of memory for the other datasets. The minimum value of stocks used among the compared algorithms is highlighted in boldface. For the G&S algorithm, we provide the best solution found within the 10 executions, and likewise, for the FBS-PN-CSP algorithm, we include the best combinations of beam width and filter width where the best stock result was found.

For the hinterding dataset, the three algorithms obtained the same optimal theoretical value for the 1a-5a instances and the 7a instance. The FBS-PN-CSP and the G&S algorithms also obtained the same best result for the 6a instance. Only the G&S had the optimal theoretical value for the 8a-10a instances. The FBS-PN-CSP algorithm had better results than the LLA algorithm for the 6a, 8a, and 9a instances, while for the 10a instance, both had the same result. Although the FBS-PN-CSP algorithm was not better than the G&S algorithm for instances 8a-10a, it obtained a percentage above the optimal theoretical value of less than 1%. The FBS-PN-CSP algorithm was more efficient than the G&S algorithm in 6 out of 7 instances where they are equal, namely 2a-7a. However, the LLA algorithm had better times in all the instances where it had the same results as the other two algorithms. The beam width of the FBS-PN-CSP algorithm seems to increase as the number of items and the stock length increase. The highest beam width and filter width values were obtained in the 10a instance, which has the higher number of items and the longest stock.

The results of the FBS-PN-CSP and G&S algorithms are shown in Table 3 for the hard28 dataset. Both algorithms obtained the same best results except for the BPP900 and BPP781 instances, where the FBS-PN-CSP had the better results. Nevertheless, neither obtained the optimal theoretical value. The FBS-PN-CSP was more efficient in 15 out of 27 instances where both algorithms obtained the same result. As for the beam width of the FBS-PN-CSP, this was the minimum value of the tested values, that is, $\beta =5$ for all the instances. For the hard10 instances, the G&S algorithm was better than the FBS-PN-CSP algorithm as shown in Table 4. However, the G&S algorithm only had the optimal theoretical value in one instance, namely HARD9. Even with the smaller value of the beam width, the time required for the FBS-PN-CSP was far higher than the time of the G&S algorithm. For this reason, we could not implement larger values for the beam width since the computational cost was prohibited, although it did not run out of memory. Despite the results for stock used by the FBS-PN-CSP algorithm, it was not far from those of the G&S and could be at least equal if higher values of the beam width can be implemented.

In 12 out of 17 instances of the waegau1 dataset, the FBS-PN-CSP algorithm matched the best results of the G&S algorithm (Table 5), and in 9 of these 12 instances, it had better times. In 4 out of the 17 instances, the FBS-PN-CSP had better results, and only in one instance, the G&S was better than the FBS-PN-CSP. In 7 out of the 17 instances, the FBS-PN-CSP obtained the optimal theoretical value, while the G&S did so in only four instances. Regarding beam width, four instances required the higher tested values, namely 25 and 35, but with a small filter width.

Table 6 shows the best stock results for the binpack1 instances of the falkenauer_U dataset, where the FBS-PN-CSP, G&S, and LLA algorithms are compared. The FBS-PN-CSP and the G&S algorithms obtained the optimal theoretical values of each instance, but the FBS-PN-CSP had the best times. The LLA algorithm achieved the optimal theoretical values in 12 out of 20 instances and had better times than the FBS-PN-CSP algorithm. Only in two instances did the FBS-PN-CSP use the higher beam width values (25 and 35), but in most instances, the value was 5, which is the smallest value tested. For the binpack2 instances, the FBS-PN-CSP had the same results as the G&S algorithm in 19 instances (Table 7), and these were equal to the corresponding optimal theoretical values. All the results of the G&S algorithm were equal to the optimal theoretical value. The LLA algorithm had the optimal theoretical values in half of the instances and the best times among the three algorithms. The FBS-PN-CSP algorithm had better times than the G&S in ten of the instances where they had the same stock result. The beam width was mostly 15 throughout the binpack2 instances. The binpack3 results (Table 8) show that the G&S algorithm obtained the optimal theoretical values in all instances, and the FBS-PN-CSP algorithm did so in 17 instances. Likewise, the LLA algorithm achieved this in only 11 instances, but with better times than the FBS-PN-CSP and G&S algorithms. Only in the instance u500_14 did the FBS-PN-CSP obtain a better time than the G&S algorithm. The values of the beam width were 15 and 25 in most instances, and the filter value was higher than in the previous groups of instances of the falkenauer_U dataset. Similar results were obtained for the binpack4 instances shown in Table 9, where the G&S algorithm outperformed the other two algorithms. In all instances, this algorithm obtained the optimal theoretical value, the FBS-PN-CSP did so in 14 instances, and the LLA algorithm in 7 instances, but only with better times than the FBS-PN-CSP. The beam width was mostly higher, namely 25 throughout the instances, as was the filter width. Although in the instances of binpack2, binpack3, and binpack4 the FBS-PN-CSP was not better than the G&S in terms of the number of stocks used, the percentage above the optimal theoretical values was less than 1% with a difference of one stock in all cases.

8. Discussion

It can be observed from the results shown in the previous section that the FBS-PN-CSP algorithm is effective and efficient in most of the instances it was tested in comparison with the other two algorithms, except for the hard10 dataset. It must be taken into account that the G&S algorithm had a time limit for its execution of 600 seconds. The FBS-PN-CSP best results were obtained below such time in only one execution. Additionally, the average time from 10 executions of the G&S was longer than the FBS-PN-CSP algorithm and the LLA algorithm in most instances. On the other hand, the LLA algorithm had the best times compared to the FBS-PN-CSP and G&S algorithms, but the number of stocks used was of lower quality. We have observed that the computational time for the FBS-PN-CSP for longer stocks tends to be higher due to the PN model structure. As the stock gets longer, the number of pattern transitions increases considerably. This observation can be noted for the hard10 dataset, which has the longest stock among the datasets tested. Moreover, the beam width tends to increase as well as the computational time, as both the number of item types and the number of orders increase. Except for the hard10 dataset, in those instances where the FBS-PN-CSP algorithm did not obtain the optimal theoretical value, the solutions were just above this value by one stock.

9. Conclusions

In this paper, an algorithm based on the modeling power of Petri nets was presented for solving the one-dimensional cutting stock problem (1D-CSP). To the best of our knowledge, Petri nets had not been used to address the 1D-CSP. Firstly, we provided the Petri net general model to construct the solutions of a given instance of the 1D-CSP. Each solution consists of a series of cutting patterns with or without waste. Then, the filtered beam search algorithm is used to search for the best solution throughout the exploration of the nodes of the reachability tree. The algorithm, called FBS-PN-CSP, was compared to the Generate & Solve and the Least Lost algorithms over five datasets with a total of 145 instances and shown to be effective and efficient in most of the datasets. Therefore, either with instances of low or high complexity, the FBS-PN-CSP is shown to be competitive. Some enhancements to the Petri model and slight modifications to the algorithm implementation could improve the performance.