# Degree-Constrainedk-Minimum Spanning Tree Problem

Authors: Pablo Adasme; Ali Dehghan Firoozabadi
Journal: Complexity
 (2020)

---

## Abstract

LetGV,E be a simple undirected complete graph with vertex and edge sets V and E, respectively. In this paper, we consider the degree-constrained k-minimum spanning tree (DCkMST) problem which consists of finding a minimum cost subtree of G formed with at least k vertices of V where the degree of each vertex is less than or equal to an integer value d≤k−2. In particular, in this paper, we consider degree values of d∈2,3. Notice that DCkMST generalizes both the classical degree-constrained and k-minimum spanning tree problems simultaneously. In particular, when d=2, it reduces to a k-Hamiltonian path problem. Application domains where DCkMST can be adapted or directly utilized include backbone network structures in telecommunications, facility location, and transportation networks, to name a few. It is easy to see from the literature that the DCkMST problem has not been studied in depth so far. Thus, our main contributions in this paper can be highlighted as follows. We propose three mixed-integer linear programming (MILP) models for the DCkMST problem and derive for each one an equivalent counterpart by using the handshaking lemma. Then, we further propose ant colony optimization (ACO) and variable neighborhood search (VNS) algorithms. Each proposed ACO and VNS method is also compared with another variant of it which is obtained while embedding a Q-learning strategy. We also propose a pure Q-learning algorithm that is competitive with the ACO ones. Finally, we conduct substantial numerical experiments using benchmark input graph instances from TSPLIB and randomly generated ones with uniform and Euclidean distance costs with up to 400 nodes. Our numerical results indicate that the proposed models and algorithms allow obtaining optimal and near-optimal solutions, respectively. Moreover, we report better solutions than CPLEX for the large-size instances. Ultimately, the empirical evidence shows that the proposed Q-learning strategies can bring considerable improvements.


---

## Body

## 1. Introduction

Intelligent communication systems will be mandatorily required in the next decades to provide low-cost connectivity within many application domains involving network structures in the form of ring, tree, and star topologies [1–3], for instance, when designing network structures in telecommunications, facility location, electrical power systems, water and transportation networks, and for networks to be constructed under the Internet of Things (IoT) paradigm [2–10]. A particular example related with wireless sensor networks is due to density control methods [1, 11, 12]. These methods allow reducing energy consumption in sensor networks by means of activation or deactivation mechanisms which mainly consist of putting into sleep mode some of the nodes of the network while ensuring sensing operations, communication, and connectivity [1–3, 12–16].LetGV,E represent a simple undirected complete graph with set of nodes V=1,…,n and set of edges E=1,…,m. In this paper, we consider the degree-constrained k-minimum spanning tree problem (DCkMST), which consists of finding a minimum cost subtree of G formed with at least k≤n vertices of V where each vertex has a degree lower than or equal to d∈2,…,k−2. Notice that, for d=2, DCkMST reduces to find a minimum cost Hamiltonian path with cardinality k. The DCkMST problem generalizes two classical combinatorial optimization problems, namely, the degree-constrained and k-minimum spanning tree problems (resp., DCMST and kMST). Recently, a variant of DCkMST was studied in [17], where the authors assume that there exists a predefined root node that should not satisfy the degree constraint. In this paper, we consider the more general case where all nodes should satisfy this condition. Proofs related with the NP-hardness of DCMST and kMST can be found, respectively, in [18–20]. However, other applications related to the kMST problem can be consulted in [17, 20]. Notice that the degree constraints ensure that no overloaded nodes are present in the network.It is easy to check from the literature that the DCkMST problem has not been addressed in depth so far. Consequently, our main contributions in this paper can be highlighted as follows. First, we propose three mixed-integer linear programming (MILP) models for the DCkMST problem and derive for each one an equivalent formulation using the handshaking lemma [21–23]. More precisely, we propose four compact polynomial formulations and two exponential models. Two of the compact models are formulated based on a Miller–Tucker–Zemlin-constrained approach, whilst the two remaining ones are flow-based models. We solve our exponential models with an exact iterative method adapted from [2, 24]. In order to obtain feasible solutions alternatively, we further propose ant colony optimization and variable neighborhood search (resp., ACO and VNS for short) algorithms for d=2 and d=3, respectively [25–29]. We choose VNS and ACO metaheuristics as they are well-known techniques in the operations research community which have proved to be highly efficient in order to find feasible solutions for many hard combinatorial optimization problems [25–29]. In particular, within each iteration of the VNS algorithm, for each random k−tuple of vertices generated, we obtain degree-constrained spanning trees by using a modified penalty approach proposed in [30]. However, for the VNS algorithm, we further introduce an embedded Q-learning strategy that allows performing a random local search based on the experience of previous solutions found [31]. The Q-learning approach is a simple reinforcement learning technique initially proposed in [31] that allows an agent to interact with its environment in order to learn from the experience of previous occurrences the best action to choose from a set of actions in order to maximize its revenue. The underlying idea in doing so is to construct optimization methods with learning capabilities in order to make them robust, self-adaptive, and independent from decision-makers. Notice that recently there is a growing interest from the operation research community in developing novel machine learning-based optimization methods that allow solving hard combinatorial optimization problems [32]. Consequently, embedding a Q-learning strategy in a random local search algorithm is a novel approach to the literature. Our embedded Q-learning strategy is then compared with a traditional VNS near-far random local search procedure [28]. Another paper dealing with a similar embedded reinforcement learning strategy using VNS is proposed in [33]. More precisely, the authors propose a reactive search VNS method that allows learning which is the order in which different local search heuristics must be applied in order to obtain better solutions based on the experience of previous trials. Their method is applied and tested on the symmetric traveling salesman problem obtaining good results in terms of solution quality and speed. The embedded Q-learning strategy in our VNS approach is different as we use it as a learning mechanism in order to perform a random local search. Besides, it does not require the use of specialized local search methods, and consequently, it can be extended and used straightforwardly for any other combinatorial optimization problems. Similarly, we also compare the proposed ACO algorithm with an adapted version of the generic ant-Q method proposed in [27]. Finally, we propose a pure Q-learning-based algorithm that is competitive with both ACO methods. We report numerical results for Euclidean benchmark instances from [34] and for randomly generated ones with both uniform and Euclidean distance costs with up to 400 nodes.The paper is organized as follows: in Section2, we present some related work, whilst in Section 3, we present the MILP formulations. Then, in Section 4, we present and explain each proposed algorithm. Next, in Section 5, we conduct substantial numerical experiments in order to compare all the proposed models and algorithms. Finally, in Section 6, we give the main conclusions of the paper and provide some insights into future research.


## 2. Related Work

As mentioned in Section1, the DCkMST problem generalizes both the DCMST and kMST problems simultaneously. Consequently, previously related work is mainly focused on these graph combinatorial optimization problems. Notice that our work in this paper corresponds to an extended version of the preliminary work reported in [21]. However, now we generalize the previous formulations presented in [21] and allow each model to obtain feasible solutions with at least k vertices instead of using a unique value of k. Thus, the proposed models in this paper are more accurate with respect to the definition of the kMST problem [20]. In addition, we propose several algorithmic approaches for the DCkMST problem. Finally, we report numerical results for complete graph instances using degree values of d∈2,3. Notice that solving complete graph instances for the DCMST problem while using degree values of d∈2,3 implies solving the hardest type of instances as reported in the literature [35].Regarding the complexity of thekMST problem, it has proved to be NP-hard by reduction from the Steiner tree problem for variable values of k [19, 20]. Consequently, it is hard to obtain an approximation ratio better than 96/95 [36]. In fact, the best-known approximation ratio for this problem is 2 and is reported in [37]. Some recent metaheuristic approaches are due to [38, 39]. In particular, the authors in [39] propose a new hybrid algorithm based on tabu search and ant colony optimization leading to better numerical results than state-of-the-art values reported in the literature. Similarly, in [38], the authors propose a hybrid approach using a memetic algorithm where the genetic operator is based on a dynamic programming method. Numerical results show that their proposed method outperforms several existing algorithms in terms of solution quality and accuracy.On the contrary, the DCMST problem has been studied extensively in the literature including exact and heuristic approaches such as Lagrangian relaxations, approximation algorithms, and branch-and-cut and metaheuristics approaches [22, 30, 35, 40–44]. In particular, the authors in [35] propose a Lagrangian-based heuristic for the DCMST problem that uses a greedy construction heuristic followed by a heuristic improvement procedure. They report extensive computational experiments and indicate that their solving method is competitive with the best heuristics and metaheuristics proposed in the literature. Instances with up to 2000 nodes are reported. Similarly, the authors in [22] propose a branch-and-cut algorithm for this problem dealing with instances with as many as 800 nodes. Concerning metaheuristic approaches, a novel genetic algorithm is proposed in [41]. More precisely, the authors first propose a transformation of the problem into a preference with two objective minimum spanning tree problems. Then, new crossover, mutation, and selection operators together with a local search approach are derived based on the preference of the two objectives. Finally, they show that their algorithm converges with probability one to a globally optimal solution. A more recent VNS approach is presented in [45] where the authors develop ideas that allow the enhancement of the performance of the VNS by guiding the search in the shaking phase while using second-order algorithms and skewed functions. They conduct computational experiments on hard benchmark instances from the literature and improve upon the best-known solutions. Finally, their solutions significantly reduced the gaps with respect to tight lower bounds reported in the literature as well. Other variants in the DCMST problem such as the multiperiod DCMST and min-degree-constrained minimum spanning tree problems are presented in [46, 47].In this paper, our VNS approach for the DCkMST problem performs a traditional near-far random local search procedure as performed in the reduced VNS presented in [28, 29] and also a novel one which consists of embedding a Q-learning strategy [31] for each random generated k-tuple of vertices of V. In both cases, we use a modified penalty heuristic approach proposed in [30] in order to find spanning trees while satisfying the degree constraints. Finally, our ACO algorithms are adapted from [25–27] and incorporate an adaptive mechanism that allows obtaining consecutive k-Hamiltonian paths within each iteration in order to find feasible solutions. Practical applications related to the DCMST include the design of computer and communication networks, electric circuits, design of integrated circuits, road networks, and energy networks, among many others [17, 35]. Notice that, aside from these applications, all the proposed models and algorithms could also be extended and adapted for multiagent systems with varying graph topologies in order to deal with self-organization and congestion control in communication networks [48, 49].


## 3. Mathematical Formulations

In this section, for the sake of clarity, we first present a feasible solution for the DCkMST problem, and then, we present and explain the proposed MILP models.


### 3.1. A Feasible Solution for the DCkMST Problem

In Figure1(a), we depict an input complete graph instance with 10 nodes where the coordinates of each node are randomly generated within a square area of 1 km2, whereas in Figures 1(b) and 1(c), we show feasible solutions obtained for the instance in Figure 1(a) while using degree values of 3 and 2, respectively. For the sake of clarity, in Figures 1(b) and 1(c), we only draw the edges and highlight with green color the nodes which are part of the solutions, i.e., the nodes 1,2,3,5,6,7. Notice that the minimum cost value of an optimal solution in Figure 1(b) should be lower than or equal to the minimum cost value of an optimal solution in Figure 1(c) since the latter is also a feasible solution for the degree value of 3. On the opposite, the solution in Figure 1(b) is not feasible for the degree value of 2.Figure 1
An input complete graph instance with feasible solutions for the DCkMST problem while using degree values of 3 and 2, respectively: (a) complete graph instance composed of 10 nodes with random uniform costs; (b) feasible solution with node degree values of at most 3; (c) feasible solution with node degree values of at most 2.
(a)(b)(c)Proposition 1.
The following statements are true:(1)
Ifk=n and d≤n−2, the DCkMST problem reduces to the DCMST problem.(2)
Ifk<n and d≥k−1, the DCkMST problem reduces to the kMST problem.(3)
Ifk=n and d≥n−1, the problem reduces to the classical minimum spanning tree (MST) problem, which can be solved in polynomial time by a greedy algorithm [50, 51].
From Proposition1, it is easy to see that the DCkMST problem is NP-hard as it generalizes both the DCMST and KMST problems. Hereafter, we present the MILP models.


### 3.2. MILP Models

In order to formulate a first MILP model, letH=V,A be the directed graph obtained from GV,E, where each edge i,j∈E is replaced by the two arcs i,j and j,i∈A. Thus, we have the following proposition.Proposition 2.
ModelP1 allows to obtain an optimal solution for the DCkMST problem with minimum cost:(1)P1:minx,u,z∑i,j∈APijzij,(2)s.t.∑j∈Vxj≥k,(3)∑i,j∈Azij=∑j∈Vxj−1,(4)ui≤∑j∈Vxj,∀i∈V,(5)ui≥xi,∀i∈V,(6)∑i|i,j∈Azij≤xj,∀j∈V,(7)uj−ui−n−1zij−n−3zji≥2−n,∀i,j∈A,(8)∑i|i,j∈Azij+∑i|j,i∈Azji≤dxj,∀j∈V,(9)∑i|i,j∈Azij+∑i|j,i∈Azji≥xj,∀j∈V,(10)zij+zji≤xi,∀i,j∈A,(11)x∈0,1n,z∈0,1n2,u∈0,∞n,where P=Pij for each i,j∈A is defined as an input symmetric matrix denoting the connectivity costs between nodes i,j∈V. We also define the binary variable zij to be equal to one if and only if arc i,j∈A belongs to the solution of the resulting spanning tree and zij=0 otherwise.Proof.
Notice that the total connectivity cost is minimized in (1). Next, notice that constraints (2) and (3) ensure that at least k nodes from set V should belong to the output solution of the problem and that the number of arcs must be equal to the number of nodes minus one, respectively. For this purpose, we define the binary variable xj for each j∈V where xj=1 if and only if node j is in the solution and xj=0 otherwise. Subsequently, the constraints (4)–(7) ensure that the solution does not contain subtours. For this purpose, we also define the nonnegative variable uj for each j∈V. To see how these constraints avoid cycles, let us assume that the subset of nodes a1,a2,a3,…,at (t≤k) is part of the solution tree and that these nodes are connected according to the following sequence a1−a2−a3−,…,−at. Then, the set of arcs a1,a2,a2,a3,…,at−1,at is also part of the solution tree and za1,a2=za2,a3=za3,a4=⋯=zat−1,at=1. Also, by constraint (7), we have the following relation with the uj variables, ua1<ua2<ua3<⋯<uat. Notice that, for any node in this sequence, we cannot have any other incoming arc; otherwise, we would violate constraint (6) which ensures that each node cannot have more than one incoming arc if it belongs to the solution. In particular, if node a1 is the root node, then no incoming arc can be connected to it either since this would imply that uat<ua1, which is a contradiction. Finally, notice that any node that belongs to another branch of the resulting tree cannot be connected to any node of any other sequence of nodes forming another branch as this is prevented by constraint (6) too. Consequently, the resulting digraph has to be acyclic. Notice that the constraints (4) and (5) only apply for the nodes which are part of the solution by restricting the values of each variable uj,j∈V. Also, the degree constraints (8) and (9) ensure that the number of incoming plus outgoing arcs of each node is less than d∈2,…,k−2 where k≥4. Notice that the degree constraints are redundant if a particular node j∈V is not part of the solution, i.e., if xj=0. In fact, this also implies that ∑i|i,j∈Azij+∑i|j,i∈Azji=0. Next, constraint (10) indicates that at most one of the arcs i,j,j,i∈A can belong to the solution if and only if node i∈V belongs to the solution. Constraint (11) is a domain constraint for the decision variables. Finally, it is observed that an undirected subtree graph can be obtained by dropping the direction of each arc in the resulting digraph, as required.Proposition 3.
Constraint (7) is tighter than constraint (35) presented in [2].Proof.
Constraint (35) in [2] is written as(12)uj−ui≥1−n1−zij,∀i,j∈A.
However, these constraints can be equivalently written as(13)uj−ui−n−1zij≥2−n,∀i,j∈A.
Then, we have that(14)uj−ui−n−1zij≥uj−ui−n−1zij−n−3zji≥2−n,∀i,j∈A.
More details related to constraint (7), when used for the traveling salesman problem and extended to other types of vehicle routing problems, can be consulted in [52]. Notice that model P1 is a directed graph formulation constructed based on a Miller–Tucker–Zemlin characterization [2, 3, 13].Proposition 4 (see [23]).
In a finite simple undirected graph, the sum of the degrees of every vertexv∈V is twice the number of edges.Proposition4 is known as the handshaking lemma, and it is also valid for directed graphs. In the latter case, the degree of each vertex is simply counted as the sum of incoming plus outgoing arcs. Consequently, this proposition allows us to write the degree constraints (8) and (9) in P1 equivalently as follows:(15)∑i|i,j∈Azij+∑i|j,i∈Azji=δj,∀j∈V,(16)∑j∈Vδj=2∑j∈Vxj−1,(17)δj≤dxj,∀j∈V,(18)δj≥xj,∀j∈V,(19)δ∈0,∞n,where δj, for each j∈V, is a continuous nonnegative variable used to denote the degree of each vertex. Thus, we obtain another MILP model that we denote hereafter by P1h. Notice that the variable δj, for each j∈V, need not be defined as an integer variable. In fact, this is ensured by the left-hand side of constraint (15) and the right-hand side of constraint (16) which are both integer values. Constraints (17) and (18) restrict the value of variable δj between 1 and d if and only if node j∈V is part of the solution; otherwise, δj=0. Finally, constraint (19) is a domain constraint for each variable δj, j∈V.In order to formulate a flow-based model for the DCkMST problem, let H^=V∪r,A∪Ar be an expanded digraph obtained from H=V,A where we add to H^ an artificial root node r and a set of arcs Ar with zero costs from r to every node v∈V. The underlying idea is thus to construct an arborescence rooted at r while expanding all nodes in the solution with exactly one arc leaving r. For this purpose, we expand the vector of node and arc variables to x∈0,1n+1 and z∈0,1n+12, respectively. We also define continuous nonnegative flow variables f∈0,∞n+12, where fij denotes the amount of flow on arc i,j∈A∪Ar. Finally, if we represent the number of nodes in the solution by a nonnegative variable λ, then the idea is to send λ units of flow from r to these nodes, one unit of flow to every one of them. Consequently, a flow-based model can be verified by means of the following proposition.Proposition 5.
ModelP2 allows to obtain an optimal solution for the DCkMST problem with minimum cost:(20)P2:minx,z,f,λ∑i,j∈APijzijs.t.λ≥k,(21)∑j∈Vxj=λ,(22)∑i,j∈Azij=λ−1,(23)∑j|r,j∈Arfrj=λ,(24)∑j|r,j∈Arzrj=1,(25)∑i|i,j∈A∪Arfij−∑i|j,i∈A∪Arfji=xj,∀j∈V,(26)fij≤n−1zij,∀i,j∈A∪Ar,(27)∑i|i,j∈A∪Arzij+∑i|j,i∈A∪Arzji≤dxj,∀j∈V,(28)∑i|i,j∈A∪Arzij+∑i|j,i∈A∪Arzji≥xj,∀j∈V,(29)zij+zji≤xi,∀i,j∈A∪Ar,xr=1,(30)x∈0,1n+1,z∈0,1n+12,f∈0,∞n+12,λ≥0,where the objective function is the same as for P1.Proof.
In order to prove the correctness of modelP2, first we make the following observation.Observation 1.
The number of leaf nodes (nodes with degree equal to one) in any spanning tree graph withn≥3 vertices is at least 2 and at most n−1, corresponding to the cases of a line and a star graph, respectively.
Next, notice that constraints (20) and (21) guarantee that λ is at least k units of flow and that it is equal to the number of active nodes in the resulting subtree, respectively, whilst constraint (22) states that the total number of arcs in the solution should be equal to the number of active nodes minus one. Constraints (23) and (24) ensure that the total amount of flow going out from node r equals λ and that it must be moved through a unique arc r,j∈Ar,j∈V, respectively. Similarly, constraint (25) states that the incoming minus the outgoing flow equals one if node v∈V is part of the solution and equals zero otherwise. Constraint (26) ensures that the flow going out from i to j equals zero if and only if arc i,j∈A∪Ar is not part of the solution. Otherwise, it should be at most n−1 units. Constraints (27) and (28) are the degree constraints. Notice that the node index in each sum of these constraints applies for all arcs i,j∈A∪Ar in contrast to the node indexes in constraints (8) and (9) which only apply for the arcs i,j∈A. These degree constraints are valid since Observation 1 implies that the artificial node r can always be connected to a leaf node of the resulting subtree without affecting the maximum degree value d imposed to each node v∈V. Next, constraint (29) forces the fact that the root node r is always active. The remaining constraints are the same as for P1 including constraint (30) which is a domain constraint for the decision variables. Finally, notice that the simultaneous conditions imposed by constraints (22)–(26) ensure that the resulting digraph obtained is a tree [3, 13]. Again, an undirected subtree can be obtained by dropping the directions of the arcs.Proposition 6.
Constraints (27) and (9) are tighter than constraints (8) and (28), respectively.Proof.
Let us assume that nodet∈V is part of the resulting subtree. Then, we have that(31)d≥∑i|i,t∈A∪Arzit+∑i|t,i∈A∪Arzti=∑i|i,t∈Azit+zrt+∑i|t,i∈Azti+ztr≥∑i|i,t∈Azit+∑i|t,i∈Azti≥1.
Analogously as forP1, we can replace the degree constraints (27) and (28) in P2 by the set of constraints (15)–(19), leading to another MILP formulation we denote hereafter by P2h. Notice that, from the constructions of the above formulations, an exponential model can also be written as follows:(32)P3:minx,z∑i,j∈APijzijs.t.∑j∈Vxj≥k∑i,j∈Azij=∑j∈Vxj−1,(33)∑i|i,j∈Azij≤xj,∀j∈V,∑i,j∈ASzij≤S−1,∀S⊆V,(34)∑i|i,j∈Azij+∑i|j,i∈Azji≤dxj,∀j∈V,(35)∑i|i,j∈Azij+∑i|j,i∈Azji≥xj,∀j∈V,(36)zij+zji≤xi,∀i,j∈A,x∈0,1n,z∈0,1n2,where inequalities (33) represent subtour elimination constraints. Similarly, as for the above models P1 and P2, we replace the degree constraints (34) and (35) in P3 by the set of constraints (15)–(19). We denote by P3h this alternative exponential model. Hereafter, we also denote the corresponding linear programming (LP) relaxations of P1, P1h, P2, and P2h by LP1, LP1h, LP2, and LP2h, respectively.


## 3.1. A Feasible Solution for the DCkMST Problem

In Figure1(a), we depict an input complete graph instance with 10 nodes where the coordinates of each node are randomly generated within a square area of 1 km2, whereas in Figures 1(b) and 1(c), we show feasible solutions obtained for the instance in Figure 1(a) while using degree values of 3 and 2, respectively. For the sake of clarity, in Figures 1(b) and 1(c), we only draw the edges and highlight with green color the nodes which are part of the solutions, i.e., the nodes 1,2,3,5,6,7. Notice that the minimum cost value of an optimal solution in Figure 1(b) should be lower than or equal to the minimum cost value of an optimal solution in Figure 1(c) since the latter is also a feasible solution for the degree value of 3. On the opposite, the solution in Figure 1(b) is not feasible for the degree value of 2.Figure 1
An input complete graph instance with feasible solutions for the DCkMST problem while using degree values of 3 and 2, respectively: (a) complete graph instance composed of 10 nodes with random uniform costs; (b) feasible solution with node degree values of at most 3; (c) feasible solution with node degree values of at most 2.
(a)(b)(c)Proposition 1.
The following statements are true:(1)
Ifk=n and d≤n−2, the DCkMST problem reduces to the DCMST problem.(2)
Ifk<n and d≥k−1, the DCkMST problem reduces to the kMST problem.(3)
Ifk=n and d≥n−1, the problem reduces to the classical minimum spanning tree (MST) problem, which can be solved in polynomial time by a greedy algorithm [50, 51].
From Proposition1, it is easy to see that the DCkMST problem is NP-hard as it generalizes both the DCMST and KMST problems. Hereafter, we present the MILP models.


## 3.2. MILP Models

In order to formulate a first MILP model, letH=V,A be the directed graph obtained from GV,E, where each edge i,j∈E is replaced by the two arcs i,j and j,i∈A. Thus, we have the following proposition.Proposition 2.
ModelP1 allows to obtain an optimal solution for the DCkMST problem with minimum cost:(1)P1:minx,u,z∑i,j∈APijzij,(2)s.t.∑j∈Vxj≥k,(3)∑i,j∈Azij=∑j∈Vxj−1,(4)ui≤∑j∈Vxj,∀i∈V,(5)ui≥xi,∀i∈V,(6)∑i|i,j∈Azij≤xj,∀j∈V,(7)uj−ui−n−1zij−n−3zji≥2−n,∀i,j∈A,(8)∑i|i,j∈Azij+∑i|j,i∈Azji≤dxj,∀j∈V,(9)∑i|i,j∈Azij+∑i|j,i∈Azji≥xj,∀j∈V,(10)zij+zji≤xi,∀i,j∈A,(11)x∈0,1n,z∈0,1n2,u∈0,∞n,where P=Pij for each i,j∈A is defined as an input symmetric matrix denoting the connectivity costs between nodes i,j∈V. We also define the binary variable zij to be equal to one if and only if arc i,j∈A belongs to the solution of the resulting spanning tree and zij=0 otherwise.Proof.
Notice that the total connectivity cost is minimized in (1). Next, notice that constraints (2) and (3) ensure that at least k nodes from set V should belong to the output solution of the problem and that the number of arcs must be equal to the number of nodes minus one, respectively. For this purpose, we define the binary variable xj for each j∈V where xj=1 if and only if node j is in the solution and xj=0 otherwise. Subsequently, the constraints (4)–(7) ensure that the solution does not contain subtours. For this purpose, we also define the nonnegative variable uj for each j∈V. To see how these constraints avoid cycles, let us assume that the subset of nodes a1,a2,a3,…,at (t≤k) is part of the solution tree and that these nodes are connected according to the following sequence a1−a2−a3−,…,−at. Then, the set of arcs a1,a2,a2,a3,…,at−1,at is also part of the solution tree and za1,a2=za2,a3=za3,a4=⋯=zat−1,at=1. Also, by constraint (7), we have the following relation with the uj variables, ua1<ua2<ua3<⋯<uat. Notice that, for any node in this sequence, we cannot have any other incoming arc; otherwise, we would violate constraint (6) which ensures that each node cannot have more than one incoming arc if it belongs to the solution. In particular, if node a1 is the root node, then no incoming arc can be connected to it either since this would imply that uat<ua1, which is a contradiction. Finally, notice that any node that belongs to another branch of the resulting tree cannot be connected to any node of any other sequence of nodes forming another branch as this is prevented by constraint (6) too. Consequently, the resulting digraph has to be acyclic. Notice that the constraints (4) and (5) only apply for the nodes which are part of the solution by restricting the values of each variable uj,j∈V. Also, the degree constraints (8) and (9) ensure that the number of incoming plus outgoing arcs of each node is less than d∈2,…,k−2 where k≥4. Notice that the degree constraints are redundant if a particular node j∈V is not part of the solution, i.e., if xj=0. In fact, this also implies that ∑i|i,j∈Azij+∑i|j,i∈Azji=0. Next, constraint (10) indicates that at most one of the arcs i,j,j,i∈A can belong to the solution if and only if node i∈V belongs to the solution. Constraint (11) is a domain constraint for the decision variables. Finally, it is observed that an undirected subtree graph can be obtained by dropping the direction of each arc in the resulting digraph, as required.Proposition 3.
Constraint (7) is tighter than constraint (35) presented in [2].Proof.
Constraint (35) in [2] is written as(12)uj−ui≥1−n1−zij,∀i,j∈A.
However, these constraints can be equivalently written as(13)uj−ui−n−1zij≥2−n,∀i,j∈A.
Then, we have that(14)uj−ui−n−1zij≥uj−ui−n−1zij−n−3zji≥2−n,∀i,j∈A.
More details related to constraint (7), when used for the traveling salesman problem and extended to other types of vehicle routing problems, can be consulted in [52]. Notice that model P1 is a directed graph formulation constructed based on a Miller–Tucker–Zemlin characterization [2, 3, 13].Proposition 4 (see [23]).
In a finite simple undirected graph, the sum of the degrees of every vertexv∈V is twice the number of edges.Proposition4 is known as the handshaking lemma, and it is also valid for directed graphs. In the latter case, the degree of each vertex is simply counted as the sum of incoming plus outgoing arcs. Consequently, this proposition allows us to write the degree constraints (8) and (9) in P1 equivalently as follows:(15)∑i|i,j∈Azij+∑i|j,i∈Azji=δj,∀j∈V,(16)∑j∈Vδj=2∑j∈Vxj−1,(17)δj≤dxj,∀j∈V,(18)δj≥xj,∀j∈V,(19)δ∈0,∞n,where δj, for each j∈V, is a continuous nonnegative variable used to denote the degree of each vertex. Thus, we obtain another MILP model that we denote hereafter by P1h. Notice that the variable δj, for each j∈V, need not be defined as an integer variable. In fact, this is ensured by the left-hand side of constraint (15) and the right-hand side of constraint (16) which are both integer values. Constraints (17) and (18) restrict the value of variable δj between 1 and d if and only if node j∈V is part of the solution; otherwise, δj=0. Finally, constraint (19) is a domain constraint for each variable δj, j∈V.In order to formulate a flow-based model for the DCkMST problem, let H^=V∪r,A∪Ar be an expanded digraph obtained from H=V,A where we add to H^ an artificial root node r and a set of arcs Ar with zero costs from r to every node v∈V. The underlying idea is thus to construct an arborescence rooted at r while expanding all nodes in the solution with exactly one arc leaving r. For this purpose, we expand the vector of node and arc variables to x∈0,1n+1 and z∈0,1n+12, respectively. We also define continuous nonnegative flow variables f∈0,∞n+12, where fij denotes the amount of flow on arc i,j∈A∪Ar. Finally, if we represent the number of nodes in the solution by a nonnegative variable λ, then the idea is to send λ units of flow from r to these nodes, one unit of flow to every one of them. Consequently, a flow-based model can be verified by means of the following proposition.Proposition 5.
ModelP2 allows to obtain an optimal solution for the DCkMST problem with minimum cost:(20)P2:minx,z,f,λ∑i,j∈APijzijs.t.λ≥k,(21)∑j∈Vxj=λ,(22)∑i,j∈Azij=λ−1,(23)∑j|r,j∈Arfrj=λ,(24)∑j|r,j∈Arzrj=1,(25)∑i|i,j∈A∪Arfij−∑i|j,i∈A∪Arfji=xj,∀j∈V,(26)fij≤n−1zij,∀i,j∈A∪Ar,(27)∑i|i,j∈A∪Arzij+∑i|j,i∈A∪Arzji≤dxj,∀j∈V,(28)∑i|i,j∈A∪Arzij+∑i|j,i∈A∪Arzji≥xj,∀j∈V,(29)zij+zji≤xi,∀i,j∈A∪Ar,xr=1,(30)x∈0,1n+1,z∈0,1n+12,f∈0,∞n+12,λ≥0,where the objective function is the same as for P1.Proof.
In order to prove the correctness of modelP2, first we make the following observation.Observation 1.
The number of leaf nodes (nodes with degree equal to one) in any spanning tree graph withn≥3 vertices is at least 2 and at most n−1, corresponding to the cases of a line and a star graph, respectively.
Next, notice that constraints (20) and (21) guarantee that λ is at least k units of flow and that it is equal to the number of active nodes in the resulting subtree, respectively, whilst constraint (22) states that the total number of arcs in the solution should be equal to the number of active nodes minus one. Constraints (23) and (24) ensure that the total amount of flow going out from node r equals λ and that it must be moved through a unique arc r,j∈Ar,j∈V, respectively. Similarly, constraint (25) states that the incoming minus the outgoing flow equals one if node v∈V is part of the solution and equals zero otherwise. Constraint (26) ensures that the flow going out from i to j equals zero if and only if arc i,j∈A∪Ar is not part of the solution. Otherwise, it should be at most n−1 units. Constraints (27) and (28) are the degree constraints. Notice that the node index in each sum of these constraints applies for all arcs i,j∈A∪Ar in contrast to the node indexes in constraints (8) and (9) which only apply for the arcs i,j∈A. These degree constraints are valid since Observation 1 implies that the artificial node r can always be connected to a leaf node of the resulting subtree without affecting the maximum degree value d imposed to each node v∈V. Next, constraint (29) forces the fact that the root node r is always active. The remaining constraints are the same as for P1 including constraint (30) which is a domain constraint for the decision variables. Finally, notice that the simultaneous conditions imposed by constraints (22)–(26) ensure that the resulting digraph obtained is a tree [3, 13]. Again, an undirected subtree can be obtained by dropping the directions of the arcs.Proposition 6.
Constraints (27) and (9) are tighter than constraints (8) and (28), respectively.Proof.
Let us assume that nodet∈V is part of the resulting subtree. Then, we have that(31)d≥∑i|i,t∈A∪Arzit+∑i|t,i∈A∪Arzti=∑i|i,t∈Azit+zrt+∑i|t,i∈Azti+ztr≥∑i|i,t∈Azit+∑i|t,i∈Azti≥1.
Analogously as forP1, we can replace the degree constraints (27) and (28) in P2 by the set of constraints (15)–(19), leading to another MILP formulation we denote hereafter by P2h. Notice that, from the constructions of the above formulations, an exponential model can also be written as follows:(32)P3:minx,z∑i,j∈APijzijs.t.∑j∈Vxj≥k∑i,j∈Azij=∑j∈Vxj−1,(33)∑i|i,j∈Azij≤xj,∀j∈V,∑i,j∈ASzij≤S−1,∀S⊆V,(34)∑i|i,j∈Azij+∑i|j,i∈Azji≤dxj,∀j∈V,(35)∑i|i,j∈Azij+∑i|j,i∈Azji≥xj,∀j∈V,(36)zij+zji≤xi,∀i,j∈A,x∈0,1n,z∈0,1n2,where inequalities (33) represent subtour elimination constraints. Similarly, as for the above models P1 and P2, we replace the degree constraints (34) and (35) in P3 by the set of constraints (15)–(19). We denote by P3h this alternative exponential model. Hereafter, we also denote the corresponding linear programming (LP) relaxations of P1, P1h, P2, and P2h by LP1, LP1h, LP2, and LP2h, respectively.


## 4. Proposed Algorithms

In this section, first we present and explain the pseudocode for the exact iterative approach used to solve the exponential modelsP3 and P3h. Then, for the sake of clarity, we present the modified penalty approach proposed in [30] which allows to obtain feasible solutions for the DCMST problem. Then, we present our VNS and ACO algorithms while using traditional and embedded Q-learning random local search strategies. Finally, we present a pure Q-learning-based algorithm.


### 4.1. Exact Iterative Algorithm

As mentioned in Section1, we solve our models P3 and P3h with an exact iterative method adapted from [2, 24]. The method is simple and can be described as follows. It consists of solving first the MILP model P3 (or P3h) without subtour elimination constraints (SECs). Then, new cycles are obtained from the current optimal solution found with a depth-first search procedure [53]. Next, for each cycle found, we write a new constraint of the form (33) and add it to the feasible set of P3 (or P3h). Finally, the model P3 (or P3h) including all the new added inequalities is reoptimized. This process goes on until the current solution found does not contain any cycle. When this condition is met, it means the optimal solution has been reached. This procedure is depicted in Algorithm 1.Algorithm 1: Iterative procedure for solving P3 (or P3h)
Data: a problem instance of P3 (or P3h)Result: the optimal solution of P3 (or P3h)SolveP3 (resp., P3h) without using any constraint of the form (33)Find cycles on the resulting digraph obtained by using a depth-first search algorithmWhile (no cycles remain in the optimal solution of P3 (resp., P3h)) doFor each cycle found, write a new constraint of the form (33) and add it to the feasible set of P3 (resp., P3h)SolveP3 (resp., P3h)Return: the optimal solution found and the objective function valueThe convergence proof of this method is reported in Theorem 2 in [24]. Notice that, according to Theorem 2 in [24], the solutions obtained within each iteration of Algorithm 1 are lower bounds for the optimal solution of the problem. Further notice that a depth-first search algorithm will always find cycles within each iteration of Algorithm 1 if there exists at least one in the resulting digraph. This fact ensures the convergence of the method in a finite number of iterations (see Property 1 in [24]). We refer the reader to the works in [2, 24] for further details on how we obtain cycles.


### 4.2. Modified Penalty Approach

In particular, within each iteration of our VNS algorithms, for each randomk−tuple of vertices generated, we obtain degree-constrained spanning trees by using a modified penalty approach proposed in [30]. The procedure of this method is depicted in Algorithm 2, and it is explained as follows.Algorithm 2: Modified penalty approach to obtain feasible degree-constrained spanning trees.
Data: a problem instance of the DCMST problem.Result: a feasible solution for the DCMST problem. Find a minimum spanning tree using the Kruskal method [53]SetSW=trueWhile (SW) doSW=falseForEach (j∈V) doCompute the degree of nodej and save this value in variable djIf (dj>d) thenSW=trueUpdate all incident arcsi,j∈A for node j as follows:Pij=Pij+θPij−eMin/eMax−eMineMaxPij=PjiIf (SW) thenFind a minimum spanning tree using the Kruskal method [53]Return feasible solution obtained and its objective function valueInitially, we find a minimum spanning tree using the Kruskal method [53]. Then, we enter into a while loop doing the following. First, we check if each node satisfies the degree condition. If it is not the case for a particular node, then we update all the weights of the arcs which are incident to it using the formula Pij=Pij+θPij−eMin/eMax−eMineMax, where θ represents an arbitrary parameter which can be drawn from the interval 0;1 and eMin and eMax are the smallest and largest weights of the current spanning tree obtained, respectively. Subsequently, we symmetrize the updated arcs of matrix Pij in order to avoid choosing previous arcs with opposite directions. Finally, if there still exists a node violating the degree condition, we continue with the process and find a new minimum spanning tree with the Kruskal method. This process continues until a feasible spanning tree is obtained where all nodes satisfy the degree condition. As it can be observed, Algorithm 2 mainly consists of penalizing iteratively the edges (arcs) which are incident to the vertices with degree violations until a feasible spanning tree solution for the problem is obtained.


### 4.3. VNS Algorithms

We now present our VNS algorithms with the classical near-far random local search strategy and then using an embedded Q-learning strategy.


#### 4.3.1. Classical Random Local Search Strategy

Our VNS procedure is depicted in Algorithm3. As it can be observed, the method requires an instance of the DCkMST problem and a particular degree with value d≥3. As in this paper, we only consider degree values of d∈2,3, then we only solve instances for a degree value of d=3 with this VNS algorithm. The method is simple and can be described as follows. First, it checks if k equals n; if this is the case, then we execute Algorithm 2 using parameters n,P, and d. For this particular case, it means the algorithm finds a feasible solution for the DCMST problem including all nodes of set V. On the opposite, if k<n, then we randomly choose an initial k-tuple of vertices from V in order to form the set As. Let Ac be the complement of As.Algorithm 3: VNS algorithm for the DCkMST problem using classical random local search strategy.
Data: an instance of the DCkMST problem using degree d=3.Result: a feasible solution and its objective function value.If (k=n) thenExecute Algorithm2 using parameters n,P,dSave feasible solution found and its objective function valueElseClassical random local search strategy:Generate an initial randomk-tuple of vertices. Let As denote this set of vertices and Ac its complement, P′=PAs,AsExecute Algorithm2 using parameters k,P′,dSave the initial feasible solution found and its objective function valueAsOP=As, AcOP=AcindK=1, cont=0, CpuTime=0While (CpuTime≤maxTime) doiter=iter+1Fori=1 to indKdoInterchange randomly an element ofAs with an element of AcP′=PAs,As, execute Algorithm 2 using parameters k,P′,dIf (a better solution is obtained) thenSave the new solution and setAsOP=As, AcOP=AciterOp=iter, cont=0, indK=1, Optime=Optime+CpuTime, CpuTime=0Elsecont=cont+1, As=AsOP, Ac=AcOPIf (cont≥1) thencont=0If (indK<indKMax) thenindK=indK+1ElseindK=1Return best feasible solution obtained and its objective function valueNext, we construct the submatrixP′ with the rows and columns of P corresponding to the nodes in As and execute Algorithm 2 using parameters k, P′, and d. This allows us to save an initial feasible solution together with its objective function value. Notice that the cardinality of set As is k, whilst the cardinality of set Ac is n−k. The current sets As and Ac are also saved in AsOP and AcOP, respectively. Subsequently, we perform the following steps iteratively, while the CpuTime variable is less than or equal to the maximum CPU time allowed which is controlled with parameter maxTime. Inside this loop, we randomly interchange an element of As with an element of Ac a predefined number of times which is controlled with parameter indK that is initially set to a value of one. Then, we check if a better solution is obtained. If so, we save the new solution found, its objective function value, update sets AsOP and AcOP, and reset the CpuTime variable to zero. Otherwise, we keep the previous best solution found and go back to the interchange step. Notice that parameter indK controls the number of swap moves the algorithm performs between sets As and Ac. It turns out that this step corresponds to a classical near-far random local search strategy related to VNS algorithm [28, 29]. The larger the number of swap moves, the wider the feasible space being explored. On the opposite, the smaller the number of moves, the closer the solutions which are being tested. In particular, in Algorithm 3, this is handled with indK variable which is increased by one unit each time a worse solution is obtained. This variable can be increased with up to a predefined value of indKMax−1. And we reset this parameter to a value of one each time a better solution is obtained in order to restart the search from a local to a wider feasible region. Finally, if indK equals indKMax, then we also reset indK to a value of one.


#### 4.3.2. Embedded Q-Learning Strategy

Our Q-learning strategy is simple and consists of embedding a reinforcement learning approach in Algorithm3 instead of using classical random local search. Notice that the Q-learning approach was initially proposed in [31]. This approach mainly consists of an agent, a set of states S, and a set A of actions per state whereby performing an action a∈A; the agent moves from one state to another one according to a reward value, and the goal of the agent is to maximize its total future reward. The potential reward that an agent can obtain is computed as a weighted sum of the expected values of the rewards of all future steps starting from any of the current states. In order to embed this Q-learning approach in Algorithm 3, before entering the inner while loop of Algorithm 3, we need to initialize the matrix variable Q=QindKM,indKM with zero entries and also to declare the two scalar variables indKA and indicator.In particular, the value ofindicator variable, which can be true or false, is used to perform or not a piece of the code, whereas variables indK and indKA represent the number of swap moves from one to indKM to be performed between sets As and Ac. Thus, the underlying idea is to move from a particular number of swap moves to another one by using the aforementioned reinforcement learning strategy. Under this setting, the values of indKA and indK represent two consecutive states where an agent can be in different moments. For example, consider the case when indKA=2 and indK=8, then we should read it as an agent that is moving from state 2 to state 8. In our problem, it means we are first interchanging 2 elements and then 8 elements between sets As and Ac. Consequently, we can replace the inner code of the Else statement corresponding to the classical random local search strategy step of Algorithm 3 with the code of Algorithm 4. By doing so, we provide learning capabilities to Algorithm 3 in order to perform the random local search. In Algorithm 4, the parameters λ, μ, and Prob∈0;1 are the learning rate, discount factor, and probability of moving randomly from one state to another one [31]. Finally, the function expressions Rand, RandintindKM, and argmaxj∈1,…,indKMQindKA,j return a random fractional value between zero and one, a random integer value between one and indKM, and the jth column position with current maximum value in matrix Q while moving from state indKA, respectively.Algorithm 4: Embedded Q-learning code.
Generate an initial randomk-tuple of vertices. Let As denote this set of vertices and Ac its complement P′=PAs,AsExecute Algorithm2 using parameters k,P′,dSave the initial feasible solution found and its objective function valueAsOP=As, AcOP=AcindK=1, indKA=1, Q=zerosindKM,indKM, indicator=false, CpuTime=0While (CpuTime≤maxTime) doiter=iter+1Fori=1 to indKdoInterchange randomly an element ofAs with an element of AcSetP′=PAs,As and run Algorithm 2 using parameters k,P′,dIf (a better solution is obtained) thenSave the new solution and setAsOP=As, AcOP=Ac, indicator=trueiterOp=iter, Optime=Optime+CpuTime, CpuTime=0ElseAs=AsOP, Ac=AcOPIf (indicator) thenQindKA,indK=1−λQindKA,indK+λReward+μmaxj∈1,…,indKMQindK,jindKA=indKIf (Rand<Prob) thenindK=RandintindKMElseindK=argmaxj∈1,…,indKMQindKA,jindicator=false


### 4.4. ACO Algorithms

We now present the ACO algorithms. More precisely, first we present our classical ACO based strategy and then we present an adapted version of the generic ant-Q algorithm which is based on Q-learning strategy [26, 27]. Finally, we present a pure Q-learning-based algorithm. All these methods allow finding feasible solutions for the DCkMST problem while using a degree value of d=2.


#### 4.4.1. Classical Ant Colony Optimization Strategy

Our first procedure is depicted in Algorithm5.Algorithm 5: Classical ACO algorithm for the DC kMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersρ, α, β, q and the set of ants (Ants)CpuTime=0, iter=0, LBest=∞, BestCPU=0τij=rand, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Δτija=0,∀i,j∈V, a∈AntsForEacha∈AntsdoConstruction of each ant tourRandomly assign an initial positionj∈V for the ant a according to graph G and construct a tour Ta using all remaining nodes of G.Find best solutionCompute the lengthLa of TaDenote byLakt the length of k consecutive nodes in Ta starting from node tComputeLak=minLak1,…,Lakn−k+1If (LBest>Lak) thenSetLBest=Lak and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Update accumulated pheromone matrixForEachi,j∈TadoΔτija=q/La, Δτjia=q/LaUpdate pheromone matrixForEachi,j∈Adoτi,j=1−ρτi,j+∑a∈AntsΔτijaReturn best feasible solution found and its objective function valueThe method requires an instance of the DCkMST problem. In step one, it first initializes all required parameters and variables. In particular, ρ, α, β, and q denote the pheromone evaporation coefficient, the influence of the pheromone, influence of visibility matrix, and a constant real value, respectively. Each entry of the pheromone matrix τ=τij for all i,j∈V represents the amount of pheromone deposited by the ants in the edge i,j for a state transition going from node i to j. Notice that each of these entries is initialized to a random value between zero and one, which is performed by using the Rand function. Similarly, the visibility matrix η=ηij for all i,j∈V represents the desirability value of edge i,j while choosing the state transition from i to j. In particular, this value is computed by ηij=1/Pij. The set of ants is denoted by Ants. The variables CpuTime, iter, LBest, and BestCPU denote the current CPU time, the number of iterations, and best current objective function value found so far and its exact CPU time in which this solution is obtained. Recall that, within ACO algorithms, each ant constructs a tour solution starting from a particular initial node while covering all nodes of the graph. For this purpose, each ant selects the next edge in its tour according to the following probability:(37)pija=τijαηijβ∑z∈Raτizαηizβ,where Ra denotes the unvisited set of nodes of ant a. As it can be observed, this probability depends on both the pheromone and visibility values [26].Subsequently, in step two of Algorithm5, we enter into a while loop in which we perform the following steps. First, we initialize to a value of zero each entry of matrix Δτ=Δτija, i,j∈V, a∈Ants, where each entry in this matrix represents the amount of pheromone deposited by ant a∈Ants. Next, for each ant, we construct a tour Ta starting from a randomly assigned node and compute its length La. Then, for each arc in Ta, we deposit the amount of pheromone of ant a. This value is computed by Δτija=q/La. Then, we compute the shortest tour of length k according to Ta and save it as the best tour found so far if its objective function value is less than LBest. For this purpose, we denote by Lakt the length of k consecutive nodes in Ta starting from node t. Notice that t=1,…,n−k+1. Also notice that the variable CpuTime is reset to a value of zero each time a better solution is obtained. This allows Algorithm 5 to run for another maxTime unit of time with the hope of finding better solutions. Next, we update each arc of the pheromone matrix τ according to parameter ρ and matrix Δτ. Finally, when the while loop is finished, we return the best feasible solution obtained together with its objective function value.


#### 4.4.2. Ant-Q-Based Approach

Now, we present our adapted version of the generic Ant-Q algorithm which is based on Q-learning strategy [25, 27]. This procedure is depicted in Algorithm 6.Algorithm 6: Ant-Q algorithm for the DCkMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersα, β, q, Prob, λ, μ, and the set of ants (Ants)CpuTime=0, iter=0, LBest=∞, BestCPU=0AQi,j=0, Δi,j=0, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Construct ant toursForEacha∈AntsdoRandomly choose a nodej∈VTa=∅, Ta=Ta∪jFori=1 to n−1doForEacha∈AntsdoIfRand<ProbthenRandomly choose an unvisited node of the ant and add it toTa.Elsej=argmaxrallowedAQTai,rα1/PTai,rβTa=Ta∪jUpdate AQ-valuesAQTai,Tai+1=1−λAQTai,Tai+1+λμmaxjallowedAQTai+1,jAQTai+1,Tai=AQTai,Tai+1Find best solutionForEacha∈AntsdoDenote byLakt the length of k consecutive nodes in Ta starting from node tComputeLak=minLak1,…,Lakn−k+1If (LBest>Lak) thenSetLBest=Lak and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Reinforcement of ant toursForEachi,j∈TadoΔi,j=q/Lak, Δj,i=q/LakReinforcement of best solutionForEachi,j∈TkdoΔi,j=q/LBest, Δj,i=q/LBestDelayed reinforcement of AQ valuesForEacha∈AntsdoV=1,…,nFori=1 to n−1doV=V−Tai, V=V−Tai+1IfV≠∅thenAQTai,Tai+1=1−λAQTai,Tai+1+λΔTai,Tai+1+μmaxAQTai+1,VAQTai+1,Tai=AQTai,Tai+1Return best feasible solution found and its objective function valueSimilarly to Algorithm5, in step one of the Ant-Q Algorithm 6, we first initialize the required parameters and variables. In this case, the parameters α and β allow us to weigh the relative importance of learned AQ and visibility values, respectively. Parameter q is a nonnegative constant value, whilst parameters Prob, λ, and μ represent a probability value, a learning step, and a discount factor, respectively. Variables CpuTime, iter, LBest, BestCPU are analogously defined as for Algorithm 5. In particular, matrices AQ=AQi,j and Δ=Δi,j, for all i,j∈V, denote the Q-learning matrix and the delayed reinforcement matrices associated with the Ant-Q learning method [31].In step two of Algorithm6, we enter into a while loop and perform the following substeps. First, we construct a tour for each ant a∈Ants. This is performed iteratively while choosing a new unvisited node j randomly or according to the following expression:(38)j=argmaxrallowedAQTai,rα1PTai,rβ,where Tai denotes the current ith position of ant a in its tour Ta. Then, we update the learning matrix AQ according to all the edges of each ant tour and save the best solution obtained with cardinality k. Next, we proceed with reinforcement steps that are applied to all the edges belonging to each ant tour and best solution found so far [25, 27]. Finally, the algorithm returns the best feasible solution obtained and its objective function value.


### 4.5. A Pure Q-Learning-Based Algorithm

Now, we present a pure Q-Learning based approach that allows obtaining feasible solutions for the DCkMST problem while using d=2. This method is depicted in Algorithm 7, and it is described as follows.Algorithm 7: Pure Q-learning-based approach for the DCkMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersλ, μ, ProbCpuTime=0, iter=0, LBest=∞, BestCPU=0, Qi,j=0, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Construct a unique tourRandomly choose a nodej∈VT=∅, T=T∪j, V=V−jFori=1 to n−1doIfRand<ProbthenRandomly choose a nodeu∈VElseu=argmaxr∈VQj,rT=T∪u, V=V−uIfV≠∅thenQj,u=1−λQj,u+λ1/Pj,u+μmaxr∈VQu,rj=uFind best solutionDenote byLkt the length of k consecutive nodes in T starting from node tComputeLk=minLk1,…,Lkn−k+1If (LBest>Lk) thenSetLBest=Lk and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Return best feasible solution found and its objective function valueSimilarly to Algorithm6, the first step of Algorithm 7 consists of initializing parameters λ, μ, and Prob, which represent the learning rate, discount factor, and a given probability value, respectively. We also define and initialize the required variables CpuTime, iter, LBest, BestCPU, and Q=Qi,j for all i,j∈V. These variables allow us to handle the current CPU time, the number of iterations, the best objective function value obtained so far, CPU time of the best current solution found, and the Q-learning matrix, respectively. Next, in step two of Algorithm 7, we iteratively construct a unique tour T with cardinality n and evaluate each subtour composed of k consecutive nodes in T. These steps are performed, while the current CPU time value is lower than the maximum value allowed which is denoted by maxTime. In order to obtain the best subtour of cardinality k from each tour T generated, we let Lkt denote the length of k consecutive nodes in T, starting from node t, and let Lk=minLk1,…,Lkn−k+1 be the minimum length value obtained. Notice that index t=1,…,n−k+1; otherwise, no subtour of cardinality k can be obtained from T. If a better solution is obtained, we save its length value in LBest, the subtour solution in Tk, the number of iterations in which this new solution is obtained, and reset the current CPU time value to a value of zero. The latter allows the algorithm to search for another maxTime unit of time. Also, notice that the construction of tour T requires to add at each step a new unvisited node u randomly or with maximum value Qj,u where j denotes the last node added to T. For this purpose, we remove node u from V at each step. The unvisited node u is chosen randomly if the Rand function generates a value in the interval 0;1 which is lower than Prob; otherwise, it is chosen according to the maximum value in the Q-learning matrix. Next, if set V is not empty, we update the corresponding entry in the Q-learning matrix as follows:(39)Qj,u=1−λQj,u+λ1Pj,u+μmaxr∈VQu,r.Notice that equation (39) is analog to the classical Q-learning algorithm [31]. In particular, the term 1/Pj,u represents the reward value obtained when going from node j to node u∈V. Finally, the algorithm returns the best feasible solution obtained and its objective function value.


## 4.1. Exact Iterative Algorithm

As mentioned in Section1, we solve our models P3 and P3h with an exact iterative method adapted from [2, 24]. The method is simple and can be described as follows. It consists of solving first the MILP model P3 (or P3h) without subtour elimination constraints (SECs). Then, new cycles are obtained from the current optimal solution found with a depth-first search procedure [53]. Next, for each cycle found, we write a new constraint of the form (33) and add it to the feasible set of P3 (or P3h). Finally, the model P3 (or P3h) including all the new added inequalities is reoptimized. This process goes on until the current solution found does not contain any cycle. When this condition is met, it means the optimal solution has been reached. This procedure is depicted in Algorithm 1.Algorithm 1: Iterative procedure for solving P3 (or P3h)
Data: a problem instance of P3 (or P3h)Result: the optimal solution of P3 (or P3h)SolveP3 (resp., P3h) without using any constraint of the form (33)Find cycles on the resulting digraph obtained by using a depth-first search algorithmWhile (no cycles remain in the optimal solution of P3 (resp., P3h)) doFor each cycle found, write a new constraint of the form (33) and add it to the feasible set of P3 (resp., P3h)SolveP3 (resp., P3h)Return: the optimal solution found and the objective function valueThe convergence proof of this method is reported in Theorem 2 in [24]. Notice that, according to Theorem 2 in [24], the solutions obtained within each iteration of Algorithm 1 are lower bounds for the optimal solution of the problem. Further notice that a depth-first search algorithm will always find cycles within each iteration of Algorithm 1 if there exists at least one in the resulting digraph. This fact ensures the convergence of the method in a finite number of iterations (see Property 1 in [24]). We refer the reader to the works in [2, 24] for further details on how we obtain cycles.


## 4.2. Modified Penalty Approach

In particular, within each iteration of our VNS algorithms, for each randomk−tuple of vertices generated, we obtain degree-constrained spanning trees by using a modified penalty approach proposed in [30]. The procedure of this method is depicted in Algorithm 2, and it is explained as follows.Algorithm 2: Modified penalty approach to obtain feasible degree-constrained spanning trees.
Data: a problem instance of the DCMST problem.Result: a feasible solution for the DCMST problem. Find a minimum spanning tree using the Kruskal method [53]SetSW=trueWhile (SW) doSW=falseForEach (j∈V) doCompute the degree of nodej and save this value in variable djIf (dj>d) thenSW=trueUpdate all incident arcsi,j∈A for node j as follows:Pij=Pij+θPij−eMin/eMax−eMineMaxPij=PjiIf (SW) thenFind a minimum spanning tree using the Kruskal method [53]Return feasible solution obtained and its objective function valueInitially, we find a minimum spanning tree using the Kruskal method [53]. Then, we enter into a while loop doing the following. First, we check if each node satisfies the degree condition. If it is not the case for a particular node, then we update all the weights of the arcs which are incident to it using the formula Pij=Pij+θPij−eMin/eMax−eMineMax, where θ represents an arbitrary parameter which can be drawn from the interval 0;1 and eMin and eMax are the smallest and largest weights of the current spanning tree obtained, respectively. Subsequently, we symmetrize the updated arcs of matrix Pij in order to avoid choosing previous arcs with opposite directions. Finally, if there still exists a node violating the degree condition, we continue with the process and find a new minimum spanning tree with the Kruskal method. This process continues until a feasible spanning tree is obtained where all nodes satisfy the degree condition. As it can be observed, Algorithm 2 mainly consists of penalizing iteratively the edges (arcs) which are incident to the vertices with degree violations until a feasible spanning tree solution for the problem is obtained.


## 4.3. VNS Algorithms

We now present our VNS algorithms with the classical near-far random local search strategy and then using an embedded Q-learning strategy.


### 4.3.1. Classical Random Local Search Strategy

Our VNS procedure is depicted in Algorithm3. As it can be observed, the method requires an instance of the DCkMST problem and a particular degree with value d≥3. As in this paper, we only consider degree values of d∈2,3, then we only solve instances for a degree value of d=3 with this VNS algorithm. The method is simple and can be described as follows. First, it checks if k equals n; if this is the case, then we execute Algorithm 2 using parameters n,P, and d. For this particular case, it means the algorithm finds a feasible solution for the DCMST problem including all nodes of set V. On the opposite, if k<n, then we randomly choose an initial k-tuple of vertices from V in order to form the set As. Let Ac be the complement of As.Algorithm 3: VNS algorithm for the DCkMST problem using classical random local search strategy.
Data: an instance of the DCkMST problem using degree d=3.Result: a feasible solution and its objective function value.If (k=n) thenExecute Algorithm2 using parameters n,P,dSave feasible solution found and its objective function valueElseClassical random local search strategy:Generate an initial randomk-tuple of vertices. Let As denote this set of vertices and Ac its complement, P′=PAs,AsExecute Algorithm2 using parameters k,P′,dSave the initial feasible solution found and its objective function valueAsOP=As, AcOP=AcindK=1, cont=0, CpuTime=0While (CpuTime≤maxTime) doiter=iter+1Fori=1 to indKdoInterchange randomly an element ofAs with an element of AcP′=PAs,As, execute Algorithm 2 using parameters k,P′,dIf (a better solution is obtained) thenSave the new solution and setAsOP=As, AcOP=AciterOp=iter, cont=0, indK=1, Optime=Optime+CpuTime, CpuTime=0Elsecont=cont+1, As=AsOP, Ac=AcOPIf (cont≥1) thencont=0If (indK<indKMax) thenindK=indK+1ElseindK=1Return best feasible solution obtained and its objective function valueNext, we construct the submatrixP′ with the rows and columns of P corresponding to the nodes in As and execute Algorithm 2 using parameters k, P′, and d. This allows us to save an initial feasible solution together with its objective function value. Notice that the cardinality of set As is k, whilst the cardinality of set Ac is n−k. The current sets As and Ac are also saved in AsOP and AcOP, respectively. Subsequently, we perform the following steps iteratively, while the CpuTime variable is less than or equal to the maximum CPU time allowed which is controlled with parameter maxTime. Inside this loop, we randomly interchange an element of As with an element of Ac a predefined number of times which is controlled with parameter indK that is initially set to a value of one. Then, we check if a better solution is obtained. If so, we save the new solution found, its objective function value, update sets AsOP and AcOP, and reset the CpuTime variable to zero. Otherwise, we keep the previous best solution found and go back to the interchange step. Notice that parameter indK controls the number of swap moves the algorithm performs between sets As and Ac. It turns out that this step corresponds to a classical near-far random local search strategy related to VNS algorithm [28, 29]. The larger the number of swap moves, the wider the feasible space being explored. On the opposite, the smaller the number of moves, the closer the solutions which are being tested. In particular, in Algorithm 3, this is handled with indK variable which is increased by one unit each time a worse solution is obtained. This variable can be increased with up to a predefined value of indKMax−1. And we reset this parameter to a value of one each time a better solution is obtained in order to restart the search from a local to a wider feasible region. Finally, if indK equals indKMax, then we also reset indK to a value of one.


### 4.3.2. Embedded Q-Learning Strategy

Our Q-learning strategy is simple and consists of embedding a reinforcement learning approach in Algorithm3 instead of using classical random local search. Notice that the Q-learning approach was initially proposed in [31]. This approach mainly consists of an agent, a set of states S, and a set A of actions per state whereby performing an action a∈A; the agent moves from one state to another one according to a reward value, and the goal of the agent is to maximize its total future reward. The potential reward that an agent can obtain is computed as a weighted sum of the expected values of the rewards of all future steps starting from any of the current states. In order to embed this Q-learning approach in Algorithm 3, before entering the inner while loop of Algorithm 3, we need to initialize the matrix variable Q=QindKM,indKM with zero entries and also to declare the two scalar variables indKA and indicator.In particular, the value ofindicator variable, which can be true or false, is used to perform or not a piece of the code, whereas variables indK and indKA represent the number of swap moves from one to indKM to be performed between sets As and Ac. Thus, the underlying idea is to move from a particular number of swap moves to another one by using the aforementioned reinforcement learning strategy. Under this setting, the values of indKA and indK represent two consecutive states where an agent can be in different moments. For example, consider the case when indKA=2 and indK=8, then we should read it as an agent that is moving from state 2 to state 8. In our problem, it means we are first interchanging 2 elements and then 8 elements between sets As and Ac. Consequently, we can replace the inner code of the Else statement corresponding to the classical random local search strategy step of Algorithm 3 with the code of Algorithm 4. By doing so, we provide learning capabilities to Algorithm 3 in order to perform the random local search. In Algorithm 4, the parameters λ, μ, and Prob∈0;1 are the learning rate, discount factor, and probability of moving randomly from one state to another one [31]. Finally, the function expressions Rand, RandintindKM, and argmaxj∈1,…,indKMQindKA,j return a random fractional value between zero and one, a random integer value between one and indKM, and the jth column position with current maximum value in matrix Q while moving from state indKA, respectively.Algorithm 4: Embedded Q-learning code.
Generate an initial randomk-tuple of vertices. Let As denote this set of vertices and Ac its complement P′=PAs,AsExecute Algorithm2 using parameters k,P′,dSave the initial feasible solution found and its objective function valueAsOP=As, AcOP=AcindK=1, indKA=1, Q=zerosindKM,indKM, indicator=false, CpuTime=0While (CpuTime≤maxTime) doiter=iter+1Fori=1 to indKdoInterchange randomly an element ofAs with an element of AcSetP′=PAs,As and run Algorithm 2 using parameters k,P′,dIf (a better solution is obtained) thenSave the new solution and setAsOP=As, AcOP=Ac, indicator=trueiterOp=iter, Optime=Optime+CpuTime, CpuTime=0ElseAs=AsOP, Ac=AcOPIf (indicator) thenQindKA,indK=1−λQindKA,indK+λReward+μmaxj∈1,…,indKMQindK,jindKA=indKIf (Rand<Prob) thenindK=RandintindKMElseindK=argmaxj∈1,…,indKMQindKA,jindicator=false


## 4.3.1. Classical Random Local Search Strategy

Our VNS procedure is depicted in Algorithm3. As it can be observed, the method requires an instance of the DCkMST problem and a particular degree with value d≥3. As in this paper, we only consider degree values of d∈2,3, then we only solve instances for a degree value of d=3 with this VNS algorithm. The method is simple and can be described as follows. First, it checks if k equals n; if this is the case, then we execute Algorithm 2 using parameters n,P, and d. For this particular case, it means the algorithm finds a feasible solution for the DCMST problem including all nodes of set V. On the opposite, if k<n, then we randomly choose an initial k-tuple of vertices from V in order to form the set As. Let Ac be the complement of As.Algorithm 3: VNS algorithm for the DCkMST problem using classical random local search strategy.
Data: an instance of the DCkMST problem using degree d=3.Result: a feasible solution and its objective function value.If (k=n) thenExecute Algorithm2 using parameters n,P,dSave feasible solution found and its objective function valueElseClassical random local search strategy:Generate an initial randomk-tuple of vertices. Let As denote this set of vertices and Ac its complement, P′=PAs,AsExecute Algorithm2 using parameters k,P′,dSave the initial feasible solution found and its objective function valueAsOP=As, AcOP=AcindK=1, cont=0, CpuTime=0While (CpuTime≤maxTime) doiter=iter+1Fori=1 to indKdoInterchange randomly an element ofAs with an element of AcP′=PAs,As, execute Algorithm 2 using parameters k,P′,dIf (a better solution is obtained) thenSave the new solution and setAsOP=As, AcOP=AciterOp=iter, cont=0, indK=1, Optime=Optime+CpuTime, CpuTime=0Elsecont=cont+1, As=AsOP, Ac=AcOPIf (cont≥1) thencont=0If (indK<indKMax) thenindK=indK+1ElseindK=1Return best feasible solution obtained and its objective function valueNext, we construct the submatrixP′ with the rows and columns of P corresponding to the nodes in As and execute Algorithm 2 using parameters k, P′, and d. This allows us to save an initial feasible solution together with its objective function value. Notice that the cardinality of set As is k, whilst the cardinality of set Ac is n−k. The current sets As and Ac are also saved in AsOP and AcOP, respectively. Subsequently, we perform the following steps iteratively, while the CpuTime variable is less than or equal to the maximum CPU time allowed which is controlled with parameter maxTime. Inside this loop, we randomly interchange an element of As with an element of Ac a predefined number of times which is controlled with parameter indK that is initially set to a value of one. Then, we check if a better solution is obtained. If so, we save the new solution found, its objective function value, update sets AsOP and AcOP, and reset the CpuTime variable to zero. Otherwise, we keep the previous best solution found and go back to the interchange step. Notice that parameter indK controls the number of swap moves the algorithm performs between sets As and Ac. It turns out that this step corresponds to a classical near-far random local search strategy related to VNS algorithm [28, 29]. The larger the number of swap moves, the wider the feasible space being explored. On the opposite, the smaller the number of moves, the closer the solutions which are being tested. In particular, in Algorithm 3, this is handled with indK variable which is increased by one unit each time a worse solution is obtained. This variable can be increased with up to a predefined value of indKMax−1. And we reset this parameter to a value of one each time a better solution is obtained in order to restart the search from a local to a wider feasible region. Finally, if indK equals indKMax, then we also reset indK to a value of one.


## 4.3.2. Embedded Q-Learning Strategy

Our Q-learning strategy is simple and consists of embedding a reinforcement learning approach in Algorithm3 instead of using classical random local search. Notice that the Q-learning approach was initially proposed in [31]. This approach mainly consists of an agent, a set of states S, and a set A of actions per state whereby performing an action a∈A; the agent moves from one state to another one according to a reward value, and the goal of the agent is to maximize its total future reward. The potential reward that an agent can obtain is computed as a weighted sum of the expected values of the rewards of all future steps starting from any of the current states. In order to embed this Q-learning approach in Algorithm 3, before entering the inner while loop of Algorithm 3, we need to initialize the matrix variable Q=QindKM,indKM with zero entries and also to declare the two scalar variables indKA and indicator.In particular, the value ofindicator variable, which can be true or false, is used to perform or not a piece of the code, whereas variables indK and indKA represent the number of swap moves from one to indKM to be performed between sets As and Ac. Thus, the underlying idea is to move from a particular number of swap moves to another one by using the aforementioned reinforcement learning strategy. Under this setting, the values of indKA and indK represent two consecutive states where an agent can be in different moments. For example, consider the case when indKA=2 and indK=8, then we should read it as an agent that is moving from state 2 to state 8. In our problem, it means we are first interchanging 2 elements and then 8 elements between sets As and Ac. Consequently, we can replace the inner code of the Else statement corresponding to the classical random local search strategy step of Algorithm 3 with the code of Algorithm 4. By doing so, we provide learning capabilities to Algorithm 3 in order to perform the random local search. In Algorithm 4, the parameters λ, μ, and Prob∈0;1 are the learning rate, discount factor, and probability of moving randomly from one state to another one [31]. Finally, the function expressions Rand, RandintindKM, and argmaxj∈1,…,indKMQindKA,j return a random fractional value between zero and one, a random integer value between one and indKM, and the jth column position with current maximum value in matrix Q while moving from state indKA, respectively.Algorithm 4: Embedded Q-learning code.
Generate an initial randomk-tuple of vertices. Let As denote this set of vertices and Ac its complement P′=PAs,AsExecute Algorithm2 using parameters k,P′,dSave the initial feasible solution found and its objective function valueAsOP=As, AcOP=AcindK=1, indKA=1, Q=zerosindKM,indKM, indicator=false, CpuTime=0While (CpuTime≤maxTime) doiter=iter+1Fori=1 to indKdoInterchange randomly an element ofAs with an element of AcSetP′=PAs,As and run Algorithm 2 using parameters k,P′,dIf (a better solution is obtained) thenSave the new solution and setAsOP=As, AcOP=Ac, indicator=trueiterOp=iter, Optime=Optime+CpuTime, CpuTime=0ElseAs=AsOP, Ac=AcOPIf (indicator) thenQindKA,indK=1−λQindKA,indK+λReward+μmaxj∈1,…,indKMQindK,jindKA=indKIf (Rand<Prob) thenindK=RandintindKMElseindK=argmaxj∈1,…,indKMQindKA,jindicator=false


## 4.4. ACO Algorithms

We now present the ACO algorithms. More precisely, first we present our classical ACO based strategy and then we present an adapted version of the generic ant-Q algorithm which is based on Q-learning strategy [26, 27]. Finally, we present a pure Q-learning-based algorithm. All these methods allow finding feasible solutions for the DCkMST problem while using a degree value of d=2.


### 4.4.1. Classical Ant Colony Optimization Strategy

Our first procedure is depicted in Algorithm5.Algorithm 5: Classical ACO algorithm for the DC kMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersρ, α, β, q and the set of ants (Ants)CpuTime=0, iter=0, LBest=∞, BestCPU=0τij=rand, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Δτija=0,∀i,j∈V, a∈AntsForEacha∈AntsdoConstruction of each ant tourRandomly assign an initial positionj∈V for the ant a according to graph G and construct a tour Ta using all remaining nodes of G.Find best solutionCompute the lengthLa of TaDenote byLakt the length of k consecutive nodes in Ta starting from node tComputeLak=minLak1,…,Lakn−k+1If (LBest>Lak) thenSetLBest=Lak and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Update accumulated pheromone matrixForEachi,j∈TadoΔτija=q/La, Δτjia=q/LaUpdate pheromone matrixForEachi,j∈Adoτi,j=1−ρτi,j+∑a∈AntsΔτijaReturn best feasible solution found and its objective function valueThe method requires an instance of the DCkMST problem. In step one, it first initializes all required parameters and variables. In particular, ρ, α, β, and q denote the pheromone evaporation coefficient, the influence of the pheromone, influence of visibility matrix, and a constant real value, respectively. Each entry of the pheromone matrix τ=τij for all i,j∈V represents the amount of pheromone deposited by the ants in the edge i,j for a state transition going from node i to j. Notice that each of these entries is initialized to a random value between zero and one, which is performed by using the Rand function. Similarly, the visibility matrix η=ηij for all i,j∈V represents the desirability value of edge i,j while choosing the state transition from i to j. In particular, this value is computed by ηij=1/Pij. The set of ants is denoted by Ants. The variables CpuTime, iter, LBest, and BestCPU denote the current CPU time, the number of iterations, and best current objective function value found so far and its exact CPU time in which this solution is obtained. Recall that, within ACO algorithms, each ant constructs a tour solution starting from a particular initial node while covering all nodes of the graph. For this purpose, each ant selects the next edge in its tour according to the following probability:(37)pija=τijαηijβ∑z∈Raτizαηizβ,where Ra denotes the unvisited set of nodes of ant a. As it can be observed, this probability depends on both the pheromone and visibility values [26].Subsequently, in step two of Algorithm5, we enter into a while loop in which we perform the following steps. First, we initialize to a value of zero each entry of matrix Δτ=Δτija, i,j∈V, a∈Ants, where each entry in this matrix represents the amount of pheromone deposited by ant a∈Ants. Next, for each ant, we construct a tour Ta starting from a randomly assigned node and compute its length La. Then, for each arc in Ta, we deposit the amount of pheromone of ant a. This value is computed by Δτija=q/La. Then, we compute the shortest tour of length k according to Ta and save it as the best tour found so far if its objective function value is less than LBest. For this purpose, we denote by Lakt the length of k consecutive nodes in Ta starting from node t. Notice that t=1,…,n−k+1. Also notice that the variable CpuTime is reset to a value of zero each time a better solution is obtained. This allows Algorithm 5 to run for another maxTime unit of time with the hope of finding better solutions. Next, we update each arc of the pheromone matrix τ according to parameter ρ and matrix Δτ. Finally, when the while loop is finished, we return the best feasible solution obtained together with its objective function value.


### 4.4.2. Ant-Q-Based Approach

Now, we present our adapted version of the generic Ant-Q algorithm which is based on Q-learning strategy [25, 27]. This procedure is depicted in Algorithm 6.Algorithm 6: Ant-Q algorithm for the DCkMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersα, β, q, Prob, λ, μ, and the set of ants (Ants)CpuTime=0, iter=0, LBest=∞, BestCPU=0AQi,j=0, Δi,j=0, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Construct ant toursForEacha∈AntsdoRandomly choose a nodej∈VTa=∅, Ta=Ta∪jFori=1 to n−1doForEacha∈AntsdoIfRand<ProbthenRandomly choose an unvisited node of the ant and add it toTa.Elsej=argmaxrallowedAQTai,rα1/PTai,rβTa=Ta∪jUpdate AQ-valuesAQTai,Tai+1=1−λAQTai,Tai+1+λμmaxjallowedAQTai+1,jAQTai+1,Tai=AQTai,Tai+1Find best solutionForEacha∈AntsdoDenote byLakt the length of k consecutive nodes in Ta starting from node tComputeLak=minLak1,…,Lakn−k+1If (LBest>Lak) thenSetLBest=Lak and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Reinforcement of ant toursForEachi,j∈TadoΔi,j=q/Lak, Δj,i=q/LakReinforcement of best solutionForEachi,j∈TkdoΔi,j=q/LBest, Δj,i=q/LBestDelayed reinforcement of AQ valuesForEacha∈AntsdoV=1,…,nFori=1 to n−1doV=V−Tai, V=V−Tai+1IfV≠∅thenAQTai,Tai+1=1−λAQTai,Tai+1+λΔTai,Tai+1+μmaxAQTai+1,VAQTai+1,Tai=AQTai,Tai+1Return best feasible solution found and its objective function valueSimilarly to Algorithm5, in step one of the Ant-Q Algorithm 6, we first initialize the required parameters and variables. In this case, the parameters α and β allow us to weigh the relative importance of learned AQ and visibility values, respectively. Parameter q is a nonnegative constant value, whilst parameters Prob, λ, and μ represent a probability value, a learning step, and a discount factor, respectively. Variables CpuTime, iter, LBest, BestCPU are analogously defined as for Algorithm 5. In particular, matrices AQ=AQi,j and Δ=Δi,j, for all i,j∈V, denote the Q-learning matrix and the delayed reinforcement matrices associated with the Ant-Q learning method [31].In step two of Algorithm6, we enter into a while loop and perform the following substeps. First, we construct a tour for each ant a∈Ants. This is performed iteratively while choosing a new unvisited node j randomly or according to the following expression:(38)j=argmaxrallowedAQTai,rα1PTai,rβ,where Tai denotes the current ith position of ant a in its tour Ta. Then, we update the learning matrix AQ according to all the edges of each ant tour and save the best solution obtained with cardinality k. Next, we proceed with reinforcement steps that are applied to all the edges belonging to each ant tour and best solution found so far [25, 27]. Finally, the algorithm returns the best feasible solution obtained and its objective function value.


## 4.4.1. Classical Ant Colony Optimization Strategy

Our first procedure is depicted in Algorithm5.Algorithm 5: Classical ACO algorithm for the DC kMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersρ, α, β, q and the set of ants (Ants)CpuTime=0, iter=0, LBest=∞, BestCPU=0τij=rand, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Δτija=0,∀i,j∈V, a∈AntsForEacha∈AntsdoConstruction of each ant tourRandomly assign an initial positionj∈V for the ant a according to graph G and construct a tour Ta using all remaining nodes of G.Find best solutionCompute the lengthLa of TaDenote byLakt the length of k consecutive nodes in Ta starting from node tComputeLak=minLak1,…,Lakn−k+1If (LBest>Lak) thenSetLBest=Lak and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Update accumulated pheromone matrixForEachi,j∈TadoΔτija=q/La, Δτjia=q/LaUpdate pheromone matrixForEachi,j∈Adoτi,j=1−ρτi,j+∑a∈AntsΔτijaReturn best feasible solution found and its objective function valueThe method requires an instance of the DCkMST problem. In step one, it first initializes all required parameters and variables. In particular, ρ, α, β, and q denote the pheromone evaporation coefficient, the influence of the pheromone, influence of visibility matrix, and a constant real value, respectively. Each entry of the pheromone matrix τ=τij for all i,j∈V represents the amount of pheromone deposited by the ants in the edge i,j for a state transition going from node i to j. Notice that each of these entries is initialized to a random value between zero and one, which is performed by using the Rand function. Similarly, the visibility matrix η=ηij for all i,j∈V represents the desirability value of edge i,j while choosing the state transition from i to j. In particular, this value is computed by ηij=1/Pij. The set of ants is denoted by Ants. The variables CpuTime, iter, LBest, and BestCPU denote the current CPU time, the number of iterations, and best current objective function value found so far and its exact CPU time in which this solution is obtained. Recall that, within ACO algorithms, each ant constructs a tour solution starting from a particular initial node while covering all nodes of the graph. For this purpose, each ant selects the next edge in its tour according to the following probability:(37)pija=τijαηijβ∑z∈Raτizαηizβ,where Ra denotes the unvisited set of nodes of ant a. As it can be observed, this probability depends on both the pheromone and visibility values [26].Subsequently, in step two of Algorithm5, we enter into a while loop in which we perform the following steps. First, we initialize to a value of zero each entry of matrix Δτ=Δτija, i,j∈V, a∈Ants, where each entry in this matrix represents the amount of pheromone deposited by ant a∈Ants. Next, for each ant, we construct a tour Ta starting from a randomly assigned node and compute its length La. Then, for each arc in Ta, we deposit the amount of pheromone of ant a. This value is computed by Δτija=q/La. Then, we compute the shortest tour of length k according to Ta and save it as the best tour found so far if its objective function value is less than LBest. For this purpose, we denote by Lakt the length of k consecutive nodes in Ta starting from node t. Notice that t=1,…,n−k+1. Also notice that the variable CpuTime is reset to a value of zero each time a better solution is obtained. This allows Algorithm 5 to run for another maxTime unit of time with the hope of finding better solutions. Next, we update each arc of the pheromone matrix τ according to parameter ρ and matrix Δτ. Finally, when the while loop is finished, we return the best feasible solution obtained together with its objective function value.


## 4.4.2. Ant-Q-Based Approach

Now, we present our adapted version of the generic Ant-Q algorithm which is based on Q-learning strategy [25, 27]. This procedure is depicted in Algorithm 6.Algorithm 6: Ant-Q algorithm for the DCkMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersα, β, q, Prob, λ, μ, and the set of ants (Ants)CpuTime=0, iter=0, LBest=∞, BestCPU=0AQi,j=0, Δi,j=0, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Construct ant toursForEacha∈AntsdoRandomly choose a nodej∈VTa=∅, Ta=Ta∪jFori=1 to n−1doForEacha∈AntsdoIfRand<ProbthenRandomly choose an unvisited node of the ant and add it toTa.Elsej=argmaxrallowedAQTai,rα1/PTai,rβTa=Ta∪jUpdate AQ-valuesAQTai,Tai+1=1−λAQTai,Tai+1+λμmaxjallowedAQTai+1,jAQTai+1,Tai=AQTai,Tai+1Find best solutionForEacha∈AntsdoDenote byLakt the length of k consecutive nodes in Ta starting from node tComputeLak=minLak1,…,Lakn−k+1If (LBest>Lak) thenSetLBest=Lak and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Reinforcement of ant toursForEachi,j∈TadoΔi,j=q/Lak, Δj,i=q/LakReinforcement of best solutionForEachi,j∈TkdoΔi,j=q/LBest, Δj,i=q/LBestDelayed reinforcement of AQ valuesForEacha∈AntsdoV=1,…,nFori=1 to n−1doV=V−Tai, V=V−Tai+1IfV≠∅thenAQTai,Tai+1=1−λAQTai,Tai+1+λΔTai,Tai+1+μmaxAQTai+1,VAQTai+1,Tai=AQTai,Tai+1Return best feasible solution found and its objective function valueSimilarly to Algorithm5, in step one of the Ant-Q Algorithm 6, we first initialize the required parameters and variables. In this case, the parameters α and β allow us to weigh the relative importance of learned AQ and visibility values, respectively. Parameter q is a nonnegative constant value, whilst parameters Prob, λ, and μ represent a probability value, a learning step, and a discount factor, respectively. Variables CpuTime, iter, LBest, BestCPU are analogously defined as for Algorithm 5. In particular, matrices AQ=AQi,j and Δ=Δi,j, for all i,j∈V, denote the Q-learning matrix and the delayed reinforcement matrices associated with the Ant-Q learning method [31].In step two of Algorithm6, we enter into a while loop and perform the following substeps. First, we construct a tour for each ant a∈Ants. This is performed iteratively while choosing a new unvisited node j randomly or according to the following expression:(38)j=argmaxrallowedAQTai,rα1PTai,rβ,where Tai denotes the current ith position of ant a in its tour Ta. Then, we update the learning matrix AQ according to all the edges of each ant tour and save the best solution obtained with cardinality k. Next, we proceed with reinforcement steps that are applied to all the edges belonging to each ant tour and best solution found so far [25, 27]. Finally, the algorithm returns the best feasible solution obtained and its objective function value.


## 4.5. A Pure Q-Learning-Based Algorithm

Now, we present a pure Q-Learning based approach that allows obtaining feasible solutions for the DCkMST problem while using d=2. This method is depicted in Algorithm 7, and it is described as follows.Algorithm 7: Pure Q-learning-based approach for the DCkMST problem.
Data: an instance of the DCkMST problem using degree d=2.Result: a feasible solution and its objective function value.Step 1Initialize parametersλ, μ, ProbCpuTime=0, iter=0, LBest=∞, BestCPU=0, Qi,j=0, ∀i,j∈VStep 2While (CpuTime≤maxTime) doiter=iter+1Construct a unique tourRandomly choose a nodej∈VT=∅, T=T∪j, V=V−jFori=1 to n−1doIfRand<ProbthenRandomly choose a nodeu∈VElseu=argmaxr∈VQj,rT=T∪u, V=V−uIfV≠∅thenQj,u=1−λQj,u+λ1/Pj,u+μmaxr∈VQu,rj=uFind best solutionDenote byLkt the length of k consecutive nodes in T starting from node tComputeLk=minLk1,…,Lkn−k+1If (LBest>Lk) thenSetLBest=Lk and save the best tour of length k in TkiterOp=iter, BestCPU=BestCPU+CpuTime, CpuTime=0Return best feasible solution found and its objective function valueSimilarly to Algorithm6, the first step of Algorithm 7 consists of initializing parameters λ, μ, and Prob, which represent the learning rate, discount factor, and a given probability value, respectively. We also define and initialize the required variables CpuTime, iter, LBest, BestCPU, and Q=Qi,j for all i,j∈V. These variables allow us to handle the current CPU time, the number of iterations, the best objective function value obtained so far, CPU time of the best current solution found, and the Q-learning matrix, respectively. Next, in step two of Algorithm 7, we iteratively construct a unique tour T with cardinality n and evaluate each subtour composed of k consecutive nodes in T. These steps are performed, while the current CPU time value is lower than the maximum value allowed which is denoted by maxTime. In order to obtain the best subtour of cardinality k from each tour T generated, we let Lkt denote the length of k consecutive nodes in T, starting from node t, and let Lk=minLk1,…,Lkn−k+1 be the minimum length value obtained. Notice that index t=1,…,n−k+1; otherwise, no subtour of cardinality k can be obtained from T. If a better solution is obtained, we save its length value in LBest, the subtour solution in Tk, the number of iterations in which this new solution is obtained, and reset the current CPU time value to a value of zero. The latter allows the algorithm to search for another maxTime unit of time. Also, notice that the construction of tour T requires to add at each step a new unvisited node u randomly or with maximum value Qj,u where j denotes the last node added to T. For this purpose, we remove node u from V at each step. The unvisited node u is chosen randomly if the Rand function generates a value in the interval 0;1 which is lower than Prob; otherwise, it is chosen according to the maximum value in the Q-learning matrix. Next, if set V is not empty, we update the corresponding entry in the Q-learning matrix as follows:(39)Qj,u=1−λQj,u+λ1Pj,u+μmaxr∈VQu,r.Notice that equation (39) is analog to the classical Q-learning algorithm [31]. In particular, the term 1/Pj,u represents the reward value obtained when going from node j to node u∈V. Finally, the algorithm returns the best feasible solution obtained and its objective function value.


## 5. Numerical Experiments

In this section, we conduct substantial numerical experiments in order to compare all the proposed models and algorithms. For this purpose, we implement Matlab programs using CPLEX 12.7 solver [54] to solve the MILP and LP models. The numerical experiments have been carried out on an Intel(R) 64 bits core (TM) with 3 GHz and 8 G of RAM under Windows 10. CPLEX solver is used with default options. We consider complete graph instances with random uniform and Euclidean distance costs. More precisely, we generate five instances with dimensions of n=50,100,200,300,400 nodes where each entry in the symmetric matrix P=Pij, i,j∈A is randomly drawn from the interval 0;100. Then, we further consider four Euclidean benchmark instances from TSPLIB [34] referred to as “Berlin52,” “gr96,” “ch150,” and “a280” with dimensions of n=52,96,150,280 nodes, respectively. All these instances are solved for a different number of active nodes k≤n ranging from k=10 and up to k=400 nodes for the instances with random costs and up to k=280 for the Euclidean ones.


### 5.1. Numerical Results for the MILP Models

In Tables1–4, we present numerical results for P1, P1h, P2, and P2h. Notice that all these tables present the same column information. More precisely, column 1 shows the instance number, whilst columns 2 and 3 present the number of nodes of each graph instance and the value of k, respectively. Next, in columns 4–8 and 9–13, we present the optimal or best solution obtained with each model in at most 2 hours of CPU time, number of branch and bound nodes, CPU time in seconds, the optimal solution of each LP relaxation, and its CPU time in seconds. Finally, columns 14 and 15 present gaps that we compute by Opt−LP/Opt∗100 where Opt and LP refer to the optimal solution found with the MILP and LP models, respectively. Also, notice that each row in Tables 1 and 3 and in Tables 2 and 4 corresponds to a same instance.Table 1
Numerical results obtained withP1 and P1h for random graphs using d=3 and d=2.
#nkP1P1hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap1 %Gap1h%d=31501046.86750.8043.380.1746.861091.0643.380.197.437.43220113.974980.75109.460.20113.973620.73109.460.193.963.96330204.3514702.79199.050.19204.3519113.48199.050.172.592.59450459.0910051.76453.490.19459.0929015.74453.490.171.221.22510025115.44231.95114.630.48115.44251.79114.630.420.690.69650263.9525018.30262.040.47263.9518726.43262.040.420.730.73775454.4431865.04451.580.42454.44171524.23451.580.420.630.638100749.7536160202.21745.090.45749.7565193338.15745.090.440.620.62920050162.55103510.39160.881.73162.55146210.86160.881.641.031.0310100378.4380814506.06374.871.67378.4394564518.29374.871.650.940.9411150652.8629440936417.84648.891.73652.8626357647200648.891.950.610.61122001049.2540385372001044.742.261049.2521541572001044.741.780.430.43d=21501052.22491.5846.300.1752.221290.6146.300.1711.3411.34220134.082372.06129.710.22134.089042.57129.710.223.263.26330240.30451.26239.320.20240.3091.26239.320.200.410.41450564.3729487.38557.080.19564.3719537.47557.080.171.291.29510025123.714044.63121.540.48123.711483.99121.540.451.751.75650287.78109337.88285.310.61287.7818757.75285.310.500.860.86775529.5601.92529.310.50529.5613.99529.310.480.050.058100915.296839.75914.630.48915.29778.14914.630.450.070.07920050184.29012.28184.292.11184.2907.02184.292.030010100437.08642499.08435.452.50437.08256893.76435.452.280.370.3711150775.8927594.08774.762.40775.899042.20774.762.200.150.15122001285.3219067.751284.582.231285.3245998.281284.581.890.060.06Table 2
Numerical results obtained withP1 and P1h for Euclidean graphs using d=3 and d=2.
#nkP1P1hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap1 %Gap1h %d=31Berlin5210274.46102.25253.210.34274.46600.94253.210.287.747.74220926.77543991218.23801.280.20926.77591996300.65801.280.1913.5413.543301935.71424693336001666.190.201935.71529283236001666.190.1913.9213.924526081.63202461636005567.130.206081.63168860136005567.130.198.468.465gr962574.14488206360040.180.4174.14452709360040.180.3945.8145.81650163.115132933600117.660.55162.804629023600117.660.4227.8627.73775269.155935063600222.470.48269.236727753600222.470.4117.3417.37896436.6511923673600384.510.80436.6911687843600384.510.7011.9411.959ch150501607.2424334236001258.470.951610.9315735336001258.470.9421.7021.8810802657.2118647436002303.051.012694.8815476836002303.050.9013.3314.54111003457.8619991436003058.710.973441.5912931836003058.710.9211.5411.13121505914.1521578736005522.970.955923.7823488736005522.970.896.616.77d=21Berlin5210279.2300.52261.790.30279.2300.53261.790.226.256.252201015.3922537.08862.040.251015.39873515.91862.040.2215.1015.103302065.12944612.211871.760.282065.125439.361871.760.279.369.364526968.7716998.246545.080.276968.7712447.446545.080.236.086.085gr962574.59380670360045.190.7774.59324849360045.192.3939.4239.42650168.532704153600128.730.92174.251947523600128.738.4523.6226.12775283.141548383600248.4413.60288.493104293600248.440.5812.2513.88896474.951685393600452.149.66476.083975403600452.140.724.805.039ch150501753.839874236001356.911.041807.838631036001356.919.9122.6324.9410802917.727382936002489.252.292881.726021036002489.2510.6414.6913.62111003744.476671036003302.5415.884013.585451836003302.5414.0511.8017.72121506411.736274336006158.671.156397.127057436006158.671.013.953.73Table 3
Numerical results obtained withP2 and P2h for random graphs using d=3 and d=2.
#nkP2P2hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap2 %Gap2h %d=31501046.86422.2943.380.2346.863401.2543.380.207.437.43220113.97216116.65109.460.25113.9799911.42109.460.273.963.96330204.351161.67197.610.19204.35202414.65197.610.193.303.30450459.091231.00444.770.25459.096293.87447.870.203.122.44510025115.4435846.71114.540.53115.449754.60114.540.700.780.78650263.95124373.63261.020.51263.9566877.50261.020.751.111.11775454.44459143.16449.510.55454.44525111.53449.510.671.081.088100749.752422.56728.702.29749.75040.36734.840.592.811.99920050162.554793995.62160.362.21162.559202835.21160.363.291.351.3510100378.434792671.16373.283.81——7200373.283.681.36—11150652.864792097.25642.844.71652.864651663.29642.843.491.531.53122001049.25607505.271030.0430.111049.255021490.461034.233.211.831.43d=21501052.2251.0146.420.2552.221901.7346.300.2011.1211.34220134.083793.43129.900.33134.08151521.95129.710.203.123.26330240.302357.35239.320.28240.3073619.08239.320.200.410.41450564.3748820.20557.480.28564.37164742.00557.080.191.221.29510025123.711067248.12121.551.50123.711933237.10121.540.701.741.75650287.781676632.83285.352.84287.78857356.21285.310.810.840.86775529.5608.88529.313.23529.561050639.40529.310.760.050.058100——7200914.604.13——7200914.630.89——920050184.2910228.02184.2913.13——7200184.293.810—10100——7200435.4526.49——7200435.456.61——11150——7200774.7655.04——7200774.766.62——122001285.3210451408.711284.5894.68——72001284.584.150.06——: no solution found.Table 4
Numerical results obtained withP2 and P2h for Euclidean graphs using d=3 and d=2.
#nkP2P2hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap2 %Gap2h %d=31Berlin5210274.4600.97253.210.21274.4600.82253.210.217.747.74220926.7702.25697.510.25926.7702.15697.510.2124.7424.743301935.711752.311414.090.221935.714276.121414.090.2026.9526.954526081.635562.004822.430.286081.63952.204923.770.1920.7119.045gr962572.402518268.8738.470.5172.40913301.4238.470.6146.8746.87650162.803274316.75110.480.56162.801816250.78110.460.6032.1432.15775267.57591109.28210.191.86267.57615157.00210.180.6821.4421.45896436.2348070.35354.603.72436.2398542.30357.830.6018.7117.979ch150501567.771734921.271219.691.21——36001219.691.4622.20—1080——36002237.221.25——36002237.221.85——11100——36002974.081.603399.1710371009.762974.081.97—12.51121505882.45982212.275339.1827.695882.45974245.395353.541.619.248.99d=21Berlin5210279.2302.95262.400.25279.2300.67261.790.206.036.252201015.39148721.75868.220.371015.39843.21862.040.1914.4915.103302065.1204.551874.560.342065.1202.871871.760.229.239.364526968.7797231.656554.140.286968.7702.896545.080.205.956.085gr962573.561592845.8245.410.7773.562525621.6445.190.6238.2738.57650168.53862299.41129.002.81168.53601350.55128.730.6723.4623.62775283.141235408.99248.894.12283.141672447.78248.440.6412.1012.25896474.952874329.41456.192.50474.952482315.00452.140.563.954.809ch15050——36001358.7316.52——36001356.912.46——10803167.66124136002491.5222.79——36002489.252.3421.35—111005425.41122436003306.7420.47——36003302.542.7939.05—121507176.67369836006166.4929.53——36006158.671.7814.08——: no solution found.We limit CPLEX to 2 hours of CPU time for the random input graphs, whilst for the Euclidean ones, we limit CPLEX to 1 hour in order to avoid CPLEX shortage of memory events. Consequently, each reported solution is optimal if its CPU time is less than 2 hours. Otherwise, it corresponds to the best solution obtained with CPLEX in at most 2 hours. Subsequently, in Tables5 and 6, we present numerical results for both P3 and P3h. These two tables also present the same column information. In particular, columns 1–3 present the same information as in Tables 1–4. Next, in columns 4–8 and 9–13, we report the optimal or best solution obtained with Algorithm 1 in at most 1 hour of CPU time, number of branch and bound nodes, CPU time in seconds, number of cycles added to each model, and number of iterations as well. The number of branch and bound nodes, in this case, corresponds to the sum of all branched nodes within each iteration of Algorithm 1. Notice that Theorem 2 in [24] ensures that the solutions obtained within each iteration of Algorithm 1 are lower bounds for the optimal solution of the problem. Consequently, an optimal solution is obtained when the corresponding CPU time is less than 1 hour. Notice that obtaining tight lower bounds is of crucial importance when developing exact methods as it allows to speed up the process significantly. Further notice that models P1, P1h, P2, and P2h obtain upper bounds for the problem if the optimal solution cannot be reached within 2 hours. Thus, we provide an interval where the optimal solution lies.Table 5
Lower bounds obtained withP3 and P3h for random graphs using d=3 and d=2.
#nkP3P3hOptB&BnTime#Cycles#IterOptB&BnTime#Cycles#Iterd=31501046.86531.753446.86362.1534220113.9701.3323113.9721.0123330204.3510215.621414204.357555.631414450459.099556.402322459.096115.992322510025115.44534.0123115.44222.9023650263.95337324.621920263.95339426.101920775454.44715458.034344454.44755859.7543448100749.75214640441.97211211749.75157118446.80211211920050162.55823.4545162.551523.784510100378.2414387063604.81252252378.2012412773624.1123723711150651.4915097493626.19263263651.4312140613614.95239239122001046.968528083602.983353341046.9510365043615.80320319d=21501052.2260.971252.2200.8712220134.08111.6144134.08171.5644330240.30261.3112240.30101.0812450564.3700.5232564.3700.4732510025123.711663.5612123.71362.9512650287.784997.1965287.782226.0465775529.5600.9001529.5600.84018100915.2904.0745915.2903.3445920050184.2904.3201184.2904.040110100437.08327431.7344437.08107530.844411150775.89018.5312775.89019.2012122001285.3207.05121285.3206.4312Table 6
Lower bounds obtained withP3 and P3h for Euclidean graphs using d=3 and d=2.
#nkP3P3hOptB&BnTime#Cycles#IterOptB&BnTime#Cycles#Iterd=31Berlin5210274.4601.4444274.4621.1244220926.772858450.958180926.774293162.0981803301935.71432870650.352982941935.71384635710.192982944526081.6315228434720.069379326076.5712618693611.008608545gr962568.1427826393619.001699268.0725079133621.4816184650154.0618782433716.3219665153.8917406253616.4219062775265.1111887453614.37297152264.8611130793626.85274129896427.944341603602.41734638427.943802493606.057336379ch150501545.7514873463642.21172791545.2313079163697.271677510802599.8014517303659.13212772598.4013279313776.4720272111003388.098590903682.262261083387.087975713617.0121698121505882.451334622909.894994195882.451270322874.88499419d=21Berlin5210279.2300.5122279.2300.67222201015.39642.471551015.39622.091553302065.121342.231542065.12632.421544526968.7700.731036968.7700.751035gr962573.566891771525.831255173.567645551832.4012551650168.53189277981.3210936168.533036961584.5610936775283.1424803100.756318283.141588096.286318896474.951195123.933410474.95480217.5834109ch150501637.179042113620.84112351636.059899404001.451063110802788.875649183790.60127322786.685396993697.0312330111003627.9328619279.8866143627.9336800341.066614121506368.79815243.7138106368.79904645.693810From Table1, first we observe that both P1 and P1h allow obtaining the optimal solution of the problem for most of the instances in less than 2 hours. In particular, for the degree value of d=3, we see that P1 and P1h cannot solve the instances #12 and #11-#12, respectively, whilst for d=2, all the instances are solved to optimality. Regarding the CPU times, for d=3, we see that P1h requires a higher amount of CPU time, whilst for d=2, we obtain similar values with both models. Next, we observe that the number of branch and bound nodes is slightly lower for P1 than for P1h. We also see that the LP bounds are near-optimal and exactly the same for both models which are confirmed by the gaps. Notice that all the LP relaxations are solved in less than 3 seconds. Finally, we observe that the CPU time values required to solve both models increase with k.Regarding the Euclidean graph instances presented in Table2, we observe that most of them cannot be solved to optimality in one hour. We mention that these instances proved to be very difficult to solve, and this is the main reason we limit CPLEX to one hour in order to avoid CPLEX shortages of memory. The difficulty in solving these instances is also reflected in the number of branch and bound nodes and in the gaps obtained which are significantly higher compared to Table 1. Notice that these gaps decrease with k in which evidences solving the DCkMST problem is harder to solve than the classical DCMST problem, at least for these instances. Next, we observe that the CPU time values of the LP relaxations are higher than those reported in Table 1 and that the degree values do not seem to affect the performance of each model. We further see that the optimal or best objective function values obtained for the instance “ch150” are significantly higher for d=2. Finally, we observe that the objective function values obtained with P1 for the instances #7–#10, #12 and #6–#9, and #11 are lower than those obtained with P1h for d=3 and d=2, respectively. However, the opposite occurs for the instances #6, #11 and #10, #12 while using d=3 and d=2, respectively.From Table3, we mainly observe that, for d=3, only P2 allows obtaining the optimal solution of the problem for all the instances in less than 2 hours. Notice that P2h only failed to find a feasible solution for the instance #10, whilst for d=2, neither P2 nor P2h can solve all the instances to optimality. However, P2 can still solve more instances than P2h. Regarding the number of branch and bound nodes and LP times, we observe smaller and slightly higher values than in Table 1, respectively. Finally, we observe that the LP bounds are near-optimal which is confirmed by the gaps. In Table 4, we observe that P2 and P2h can solve more instances to optimality than P1 and P1h in Table 2, respectively. However, in this case, both P2 and P2h cannot find a feasible solution for some of the instances. Next, we observe smaller values for the branch and bound nodes and CPU times with similar orders of magnitude for the LP relaxations when compared to Table 2. Finally, we observe that the LP bounds are not tight which is again confirmed by the gap columns.From Table5, we mainly observe that models P3 and P3h, which are both solved with Algorithm 1, allow to solve to optimality the instances #1–#9 and #1–#12 for d=3 and d=2, respectively. In particular, we see that, for d=3, it is harder to solve these instances than for d=2. This fact is reflected in the number of branch and bound nodes, CPU times of the LP relaxations, number of iterations, and number of cycles added to each exponential model which are clearly lower when d=2. Notice that, for the instances #10–#12 and degree value of d=3, we only report the best objective function values obtained which are in fact lower bounds. Regarding the Euclidean instances reported in Table 6, we observe that the instances #1–#3, #12 and #1–#8, #11-#12 are all solved to optimality for d=3 and d=2, respectively. We also see in Table 6 that it is harder to solve the instances for d=3 than for d=2. Notice that, for d=2, we almost solve all the instances to optimality with the exception of instances #9 and #10 for which we obtain lower bounds. Again, this observation can be verified by looking at the number of branch and bound nodes, CPU times of the LP relaxations, number of iterations, and number of cycles added to each proposed model which are smaller when d=2. Finally, from the numerical results presented in Tables 1–6, we can conclude that model P1 outperforms the other ones as it allows to obtain either an optimal or a best upper bound for most of the tested instances. Notice that the flow and exponential models also have good performance in terms of optimality, but for the large-size instances, the flow models deteriorate rapidly, whilst the exponential ones cannot be handled efficiently in terms of CPU times. Consequently, in Tables 7 and 8, we present numerical results for P1 for large random and Euclidean input graph instances with up to 400 and 280 nodes, respectively. In particular, in Table 7, these numerical results are reported for d=3, whereas in Table 8, these numerical results are obtained while using a degree value of d=2.Table 7
Numerical results obtained withP1 for large input graphs using d=3.
#nkP1OptB&BnTimeLPTimeGap1%Random graph instances130050118.8212040.30117.394.401.202100283.63371964.73281.704.200.683200703.429949132022.93701.154.210.3243001334.6117753872001332.445.010.16540050110.201651193.66108.488.471.566100236.59100882316.85234.519.480.887200522.6113577259.07520.668.240.3784001459.533176372001456.3410.840.22Euclidean graph instances1a28050480.35232223600384.004.2120.0621001011.07429573600784.004.5222.4631501316.172876136001184.164.5710.0342802493.684623236002385.763.854.33Table 8
Numerical results obtained withP1 for large input graphs using d=2.
#nkP1OptB&BnTimeLPTimeGap1%Random graph instances130050132.92034.03132.876.910.042100321.64454309.74320.828.300.263200820.70272520.03820.0212.600.0843001623.8512501301.401623.1611.510.04540050123.024132003.07120.0519.052.426100265.68305873073.49262.4419.191.227200598.2025272448.10595.0834.010.5284001818.0399972001797.7616.211.12Euclidean graph instances1a28050534.97136853600384.005.4828.2221001133.49232123600784.004.9330.8331501417.201140536001186.716.1616.2642002991.961103636002513.314.6616.00From Tables7 and 8, we observe similar trends as for the above Tables 1–4. More precisely, we observe that P1 can solve almost all random input graph instances to optimality in less than 2 hours using both degree values. Only the instance number #8 could not be solved to optimality although an upper bound is reported for this particular instance. In contrast, none of the Euclidean instances can be solved to optimality in less than 1 hour of CPU time. For these instances, only upper bounds are reported too. The difficulty in solving the Euclidean instances can also be observed in the gap columns which report significantly higher values as a consequence of the LP bounds obtained. On the opposite, the LP bounds obtained for the random graph instances are near-optimal. Finally, we observe from both Tables 7 and 8 that the number of branch and bound nodes is significantly smaller for d=2 than for d=3 and that the LP times are slightly higher for d=2.


### 5.2. Numerical Results for VNS Algorithms

In Tables9 and 10, we present numerical results obtained with the proposed VNS algorithms for both random and Euclidean input graph instances while using a degree value of d=3. Recall that our second VNS approach consists of replacing the code lines of Algorithm 3 corresponding to the random local search strategy with the code of Algorithm 4 which represents the embedded Q-learning strategy. Hereafter, we denote by VNSR and VNSQ these two VNS approaches. In Algorithms 3 and 4, we arbitrarily set the maximum CPU time allowed to a value of maxTime=250 seconds, whereas the input parameters indKMax and indKM are set to a value of 10. Similarly, the required parameters of Algorithm 4 were calibrated to λ=0.25, μ=0.25, and Prob=0.15. Finally, the reward value and parameter θ of Algorithm 2 were set to Reward=1 and θ=0.1, respectively.Table 9
Numerical results obtained with VNS algorithms for random and Euclidean instances usingd=3.
#nkMILPVNSRVNSQBestTime (s)VNSiniVNSTime (s)#IterGapRVNSiniVNSTime (s)#IterGapQInput graphs with random uniform costs1501046.860.80254.7648.276.06138423.00208.7346.860.6216200220113.970.73266.68118.5322.8072134.00353.48116.078.2927961.84330204.351.67338.03206.3546.5136740.97326.75205.4330.4925030.52450459.091.00463.14463.140.0810.88463.14463.140.0810.88510025115.441.79352.14115.44115.10238260413.40118.7718.5988542.88650263.956.43564.30269.1861.8929621.98517.86265.4286.4173820.55775454.445.04600.65465.48197.6126812.42659.23456.59215.4640980.478100749.7522.56749.87749.870.2510.01749.87749.870.2410.01920050162.5510.39555.42179.981857.836364010.72463.99169.38477.84187284.2010100378.43506.06725.35395.982928.53123244.63763.28386.121657.5291942.0311150652.861663.29883.28721.071250.9191210.44926.60670.361130.828722.68122001049.25505.271056.341056.343.9810.671056.341056.343.9010.67Input graphs with Euclidean distance costs1Berlin5210274.460.821976.48274.461.41595901645.70274.460.094700220926.772.153923.04926.772.69738203161.01926.770.3682003301935.712.314961.302104.325.2951338.714484.232071.97178.341930487.034526081.6326081.636081.630.01106081.636081.630.00105gr962572.40268.87204.0079.28122.542303369.50247.0774.14271.895085832.40650162.80250.78307.92179.5074.783409410.25260.74163.702.7411630.55775267.57109.28388.13269.8820.9730490.86372.46267.57141.16211760896436.2342.30436.23436.230.0210436.23436.230.01109ch150501567.77921.273441.771716.09348.671547349.463290.911696.0637.20161078.1810802657.2136004327.972792.69115.65138895.094480.372771.08223.80216754.28111003399.171009.764722.643527.70408.36153663.784745.733499.01288.25138582.93121505882.45212.275882.455882.450.12105882.455882.450.1210Table 10
Numerical results obtained with VNS algorithms for the instances of Table7 using d=3.
#MILPVNSRVNSQBestTime (s)VNSiniVNSTime (s)#IterGapRVNSiniVNSTime (s)#IterGapQInput graphs with random uniform costs1118.8240.30534.32139.361330.834742917.28510.03127.26494.26362917.102283.6364.73762.97338.442549.43766519.32775.47298.831893.4284815.353703.422022.931088.49850.681381.5532620.931000.36750.053851.549666.6241334.6172001349.381349.3832.6811.101349.381349.3832.7011.105110.20193.66492.26126.271376.585215614.58526.08124.871015.616372013.316236.592316.85783.62306.031465.56824229.35677.53257.463628.70202968.827522.61259.071129.36696.202587.2766433.211058.37635.214130.04104521.5481459.5372001485.591485.5983.4511.781485.591485.5983.1311.78Input graphs with Euclidean distance costs1480.3536001095.00446.91314.4196332−6.961076.54440.3874.3931647−8.3221011.0736001385.72897.55425.8011551−11.221447.17890.15195.388251−11.9531316.1736001774.881304.78786.148339−0.861689.481303.48570.315417−0.9642493.6836002440.942440.947.941−2.112440.942440.947.941−2.11In particular, in Table9, we report numerical results for the instances presented in Tables 1–6, whereas in Table 10, we report numerical results for the large-size instances presented in Table 7. In both tables, columns 4–15 and 2–13 contain exactly the same column information, respectively. More precisely, in Table 10, columns 1 to 3 present the instance number, number of nodes, and the value of k. Next, in columns 4 and 5, we report the minimum objective function and CPU time values reported for each instance in Tables 1–4. Subsequently, in columns 6–10 and 11–15, we report for each VNS approach, the initial and minimum objective function values, CPU time in seconds required by VNS, number of iterations, and gaps, respectively. The gaps are computed by VNS−Best/Best∗100.From both Tables9 and 10, first we observe that the objective function values of the initial solutions are significantly higher than the best ones. The latter clearly evidences the effectiveness of both VNS approaches. Notice that when k=n, the problem reduces to the classical DCMST problem. Consequently, for these instances, the solutions are obtained only with Algorithm 2. Next, we observe that the CPU times required by both VNS methods are larger for the random instances than for the Euclidean ones. From Table 9, next, we observe that the solutions obtained with VNSQ outperform those obtained with VNSR for both random and Euclidean instances. Although, in general, we see that both VNS procedures allow obtaining near-optimal solutions which is confirmed by the gap columns. On the opposite, we see that the objective function values reported in Table 10 are not tight for the random instances when compared to those obtained by the MILP models. But still, in this case, the gaps reported for VNSQ are significantly better than for VNSR. However, for the Euclidean instances reported in Table 10, we see that both VNS approaches allow us to obtain better solutions than the MILP models. The latter can be verified by the negative gaps which show that the objective function values obtained with VNS algorithms are significantly lower than those obtained with the MILP models.


### 5.3. Numerical Results for ACO Algorithms

Now, we report numerical results obtained with Algorithms5 and 6 for random and Euclidean instances while using a degree value of d=2. More precisely, in Table 11, we report numerical results for the instances presented in Tables 1–6. However, in Table 12, we report numerical results for the large-size instances presented in Table 8. Hereafter, we denote by ACOR and ACOQ the two ant colony optimization approaches presented in Algorithms 5 and 6, respectively. Recall that the first one is constructed based on the classical ACO metaheuristic [26]. However, the latter is based on the Ant-Q algorithm proposed in [25, 27]. In Algorithm 5, we calibrated the input parameters to α=0.5, β=2, ρ=0.9, and q=1. However, in Algorithm 6, we calibrated the input parameters to α=1, β=4, q=1, Prob=0.1, λ=0.15, and μ=0.15. Finally, in both Algorithms 5 and 6, we arbitrarily set the maximum CPU time and number of ants to maxTime=250 seconds and Ants=20, respectively.Table 11
Numerical results obtained with ACO algorithms for random and Euclidean instances usingd=2.
#nkMILPACORACOQBestTime (s)ACOiniACOTime (s)#IterGapRACOiniACOTime (s)#IterGapQInput graphs with random uniform costs1501052.220.6193.6952.2254.7855760259.2652.2289.8920830220134.081.61175.53135.7524.5825511.24955.85135.75126.1721111.24330240.301.08370.28244.3126.0126481.66435.68245.24234.2544372.05450564.370.47684.50584.185.535613.512486.65588.3148.698034.24510025123.712.95163.55123.712.066701231.26123.71192.919950650287.786.04337.69293.9029.209752.122546.91298.88339.6617313.85775529.560.84650.54550.4146.0015333.933804.55563.89220.2811646.488100915.293.341093.20970.4350.6616996.025188.971021.63239.28123011.61920050184.294.04279.89193.34168.4716264.912495.55198.36698.7110937.6310100437.0830.84546.76467.7668.136577.014856.80492.90258.6340812.7711150775.8918.53941.43836.899.54907.867444.72870.56273.5346012.20122001285.326.431734.021480.7615.9115615.2010277.041614.61321.3052925.61Input graphs with Euclidean distance costs1Berlin5210279.230.512195.31279.23389.574209604581.09279.2321.4435102201015.392.091270.761015.392.0120109757.741020.0272.2411890.453302065.122.232757.232065.1258.286258012888.362069.9032.445350.234526968.770.739738.247156.6517.7118952.6922432.487234.276.471063.805gr962573.56621.64112.7874.59144.6756471.40144.7274.807.65451.68650168.53299.41219.39172.3897.1838112.28546.20180.29161.809846.97775283.1496.28379.18292.2024.979663.19609.61302.6650.223026.89896474.9517.58619.16497.0666.0225784.65944.82534.30320.22195612.499ch150501753.8336002050.391678.9338.58672−4.2718522.681722.83222.95613−1.7610802881.7236003667.842843.87169.472972−1.3127659.442929.44141.143941.65111003627.93279.884833.863788.2698.1817184.4135594.853773.76325.019084.01121506368.7943.717566.166720.3142.557425.5152740.447005.90113.5031510.00Table 12
Numerical results obtained with ACO algorithms for random and Euclidean instances of Table8 using d=2.
#MILPACORACOQBestTime (s)ACOiniACOTime (s)#IterGapRACOiniACOTime (s)#IterGapQInput graphs with random uniform costs1132.9234.03193.95146.5326.599810.232228.39158.671443.44109119.372321.64309.74419.46367.1643.3016514.154620.75396.94589.5246823.413820.70520.03980.56898.2564.562379.449596.31923.88918.3670312.5741623.851301.402122.151933.14128.8249519.0414537.582107.741130.5586529.795123.022003.07172.02136.3398.3021010.812581.64143.352080.9391716.526265.683073.49328.39298.32155.4733812.285243.13316.621187.5353019.177598.202448.10754.03683.366.101314.2310356.11724.671247.0556321.1481818.0372002370.222224.405.611122.3520229.562631.201151.9052144.72Input graphs with Euclidean distance costs1534.973600426.30401.87113.06632−24.87675.36405.20227.45196−24.2521133.493600990.04846.6796.56539−25.302944.45919.303.483−18.8931417.2036001554.941311.6242.31236−7.442909.141358.90312.40270−4.1142991.9636003477.832983.922.1811−0.266137.083239.96292.152558.28The legend of Table11 is as follows. From columns 1 to 3, we present the instance number, the number of nodes of the input graph, and the value of k. Next, in columns 4 and 5, we present the minimum objective function and CPU time values reported in Tables 1–6. We repeat this information for the sake of clarity. Next, in columns 6–10 and 11–15, we report the objective function values of the initial solutions, the objective function values of the best solutions, CPU time in seconds, number of iterations, and gaps which are computed by ACO−Best/Best∗100, respectively. Similarly, the first three columns of Table 12 report the instance number and the best objective function values and CPU times reported in Table 8. Again, this information is repeated for comparison purposes. Finally, the legends of columns 4–8 and 9–13 are exactly the same as for the columns 6–10 and 11–15 in Table 11.From Table11, we observe that the best solutions obtained with both ACO approaches are near-optimal and far from the initial solutions obtained. These facts prove the effectiveness of the proposed methods. Next, we see that the CPU times are lower for ACOR than for ACOQ for most of the instances and in particular for the larger ones. Regarding the number of iterations, in general, we observe similar orders of magnitude for both methods. Finally, we observe that the gaps obtained with ACOR are tighter than those obtained with ACOQ for both the random and Euclidean instances. In particular, for large instances, these values are significantly better. Finally, notice that we obtain negative gaps for some of the large Euclidean instances. This clearly evidences that the solutions obtained with the ACO methods outperform significantly the solutions obtained with the MILP models.From Table12, we mainly observe that the solutions obtained with both ACO algorithms are significantly worse than those reported in Table 11 for random instances. On the opposite, for all the Euclidean instances, the solutions obtained with the ACO approaches are significantly better than those obtained with P1 which is again confirmed by the negative gaps obtained.


### 5.4. Numerical Results for QL Algorithm

In Tables13 and 14, we present numerical results obtained with Algorithm 7 (denoted as QL) for both random and Euclidean instances while using a degree value of d=2. In Tables 13 and 14, we report numerical results for the same instances presented in Tables 1–6 and for the instances of Table 8, respectively. The legend of Table 13 is as follows: columns 1–3 present the instance number, number of nodes (or name of the Euclidean instance), and the value of k, respectively. Notice that the name of each Euclidean instance indicates at the end the number of nodes it contains. For example, the instance name “Berlin52” has 52 nodes. Next, columns 4 and 5 report the best solution and minimum CPU time in seconds obtained with the MILP models. Finally, columns 6–10 report the initial solution obtained with Algorithm 7, its best solution found, CPU time in seconds, number of iterations, and gaps obtained which are computed by QL−Best/Best∗100, respectively. The legend of Table 14 is exactly the same as for Table 13 and is obtained by removing columns 2 and 3 from Table 13.Table 13
Numerical results obtained with Q-learning algorithm for random and Euclidean instances usingd=2.
#nkMILPQLBestTime (s)QLiniQLTime (s)#IterGapInput graphs with random uniform costs1501052.220.61427.2652.229.03423590220134.081.61794.00138.008.76406222.92330240.301.081473.09245.622.25100792.21450564.370.472567.98601.3715.38709226.55510025123.712.951207.44123.714.7598900650287.786.042440.56291.8998.122118911.42775529.560.843473.78565.78487.9010678996.838100915.293.344975.561024.7730.106434511.96920050184.294.042543.82199.27130.441226718.1210100437.0830.845027.12476.04123.871181058.9111150775.8918.537342.28843.6038.57365128.72122001285.326.439799.771546.62227.2121410720.32Input graphs with Euclidean distance costs1Berlin5210279.230.512692.71279.2339.9320032202201015.392.0912651.131020.02151.287743860.453302065.122.2316299.432070.6545.422283630.264526968.770.7330795.847234.2785.364339013.805gr962573.56621.64858.8074.5978.731994191.40650168.53299.411774.49180.0589.152251946.83775283.1496.282500.70303.35260.646483027.13896474.9517.583352.78552.06191.6847743016.239ch150501753.83360018911.281711.33125.94186763−2.4210802881.72360028695.922949.71202.883061202.35111003627.93279.8837580.423953.79107.611606248.98121506368.7943.7154468.427122.77565.4484402011.83Table 14
Numerical results obtained with Q-learning algorithm for the random and Euclidean instances of Table8 using d = 2.
#MILPQLBestTime (s)QLiniQLTime (s)#IterGapInput graphs with random uniform costs1132.9234.032455.40147.03116.827292410.612321.64309.745203.22347.51101.90630218.043820.70520.0310356.83897.52160.04998319.3641623.851301.4015547.511990.10441.3826392222.555123.022003.072173.05137.42155.636459711.706265.683073.495346.13307.22442.6018313015.637598.202448.1010535.56691.97286.7611571715.6781818.03720020088.172287.40305.8012663125.81Input graphs with Euclidean distance costs1534.9736006193.93426.46170.69119614−20.2821133.49360011619.78909.0382.0557482−19.8031417.20360018065.401441.78275.341924371.7342991.96360033966.222800.73130.9891758−6.39From Table13, first we observe that the initial solutions obtained with Algorithm 7 are considerably worse than those obtained with the ACO ones. However, we also see that the best solutions obtained with it are near-optimal and competitive with the ACO methods. Notice that this fact is relevant as it clearly shows the effectiveness of Algorithm 7 which is mainly based on its learning capability and simplicity. Concerning the gaps reported in Table 13, we observe that, for some random and Euclidean instances, Algorithm 7 allows obtaining tighter gaps than ACOQ in less CPU time. Finally, we observe that the number of iterations required by Algorithm 7 is significantly larger than the ACO methods. This can be explained by the fact that each iteration of Algorithm 7 requires a considerable less computational effort.From Table14, we observe that the distance between the initial and best solutions obtained with Algorithm 7 is even larger than for the instances presented in Table 13. This confirms again that the learning capability of Algorithm 7 is effective. Next, we see that the CPU time required by Algorithm 7 is significantly lower than the amount required by ACOQ for random instances. We also see that the gaps reported in Table 14 are tighter than those reported in Table 12 for the ACOQ approach for all random instances. In particular, we see that, for the random instances #2 and #3 in Table 14, the gaps obtained by Algorithm 7 are smaller than those reported for both ACO methods in Table 12. Regarding the Euclidean instances, we observe that the gaps obtained are still competitive when compared to those reported for the ACO methods. In fact, we also obtain negative gaps which means we outperform the best solutions obtained with the MILP models in less than 1 h. Finally, from Table 14, we observe that the gap obtained for the Euclidean instance #4 outperforms both gaps reported for the ACO methods in Table 12.


### 5.5. Average Numerical Results for the Proposed Algorithms

In order to give more insights with respect to the behavior of the proposed algorithms, in Figures2 and 3, we report average upper bounds and CPU times in seconds obtained with both VNSR and VNSQ for different values of k. More precisely, we randomly generate 20 large-size instances using n=300 nodes for each value of k in each figure. In particular, in Figure 2, these averages are reported for instances using random costs, whereas in Figure 3, these average values are reported for Euclidean instances.Figure 2
Average upper bounds and CPU times in seconds obtained withVNSR and VNSQ while varying k for random complete graph instances of size n=300 using d=3.Figure 3
Average upper bounds and CPU times in seconds obtained withVNSR and VNSQ while varying k for Euclidean complete graph instances of size n=300 using d=3.From Figures2 and 3, we mainly observe that the average upper bounds obtained with VNSQ are significantly smaller than those obtained with VNSR. This fact clearly shows that the embedded Q-learning strategy of VNSQ outperforms the classical near-far local search approach. Next, we further notice that the average CPU times are significantly smaller for VNSQ than for VNSR when solving the Euclidean instances. On the opposite, the VNSQ approach requires a significantly higher CPU time effort when solving instances with random costs. Similarly, in Figures 4 and 5, we report average upper bounds and CPU times in seconds for ACOR, ACOQ, and QL algorithms while varying k. For this purpose, again, we randomly generate 20 large-size instances of n=300 nodes for each value of k. In particular, in Figure 4, these averages are reported for instances with random costs. However, in Figure 5, these values are reported for Euclidean instances.Figure 4
Average upper bounds and CPU times in seconds obtained withACOR, ACOQ, and QL while varying k for random complete graph instances of size n=300 using d=2.Figure 5
Average upper bounds and CPU times in seconds obtained withACOR, ACOQ, and QL while varying k for Euclidean complete graph instances of size n=300 using d=2.From Figures4 and 5, we mainly observe that the average upper bounds obtained with ACOQ are slightly larger than those obtained by QL. In turn, the average upper bounds obtained with QL are larger than those obtained with ACOR. We further notice that this trend is more evident when solving graph instances with Euclidean distance costs. On the opposite, these bounds seem to be closer for the instances with random costs. Finally, we observe that the average CPU time values are significantly smaller for ACOR than for QL and ACOQ. Similarly, QL requires less CPU time than ACOQ. In conclusion, we observe that ACOR outperforms both QL and ACOQ. However, the QL method outperforms the ACOQ approach.


## 5.1. Numerical Results for the MILP Models

In Tables1–4, we present numerical results for P1, P1h, P2, and P2h. Notice that all these tables present the same column information. More precisely, column 1 shows the instance number, whilst columns 2 and 3 present the number of nodes of each graph instance and the value of k, respectively. Next, in columns 4–8 and 9–13, we present the optimal or best solution obtained with each model in at most 2 hours of CPU time, number of branch and bound nodes, CPU time in seconds, the optimal solution of each LP relaxation, and its CPU time in seconds. Finally, columns 14 and 15 present gaps that we compute by Opt−LP/Opt∗100 where Opt and LP refer to the optimal solution found with the MILP and LP models, respectively. Also, notice that each row in Tables 1 and 3 and in Tables 2 and 4 corresponds to a same instance.Table 1
Numerical results obtained withP1 and P1h for random graphs using d=3 and d=2.
#nkP1P1hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap1 %Gap1h%d=31501046.86750.8043.380.1746.861091.0643.380.197.437.43220113.974980.75109.460.20113.973620.73109.460.193.963.96330204.3514702.79199.050.19204.3519113.48199.050.172.592.59450459.0910051.76453.490.19459.0929015.74453.490.171.221.22510025115.44231.95114.630.48115.44251.79114.630.420.690.69650263.9525018.30262.040.47263.9518726.43262.040.420.730.73775454.4431865.04451.580.42454.44171524.23451.580.420.630.638100749.7536160202.21745.090.45749.7565193338.15745.090.440.620.62920050162.55103510.39160.881.73162.55146210.86160.881.641.031.0310100378.4380814506.06374.871.67378.4394564518.29374.871.650.940.9411150652.8629440936417.84648.891.73652.8626357647200648.891.950.610.61122001049.2540385372001044.742.261049.2521541572001044.741.780.430.43d=21501052.22491.5846.300.1752.221290.6146.300.1711.3411.34220134.082372.06129.710.22134.089042.57129.710.223.263.26330240.30451.26239.320.20240.3091.26239.320.200.410.41450564.3729487.38557.080.19564.3719537.47557.080.171.291.29510025123.714044.63121.540.48123.711483.99121.540.451.751.75650287.78109337.88285.310.61287.7818757.75285.310.500.860.86775529.5601.92529.310.50529.5613.99529.310.480.050.058100915.296839.75914.630.48915.29778.14914.630.450.070.07920050184.29012.28184.292.11184.2907.02184.292.030010100437.08642499.08435.452.50437.08256893.76435.452.280.370.3711150775.8927594.08774.762.40775.899042.20774.762.200.150.15122001285.3219067.751284.582.231285.3245998.281284.581.890.060.06Table 2
Numerical results obtained withP1 and P1h for Euclidean graphs using d=3 and d=2.
#nkP1P1hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap1 %Gap1h %d=31Berlin5210274.46102.25253.210.34274.46600.94253.210.287.747.74220926.77543991218.23801.280.20926.77591996300.65801.280.1913.5413.543301935.71424693336001666.190.201935.71529283236001666.190.1913.9213.924526081.63202461636005567.130.206081.63168860136005567.130.198.468.465gr962574.14488206360040.180.4174.14452709360040.180.3945.8145.81650163.115132933600117.660.55162.804629023600117.660.4227.8627.73775269.155935063600222.470.48269.236727753600222.470.4117.3417.37896436.6511923673600384.510.80436.6911687843600384.510.7011.9411.959ch150501607.2424334236001258.470.951610.9315735336001258.470.9421.7021.8810802657.2118647436002303.051.012694.8815476836002303.050.9013.3314.54111003457.8619991436003058.710.973441.5912931836003058.710.9211.5411.13121505914.1521578736005522.970.955923.7823488736005522.970.896.616.77d=21Berlin5210279.2300.52261.790.30279.2300.53261.790.226.256.252201015.3922537.08862.040.251015.39873515.91862.040.2215.1015.103302065.12944612.211871.760.282065.125439.361871.760.279.369.364526968.7716998.246545.080.276968.7712447.446545.080.236.086.085gr962574.59380670360045.190.7774.59324849360045.192.3939.4239.42650168.532704153600128.730.92174.251947523600128.738.4523.6226.12775283.141548383600248.4413.60288.493104293600248.440.5812.2513.88896474.951685393600452.149.66476.083975403600452.140.724.805.039ch150501753.839874236001356.911.041807.838631036001356.919.9122.6324.9410802917.727382936002489.252.292881.726021036002489.2510.6414.6913.62111003744.476671036003302.5415.884013.585451836003302.5414.0511.8017.72121506411.736274336006158.671.156397.127057436006158.671.013.953.73Table 3
Numerical results obtained withP2 and P2h for random graphs using d=3 and d=2.
#nkP2P2hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap2 %Gap2h %d=31501046.86422.2943.380.2346.863401.2543.380.207.437.43220113.97216116.65109.460.25113.9799911.42109.460.273.963.96330204.351161.67197.610.19204.35202414.65197.610.193.303.30450459.091231.00444.770.25459.096293.87447.870.203.122.44510025115.4435846.71114.540.53115.449754.60114.540.700.780.78650263.95124373.63261.020.51263.9566877.50261.020.751.111.11775454.44459143.16449.510.55454.44525111.53449.510.671.081.088100749.752422.56728.702.29749.75040.36734.840.592.811.99920050162.554793995.62160.362.21162.559202835.21160.363.291.351.3510100378.434792671.16373.283.81——7200373.283.681.36—11150652.864792097.25642.844.71652.864651663.29642.843.491.531.53122001049.25607505.271030.0430.111049.255021490.461034.233.211.831.43d=21501052.2251.0146.420.2552.221901.7346.300.2011.1211.34220134.083793.43129.900.33134.08151521.95129.710.203.123.26330240.302357.35239.320.28240.3073619.08239.320.200.410.41450564.3748820.20557.480.28564.37164742.00557.080.191.221.29510025123.711067248.12121.551.50123.711933237.10121.540.701.741.75650287.781676632.83285.352.84287.78857356.21285.310.810.840.86775529.5608.88529.313.23529.561050639.40529.310.760.050.058100——7200914.604.13——7200914.630.89——920050184.2910228.02184.2913.13——7200184.293.810—10100——7200435.4526.49——7200435.456.61——11150——7200774.7655.04——7200774.766.62——122001285.3210451408.711284.5894.68——72001284.584.150.06——: no solution found.Table 4
Numerical results obtained withP2 and P2h for Euclidean graphs using d=3 and d=2.
#nkP2P2hGapsOptB&BnTimeLPTimeOptB&BnTimeLPTimeGap2 %Gap2h %d=31Berlin5210274.4600.97253.210.21274.4600.82253.210.217.747.74220926.7702.25697.510.25926.7702.15697.510.2124.7424.743301935.711752.311414.090.221935.714276.121414.090.2026.9526.954526081.635562.004822.430.286081.63952.204923.770.1920.7119.045gr962572.402518268.8738.470.5172.40913301.4238.470.6146.8746.87650162.803274316.75110.480.56162.801816250.78110.460.6032.1432.15775267.57591109.28210.191.86267.57615157.00210.180.6821.4421.45896436.2348070.35354.603.72436.2398542.30357.830.6018.7117.979ch150501567.771734921.271219.691.21——36001219.691.4622.20—1080——36002237.221.25——36002237.221.85——11100——36002974.081.603399.1710371009.762974.081.97—12.51121505882.45982212.275339.1827.695882.45974245.395353.541.619.248.99d=21Berlin5210279.2302.95262.400.25279.2300.67261.790.206.036.252201015.39148721.75868.220.371015.39843.21862.040.1914.4915.103302065.1204.551874.560.342065.1202.871871.760.229.239.364526968.7797231.656554.140.286968.7702.896545.080.205.956.085gr962573.561592845.8245.410.7773.562525621.6445.190.6238.2738.57650168.53862299.41129.002.81168.53601350.55128.730.6723.4623.62775283.141235408.99248.894.12283.141672447.78248.440.6412.1012.25896474.952874329.41456.192.50474.952482315.00452.140.563.954.809ch15050——36001358.7316.52——36001356.912.46——10803167.66124136002491.5222.79——36002489.252.3421.35—111005425.41122436003306.7420.47——36003302.542.7939.05—121507176.67369836006166.4929.53——36006158.671.7814.08——: no solution found.We limit CPLEX to 2 hours of CPU time for the random input graphs, whilst for the Euclidean ones, we limit CPLEX to 1 hour in order to avoid CPLEX shortage of memory events. Consequently, each reported solution is optimal if its CPU time is less than 2 hours. Otherwise, it corresponds to the best solution obtained with CPLEX in at most 2 hours. Subsequently, in Tables5 and 6, we present numerical results for both P3 and P3h. These two tables also present the same column information. In particular, columns 1–3 present the same information as in Tables 1–4. Next, in columns 4–8 and 9–13, we report the optimal or best solution obtained with Algorithm 1 in at most 1 hour of CPU time, number of branch and bound nodes, CPU time in seconds, number of cycles added to each model, and number of iterations as well. The number of branch and bound nodes, in this case, corresponds to the sum of all branched nodes within each iteration of Algorithm 1. Notice that Theorem 2 in [24] ensures that the solutions obtained within each iteration of Algorithm 1 are lower bounds for the optimal solution of the problem. Consequently, an optimal solution is obtained when the corresponding CPU time is less than 1 hour. Notice that obtaining tight lower bounds is of crucial importance when developing exact methods as it allows to speed up the process significantly. Further notice that models P1, P1h, P2, and P2h obtain upper bounds for the problem if the optimal solution cannot be reached within 2 hours. Thus, we provide an interval where the optimal solution lies.Table 5
Lower bounds obtained withP3 and P3h for random graphs using d=3 and d=2.
#nkP3P3hOptB&BnTime#Cycles#IterOptB&BnTime#Cycles#Iterd=31501046.86531.753446.86362.1534220113.9701.3323113.9721.0123330204.3510215.621414204.357555.631414450459.099556.402322459.096115.992322510025115.44534.0123115.44222.9023650263.95337324.621920263.95339426.101920775454.44715458.034344454.44755859.7543448100749.75214640441.97211211749.75157118446.80211211920050162.55823.4545162.551523.784510100378.2414387063604.81252252378.2012412773624.1123723711150651.4915097493626.19263263651.4312140613614.95239239122001046.968528083602.983353341046.9510365043615.80320319d=21501052.2260.971252.2200.8712220134.08111.6144134.08171.5644330240.30261.3112240.30101.0812450564.3700.5232564.3700.4732510025123.711663.5612123.71362.9512650287.784997.1965287.782226.0465775529.5600.9001529.5600.84018100915.2904.0745915.2903.3445920050184.2904.3201184.2904.040110100437.08327431.7344437.08107530.844411150775.89018.5312775.89019.2012122001285.3207.05121285.3206.4312Table 6
Lower bounds obtained withP3 and P3h for Euclidean graphs using d=3 and d=2.
#nkP3P3hOptB&BnTime#Cycles#IterOptB&BnTime#Cycles#Iterd=31Berlin5210274.4601.4444274.4621.1244220926.772858450.958180926.774293162.0981803301935.71432870650.352982941935.71384635710.192982944526081.6315228434720.069379326076.5712618693611.008608545gr962568.1427826393619.001699268.0725079133621.4816184650154.0618782433716.3219665153.8917406253616.4219062775265.1111887453614.37297152264.8611130793626.85274129896427.944341603602.41734638427.943802493606.057336379ch150501545.7514873463642.21172791545.2313079163697.271677510802599.8014517303659.13212772598.4013279313776.4720272111003388.098590903682.262261083387.087975713617.0121698121505882.451334622909.894994195882.451270322874.88499419d=21Berlin5210279.2300.5122279.2300.67222201015.39642.471551015.39622.091553302065.121342.231542065.12632.421544526968.7700.731036968.7700.751035gr962573.566891771525.831255173.567645551832.4012551650168.53189277981.3210936168.533036961584.5610936775283.1424803100.756318283.141588096.286318896474.951195123.933410474.95480217.5834109ch150501637.179042113620.84112351636.059899404001.451063110802788.875649183790.60127322786.685396993697.0312330111003627.9328619279.8866143627.9336800341.066614121506368.79815243.7138106368.79904645.693810From Table1, first we observe that both P1 and P1h allow obtaining the optimal solution of the problem for most of the instances in less than 2 hours. In particular, for the degree value of d=3, we see that P1 and P1h cannot solve the instances #12 and #11-#12, respectively, whilst for d=2, all the instances are solved to optimality. Regarding the CPU times, for d=3, we see that P1h requires a higher amount of CPU time, whilst for d=2, we obtain similar values with both models. Next, we observe that the number of branch and bound nodes is slightly lower for P1 than for P1h. We also see that the LP bounds are near-optimal and exactly the same for both models which are confirmed by the gaps. Notice that all the LP relaxations are solved in less than 3 seconds. Finally, we observe that the CPU time values required to solve both models increase with k.Regarding the Euclidean graph instances presented in Table2, we observe that most of them cannot be solved to optimality in one hour. We mention that these instances proved to be very difficult to solve, and this is the main reason we limit CPLEX to one hour in order to avoid CPLEX shortages of memory. The difficulty in solving these instances is also reflected in the number of branch and bound nodes and in the gaps obtained which are significantly higher compared to Table 1. Notice that these gaps decrease with k in which evidences solving the DCkMST problem is harder to solve than the classical DCMST problem, at least for these instances. Next, we observe that the CPU time values of the LP relaxations are higher than those reported in Table 1 and that the degree values do not seem to affect the performance of each model. We further see that the optimal or best objective function values obtained for the instance “ch150” are significantly higher for d=2. Finally, we observe that the objective function values obtained with P1 for the instances #7–#10, #12 and #6–#9, and #11 are lower than those obtained with P1h for d=3 and d=2, respectively. However, the opposite occurs for the instances #6, #11 and #10, #12 while using d=3 and d=2, respectively.From Table3, we mainly observe that, for d=3, only P2 allows obtaining the optimal solution of the problem for all the instances in less than 2 hours. Notice that P2h only failed to find a feasible solution for the instance #10, whilst for d=2, neither P2 nor P2h can solve all the instances to optimality. However, P2 can still solve more instances than P2h. Regarding the number of branch and bound nodes and LP times, we observe smaller and slightly higher values than in Table 1, respectively. Finally, we observe that the LP bounds are near-optimal which is confirmed by the gaps. In Table 4, we observe that P2 and P2h can solve more instances to optimality than P1 and P1h in Table 2, respectively. However, in this case, both P2 and P2h cannot find a feasible solution for some of the instances. Next, we observe smaller values for the branch and bound nodes and CPU times with similar orders of magnitude for the LP relaxations when compared to Table 2. Finally, we observe that the LP bounds are not tight which is again confirmed by the gap columns.From Table5, we mainly observe that models P3 and P3h, which are both solved with Algorithm 1, allow to solve to optimality the instances #1–#9 and #1–#12 for d=3 and d=2, respectively. In particular, we see that, for d=3, it is harder to solve these instances than for d=2. This fact is reflected in the number of branch and bound nodes, CPU times of the LP relaxations, number of iterations, and number of cycles added to each exponential model which are clearly lower when d=2. Notice that, for the instances #10–#12 and degree value of d=3, we only report the best objective function values obtained which are in fact lower bounds. Regarding the Euclidean instances reported in Table 6, we observe that the instances #1–#3, #12 and #1–#8, #11-#12 are all solved to optimality for d=3 and d=2, respectively. We also see in Table 6 that it is harder to solve the instances for d=3 than for d=2. Notice that, for d=2, we almost solve all the instances to optimality with the exception of instances #9 and #10 for which we obtain lower bounds. Again, this observation can be verified by looking at the number of branch and bound nodes, CPU times of the LP relaxations, number of iterations, and number of cycles added to each proposed model which are smaller when d=2. Finally, from the numerical results presented in Tables 1–6, we can conclude that model P1 outperforms the other ones as it allows to obtain either an optimal or a best upper bound for most of the tested instances. Notice that the flow and exponential models also have good performance in terms of optimality, but for the large-size instances, the flow models deteriorate rapidly, whilst the exponential ones cannot be handled efficiently in terms of CPU times. Consequently, in Tables 7 and 8, we present numerical results for P1 for large random and Euclidean input graph instances with up to 400 and 280 nodes, respectively. In particular, in Table 7, these numerical results are reported for d=3, whereas in Table 8, these numerical results are obtained while using a degree value of d=2.Table 7
Numerical results obtained withP1 for large input graphs using d=3.
#nkP1OptB&BnTimeLPTimeGap1%Random graph instances130050118.8212040.30117.394.401.202100283.63371964.73281.704.200.683200703.429949132022.93701.154.210.3243001334.6117753872001332.445.010.16540050110.201651193.66108.488.471.566100236.59100882316.85234.519.480.887200522.6113577259.07520.668.240.3784001459.533176372001456.3410.840.22Euclidean graph instances1a28050480.35232223600384.004.2120.0621001011.07429573600784.004.5222.4631501316.172876136001184.164.5710.0342802493.684623236002385.763.854.33Table 8
Numerical results obtained withP1 for large input graphs using d=2.
#nkP1OptB&BnTimeLPTimeGap1%Random graph instances130050132.92034.03132.876.910.042100321.64454309.74320.828.300.263200820.70272520.03820.0212.600.0843001623.8512501301.401623.1611.510.04540050123.024132003.07120.0519.052.426100265.68305873073.49262.4419.191.227200598.2025272448.10595.0834.010.5284001818.0399972001797.7616.211.12Euclidean graph instances1a28050534.97136853600384.005.4828.2221001133.49232123600784.004.9330.8331501417.201140536001186.716.1616.2642002991.961103636002513.314.6616.00From Tables7 and 8, we observe similar trends as for the above Tables 1–4. More precisely, we observe that P1 can solve almost all random input graph instances to optimality in less than 2 hours using both degree values. Only the instance number #8 could not be solved to optimality although an upper bound is reported for this particular instance. In contrast, none of the Euclidean instances can be solved to optimality in less than 1 hour of CPU time. For these instances, only upper bounds are reported too. The difficulty in solving the Euclidean instances can also be observed in the gap columns which report significantly higher values as a consequence of the LP bounds obtained. On the opposite, the LP bounds obtained for the random graph instances are near-optimal. Finally, we observe from both Tables 7 and 8 that the number of branch and bound nodes is significantly smaller for d=2 than for d=3 and that the LP times are slightly higher for d=2.


## 5.2. Numerical Results for VNS Algorithms

In Tables9 and 10, we present numerical results obtained with the proposed VNS algorithms for both random and Euclidean input graph instances while using a degree value of d=3. Recall that our second VNS approach consists of replacing the code lines of Algorithm 3 corresponding to the random local search strategy with the code of Algorithm 4 which represents the embedded Q-learning strategy. Hereafter, we denote by VNSR and VNSQ these two VNS approaches. In Algorithms 3 and 4, we arbitrarily set the maximum CPU time allowed to a value of maxTime=250 seconds, whereas the input parameters indKMax and indKM are set to a value of 10. Similarly, the required parameters of Algorithm 4 were calibrated to λ=0.25, μ=0.25, and Prob=0.15. Finally, the reward value and parameter θ of Algorithm 2 were set to Reward=1 and θ=0.1, respectively.Table 9
Numerical results obtained with VNS algorithms for random and Euclidean instances usingd=3.
#nkMILPVNSRVNSQBestTime (s)VNSiniVNSTime (s)#IterGapRVNSiniVNSTime (s)#IterGapQInput graphs with random uniform costs1501046.860.80254.7648.276.06138423.00208.7346.860.6216200220113.970.73266.68118.5322.8072134.00353.48116.078.2927961.84330204.351.67338.03206.3546.5136740.97326.75205.4330.4925030.52450459.091.00463.14463.140.0810.88463.14463.140.0810.88510025115.441.79352.14115.44115.10238260413.40118.7718.5988542.88650263.956.43564.30269.1861.8929621.98517.86265.4286.4173820.55775454.445.04600.65465.48197.6126812.42659.23456.59215.4640980.478100749.7522.56749.87749.870.2510.01749.87749.870.2410.01920050162.5510.39555.42179.981857.836364010.72463.99169.38477.84187284.2010100378.43506.06725.35395.982928.53123244.63763.28386.121657.5291942.0311150652.861663.29883.28721.071250.9191210.44926.60670.361130.828722.68122001049.25505.271056.341056.343.9810.671056.341056.343.9010.67Input graphs with Euclidean distance costs1Berlin5210274.460.821976.48274.461.41595901645.70274.460.094700220926.772.153923.04926.772.69738203161.01926.770.3682003301935.712.314961.302104.325.2951338.714484.232071.97178.341930487.034526081.6326081.636081.630.01106081.636081.630.00105gr962572.40268.87204.0079.28122.542303369.50247.0774.14271.895085832.40650162.80250.78307.92179.5074.783409410.25260.74163.702.7411630.55775267.57109.28388.13269.8820.9730490.86372.46267.57141.16211760896436.2342.30436.23436.230.0210436.23436.230.01109ch150501567.77921.273441.771716.09348.671547349.463290.911696.0637.20161078.1810802657.2136004327.972792.69115.65138895.094480.372771.08223.80216754.28111003399.171009.764722.643527.70408.36153663.784745.733499.01288.25138582.93121505882.45212.275882.455882.450.12105882.455882.450.1210Table 10
Numerical results obtained with VNS algorithms for the instances of Table7 using d=3.
#MILPVNSRVNSQBestTime (s)VNSiniVNSTime (s)#IterGapRVNSiniVNSTime (s)#IterGapQInput graphs with random uniform costs1118.8240.30534.32139.361330.834742917.28510.03127.26494.26362917.102283.6364.73762.97338.442549.43766519.32775.47298.831893.4284815.353703.422022.931088.49850.681381.5532620.931000.36750.053851.549666.6241334.6172001349.381349.3832.6811.101349.381349.3832.7011.105110.20193.66492.26126.271376.585215614.58526.08124.871015.616372013.316236.592316.85783.62306.031465.56824229.35677.53257.463628.70202968.827522.61259.071129.36696.202587.2766433.211058.37635.214130.04104521.5481459.5372001485.591485.5983.4511.781485.591485.5983.1311.78Input graphs with Euclidean distance costs1480.3536001095.00446.91314.4196332−6.961076.54440.3874.3931647−8.3221011.0736001385.72897.55425.8011551−11.221447.17890.15195.388251−11.9531316.1736001774.881304.78786.148339−0.861689.481303.48570.315417−0.9642493.6836002440.942440.947.941−2.112440.942440.947.941−2.11In particular, in Table9, we report numerical results for the instances presented in Tables 1–6, whereas in Table 10, we report numerical results for the large-size instances presented in Table 7. In both tables, columns 4–15 and 2–13 contain exactly the same column information, respectively. More precisely, in Table 10, columns 1 to 3 present the instance number, number of nodes, and the value of k. Next, in columns 4 and 5, we report the minimum objective function and CPU time values reported for each instance in Tables 1–4. Subsequently, in columns 6–10 and 11–15, we report for each VNS approach, the initial and minimum objective function values, CPU time in seconds required by VNS, number of iterations, and gaps, respectively. The gaps are computed by VNS−Best/Best∗100.From both Tables9 and 10, first we observe that the objective function values of the initial solutions are significantly higher than the best ones. The latter clearly evidences the effectiveness of both VNS approaches. Notice that when k=n, the problem reduces to the classical DCMST problem. Consequently, for these instances, the solutions are obtained only with Algorithm 2. Next, we observe that the CPU times required by both VNS methods are larger for the random instances than for the Euclidean ones. From Table 9, next, we observe that the solutions obtained with VNSQ outperform those obtained with VNSR for both random and Euclidean instances. Although, in general, we see that both VNS procedures allow obtaining near-optimal solutions which is confirmed by the gap columns. On the opposite, we see that the objective function values reported in Table 10 are not tight for the random instances when compared to those obtained by the MILP models. But still, in this case, the gaps reported for VNSQ are significantly better than for VNSR. However, for the Euclidean instances reported in Table 10, we see that both VNS approaches allow us to obtain better solutions than the MILP models. The latter can be verified by the negative gaps which show that the objective function values obtained with VNS algorithms are significantly lower than those obtained with the MILP models.


## 5.3. Numerical Results for ACO Algorithms

Now, we report numerical results obtained with Algorithms5 and 6 for random and Euclidean instances while using a degree value of d=2. More precisely, in Table 11, we report numerical results for the instances presented in Tables 1–6. However, in Table 12, we report numerical results for the large-size instances presented in Table 8. Hereafter, we denote by ACOR and ACOQ the two ant colony optimization approaches presented in Algorithms 5 and 6, respectively. Recall that the first one is constructed based on the classical ACO metaheuristic [26]. However, the latter is based on the Ant-Q algorithm proposed in [25, 27]. In Algorithm 5, we calibrated the input parameters to α=0.5, β=2, ρ=0.9, and q=1. However, in Algorithm 6, we calibrated the input parameters to α=1, β=4, q=1, Prob=0.1, λ=0.15, and μ=0.15. Finally, in both Algorithms 5 and 6, we arbitrarily set the maximum CPU time and number of ants to maxTime=250 seconds and Ants=20, respectively.Table 11
Numerical results obtained with ACO algorithms for random and Euclidean instances usingd=2.
#nkMILPACORACOQBestTime (s)ACOiniACOTime (s)#IterGapRACOiniACOTime (s)#IterGapQInput graphs with random uniform costs1501052.220.6193.6952.2254.7855760259.2652.2289.8920830220134.081.61175.53135.7524.5825511.24955.85135.75126.1721111.24330240.301.08370.28244.3126.0126481.66435.68245.24234.2544372.05450564.370.47684.50584.185.535613.512486.65588.3148.698034.24510025123.712.95163.55123.712.066701231.26123.71192.919950650287.786.04337.69293.9029.209752.122546.91298.88339.6617313.85775529.560.84650.54550.4146.0015333.933804.55563.89220.2811646.488100915.293.341093.20970.4350.6616996.025188.971021.63239.28123011.61920050184.294.04279.89193.34168.4716264.912495.55198.36698.7110937.6310100437.0830.84546.76467.7668.136577.014856.80492.90258.6340812.7711150775.8918.53941.43836.899.54907.867444.72870.56273.5346012.20122001285.326.431734.021480.7615.9115615.2010277.041614.61321.3052925.61Input graphs with Euclidean distance costs1Berlin5210279.230.512195.31279.23389.574209604581.09279.2321.4435102201015.392.091270.761015.392.0120109757.741020.0272.2411890.453302065.122.232757.232065.1258.286258012888.362069.9032.445350.234526968.770.739738.247156.6517.7118952.6922432.487234.276.471063.805gr962573.56621.64112.7874.59144.6756471.40144.7274.807.65451.68650168.53299.41219.39172.3897.1838112.28546.20180.29161.809846.97775283.1496.28379.18292.2024.979663.19609.61302.6650.223026.89896474.9517.58619.16497.0666.0225784.65944.82534.30320.22195612.499ch150501753.8336002050.391678.9338.58672−4.2718522.681722.83222.95613−1.7610802881.7236003667.842843.87169.472972−1.3127659.442929.44141.143941.65111003627.93279.884833.863788.2698.1817184.4135594.853773.76325.019084.01121506368.7943.717566.166720.3142.557425.5152740.447005.90113.5031510.00Table 12
Numerical results obtained with ACO algorithms for random and Euclidean instances of Table8 using d=2.
#MILPACORACOQBestTime (s)ACOiniACOTime (s)#IterGapRACOiniACOTime (s)#IterGapQInput graphs with random uniform costs1132.9234.03193.95146.5326.599810.232228.39158.671443.44109119.372321.64309.74419.46367.1643.3016514.154620.75396.94589.5246823.413820.70520.03980.56898.2564.562379.449596.31923.88918.3670312.5741623.851301.402122.151933.14128.8249519.0414537.582107.741130.5586529.795123.022003.07172.02136.3398.3021010.812581.64143.352080.9391716.526265.683073.49328.39298.32155.4733812.285243.13316.621187.5353019.177598.202448.10754.03683.366.101314.2310356.11724.671247.0556321.1481818.0372002370.222224.405.611122.3520229.562631.201151.9052144.72Input graphs with Euclidean distance costs1534.973600426.30401.87113.06632−24.87675.36405.20227.45196−24.2521133.493600990.04846.6796.56539−25.302944.45919.303.483−18.8931417.2036001554.941311.6242.31236−7.442909.141358.90312.40270−4.1142991.9636003477.832983.922.1811−0.266137.083239.96292.152558.28The legend of Table11 is as follows. From columns 1 to 3, we present the instance number, the number of nodes of the input graph, and the value of k. Next, in columns 4 and 5, we present the minimum objective function and CPU time values reported in Tables 1–6. We repeat this information for the sake of clarity. Next, in columns 6–10 and 11–15, we report the objective function values of the initial solutions, the objective function values of the best solutions, CPU time in seconds, number of iterations, and gaps which are computed by ACO−Best/Best∗100, respectively. Similarly, the first three columns of Table 12 report the instance number and the best objective function values and CPU times reported in Table 8. Again, this information is repeated for comparison purposes. Finally, the legends of columns 4–8 and 9–13 are exactly the same as for the columns 6–10 and 11–15 in Table 11.From Table11, we observe that the best solutions obtained with both ACO approaches are near-optimal and far from the initial solutions obtained. These facts prove the effectiveness of the proposed methods. Next, we see that the CPU times are lower for ACOR than for ACOQ for most of the instances and in particular for the larger ones. Regarding the number of iterations, in general, we observe similar orders of magnitude for both methods. Finally, we observe that the gaps obtained with ACOR are tighter than those obtained with ACOQ for both the random and Euclidean instances. In particular, for large instances, these values are significantly better. Finally, notice that we obtain negative gaps for some of the large Euclidean instances. This clearly evidences that the solutions obtained with the ACO methods outperform significantly the solutions obtained with the MILP models.From Table12, we mainly observe that the solutions obtained with both ACO algorithms are significantly worse than those reported in Table 11 for random instances. On the opposite, for all the Euclidean instances, the solutions obtained with the ACO approaches are significantly better than those obtained with P1 which is again confirmed by the negative gaps obtained.


## 5.4. Numerical Results for QL Algorithm

In Tables13 and 14, we present numerical results obtained with Algorithm 7 (denoted as QL) for both random and Euclidean instances while using a degree value of d=2. In Tables 13 and 14, we report numerical results for the same instances presented in Tables 1–6 and for the instances of Table 8, respectively. The legend of Table 13 is as follows: columns 1–3 present the instance number, number of nodes (or name of the Euclidean instance), and the value of k, respectively. Notice that the name of each Euclidean instance indicates at the end the number of nodes it contains. For example, the instance name “Berlin52” has 52 nodes. Next, columns 4 and 5 report the best solution and minimum CPU time in seconds obtained with the MILP models. Finally, columns 6–10 report the initial solution obtained with Algorithm 7, its best solution found, CPU time in seconds, number of iterations, and gaps obtained which are computed by QL−Best/Best∗100, respectively. The legend of Table 14 is exactly the same as for Table 13 and is obtained by removing columns 2 and 3 from Table 13.Table 13
Numerical results obtained with Q-learning algorithm for random and Euclidean instances usingd=2.
#nkMILPQLBestTime (s)QLiniQLTime (s)#IterGapInput graphs with random uniform costs1501052.220.61427.2652.229.03423590220134.081.61794.00138.008.76406222.92330240.301.081473.09245.622.25100792.21450564.370.472567.98601.3715.38709226.55510025123.712.951207.44123.714.7598900650287.786.042440.56291.8998.122118911.42775529.560.843473.78565.78487.9010678996.838100915.293.344975.561024.7730.106434511.96920050184.294.042543.82199.27130.441226718.1210100437.0830.845027.12476.04123.871181058.9111150775.8918.537342.28843.6038.57365128.72122001285.326.439799.771546.62227.2121410720.32Input graphs with Euclidean distance costs1Berlin5210279.230.512692.71279.2339.9320032202201015.392.0912651.131020.02151.287743860.453302065.122.2316299.432070.6545.422283630.264526968.770.7330795.847234.2785.364339013.805gr962573.56621.64858.8074.5978.731994191.40650168.53299.411774.49180.0589.152251946.83775283.1496.282500.70303.35260.646483027.13896474.9517.583352.78552.06191.6847743016.239ch150501753.83360018911.281711.33125.94186763−2.4210802881.72360028695.922949.71202.883061202.35111003627.93279.8837580.423953.79107.611606248.98121506368.7943.7154468.427122.77565.4484402011.83Table 14
Numerical results obtained with Q-learning algorithm for the random and Euclidean instances of Table8 using d = 2.
#MILPQLBestTime (s)QLiniQLTime (s)#IterGapInput graphs with random uniform costs1132.9234.032455.40147.03116.827292410.612321.64309.745203.22347.51101.90630218.043820.70520.0310356.83897.52160.04998319.3641623.851301.4015547.511990.10441.3826392222.555123.022003.072173.05137.42155.636459711.706265.683073.495346.13307.22442.6018313015.637598.202448.1010535.56691.97286.7611571715.6781818.03720020088.172287.40305.8012663125.81Input graphs with Euclidean distance costs1534.9736006193.93426.46170.69119614−20.2821133.49360011619.78909.0382.0557482−19.8031417.20360018065.401441.78275.341924371.7342991.96360033966.222800.73130.9891758−6.39From Table13, first we observe that the initial solutions obtained with Algorithm 7 are considerably worse than those obtained with the ACO ones. However, we also see that the best solutions obtained with it are near-optimal and competitive with the ACO methods. Notice that this fact is relevant as it clearly shows the effectiveness of Algorithm 7 which is mainly based on its learning capability and simplicity. Concerning the gaps reported in Table 13, we observe that, for some random and Euclidean instances, Algorithm 7 allows obtaining tighter gaps than ACOQ in less CPU time. Finally, we observe that the number of iterations required by Algorithm 7 is significantly larger than the ACO methods. This can be explained by the fact that each iteration of Algorithm 7 requires a considerable less computational effort.From Table14, we observe that the distance between the initial and best solutions obtained with Algorithm 7 is even larger than for the instances presented in Table 13. This confirms again that the learning capability of Algorithm 7 is effective. Next, we see that the CPU time required by Algorithm 7 is significantly lower than the amount required by ACOQ for random instances. We also see that the gaps reported in Table 14 are tighter than those reported in Table 12 for the ACOQ approach for all random instances. In particular, we see that, for the random instances #2 and #3 in Table 14, the gaps obtained by Algorithm 7 are smaller than those reported for both ACO methods in Table 12. Regarding the Euclidean instances, we observe that the gaps obtained are still competitive when compared to those reported for the ACO methods. In fact, we also obtain negative gaps which means we outperform the best solutions obtained with the MILP models in less than 1 h. Finally, from Table 14, we observe that the gap obtained for the Euclidean instance #4 outperforms both gaps reported for the ACO methods in Table 12.


## 5.5. Average Numerical Results for the Proposed Algorithms

In order to give more insights with respect to the behavior of the proposed algorithms, in Figures2 and 3, we report average upper bounds and CPU times in seconds obtained with both VNSR and VNSQ for different values of k. More precisely, we randomly generate 20 large-size instances using n=300 nodes for each value of k in each figure. In particular, in Figure 2, these averages are reported for instances using random costs, whereas in Figure 3, these average values are reported for Euclidean instances.Figure 2
Average upper bounds and CPU times in seconds obtained withVNSR and VNSQ while varying k for random complete graph instances of size n=300 using d=3.Figure 3
Average upper bounds and CPU times in seconds obtained withVNSR and VNSQ while varying k for Euclidean complete graph instances of size n=300 using d=3.From Figures2 and 3, we mainly observe that the average upper bounds obtained with VNSQ are significantly smaller than those obtained with VNSR. This fact clearly shows that the embedded Q-learning strategy of VNSQ outperforms the classical near-far local search approach. Next, we further notice that the average CPU times are significantly smaller for VNSQ than for VNSR when solving the Euclidean instances. On the opposite, the VNSQ approach requires a significantly higher CPU time effort when solving instances with random costs. Similarly, in Figures 4 and 5, we report average upper bounds and CPU times in seconds for ACOR, ACOQ, and QL algorithms while varying k. For this purpose, again, we randomly generate 20 large-size instances of n=300 nodes for each value of k. In particular, in Figure 4, these averages are reported for instances with random costs. However, in Figure 5, these values are reported for Euclidean instances.Figure 4
Average upper bounds and CPU times in seconds obtained withACOR, ACOQ, and QL while varying k for random complete graph instances of size n=300 using d=2.Figure 5
Average upper bounds and CPU times in seconds obtained withACOR, ACOQ, and QL while varying k for Euclidean complete graph instances of size n=300 using d=2.From Figures4 and 5, we mainly observe that the average upper bounds obtained with ACOQ are slightly larger than those obtained by QL. In turn, the average upper bounds obtained with QL are larger than those obtained with ACOR. We further notice that this trend is more evident when solving graph instances with Euclidean distance costs. On the opposite, these bounds seem to be closer for the instances with random costs. Finally, we observe that the average CPU time values are significantly smaller for ACOR than for QL and ACOQ. Similarly, QL requires less CPU time than ACOQ. In conclusion, we observe that ACOR outperforms both QL and ACOQ. However, the QL method outperforms the ACOQ approach.


## 6. Conclusions

In this paper, we consider the degree-constrainedk-cardinality minimum spanning tree problem which emerges as a combination of two classical combinatorial optimization problems, namely, the degree-constrained and k-minimum spanning tree problems. One can see from the literature that this problem has not been studied in depth yet. This leads us to propose three mixed-integer linear programming models for which we derive equivalent formulations by using the handshaking lemma. In order to obtain near-optimal solutions, we further propose ant colony optimization, variable neighborhood search, and a pure Q-learning-based algorithm. In particular, for each proposed metaheuristic, we obtain new algorithms while embedding a Q-learning strategy. We conduct substantial numerical experiments using benchmark input graph instances from TSPLIB and randomly generated ones with random uniform and Euclidean distance costs with up to 400 nodes. From the numerical results obtained, our main conclusions can be listed as follows:(1)
We observe that, in general, the proposed models which are constructed based on Miller–Tucker–Zemlin-constrained approach are more robust than the flow ones as they allow to obtain optimal or best upper bounds for all the instances. In particular, we see that the flow models allow us to solve to optimality Euclidean instances with up to 100 nodes, which is not possible to achieve with the other models. However, the flow models cannot provide an upper bound for some of the instances with higher dimensions. Similarly, we observe that our proposed exponential models do also show good performance in terms of optimality, but again they cannot be handled efficiently when solving large-size instances. However, they provide an interval where the optimal solution lies. We further conclude that it is not evident to decide whether the performance of the proposed models improves or deteriorates while using the handshaking lemma. Finally, we observe that the Miller–Tucker–Zemlin-based models allow us to obtain optimal solutions for instances with up to 400 nodes while using random costs and degree values ofd∈2,3. On the opposite, they cannot solve all the Euclidean instances to optimality.(2)
Concerning the VNS algorithms, we observe that the objective function values obtained with both VNS algorithms are significantly lower when compared to their initial solutions obtained. This clearly evidences the effectiveness of our VNS approaches. In general, the two VNS procedures allow obtaining near-optimal solutions and even better solutions than CPLEX for the large-size instances. Next, we see that the CPU times required by these algorithms are larger for the random instances than for the Euclidean ones. Notice that Euclidean instances are significantly harder to solve by the MILP models. We can also conclude that the embedded Q-learning strategy in our VNS algorithm allows us to obtain better solutions than the classical near-far random local search strategy. This suggests that the construction of optimization methods with learning capabilities in order to make them robust, self-adaptive, and independent from decision-makers is worth to be further investigated.(3)
Regarding the ACO algorithms, we observe that both methods allow obtaining near-optimal solutions which certainly prove their effectiveness. In general, we obtain better solutions for the Euclidean instances than for the random ones. In particular, the quality of the Euclidean solutions obtained improves significantly for large-size instances of the problem for which both ACO methods obtain better solutions than CPLEX. Next, we observe that the solutions obtained withACOR are better in terms of quality when compared to the solutions obtained with ACOQ for both random and Euclidean instances of the problem. Concerning our pure Q-learning approach proposed to solve instances with degree d=2, we observe that the solutions obtained for random and Euclidean instances are competitive with those obtained with the ACO methods. In particular, we see that the CPU time required by this method is significantly lower than the amount required by ACOQ. Then, we further see that, for some of the instances, the gaps obtained by the pure Q-learning approach are tighter than those obtained with ACOQ. We also obtain better solutions than CPLEX solver which is used to solve the MILP model. Finally, notice that this pure Q-learning method is simple and versatile, and then, it can be adapted to any other combinatorial optimization problem in a straightforward manner.As future research, we plan to propose new formulations related to the degree-constrainedk-minimum spanning tree problem with applications on network design problems.



---

*Source: 7628105-2020-11-09.xml*