## 3 Heuristic Methods
### Dynamically Changing Weight of Penalties  
To move good solutions of mBSSRP from the solution of the mBSSRP-S, we must set appropriate values of the weight parameters \(\alpha\) and \(\beta\) in the objective function (Equation (1)) of the mBSSRP-S. If we set the parameters to small values, we cannot obtain feasible solutions owing to the small values of the parameters. Because the penalty becomes weak even when the solution has many constraint violations. Meanwhile, if we set the parameters to large values, the power of the parameters will be excessively strong, which is the same situation for solving mBSSRP. Thus, we dynamically change the parameters \(\alpha\) and \(\beta\) according to the solution state. To adjust the weight parameters \(\alpha\) and \(\beta\) dynamically, the objective function of the mBSSRP-S was modified as follows:  
\[\text{minimize}\hskip 56.905512 ptT(t)+\alpha(t)\times P_{T}(t)+\beta(t)\times( P_{B+}(t)+P_{B-}(t)) \tag{2}\]  
where \(\alpha(t)\) and \(\beta(t)\) represent weight parameters for the time limit, loading and supplying constraints at the \(t\)th iteration, \(T(t)\) denotes the total travel time at the \(t\)th iteration, \(P_{T}(t)\) denotes the excess working time beyond the time limit at the \(t\)th iteration, \(P_{B+}(t)\) denotes the loading violation of the bicycles at the \(t\)th iteration and \(P_{B-}(t)\) denotes the supplying violation of the bicycles at the \(t\)th iteration. The adjustment of the parameters is determined by comparing \(\alpha(t)\times P_{T}(t)\) and \(\beta(t)\times(P_{B+}(t)+P_{B-}(t))\) at the \(t\)th iteration. If the penalty of the time limit \(\alpha(t)\times P_{T}(t)\) is larger than the sum of penalties loading and supplying constraints \(\beta(t)\times(P_{B+}(t)+P_{B-}(t))\), \(\alpha(t)\) is increased and \(\beta(t)\) is decreased. If the penalty \(\beta(t)\times(P_{B+}(t)+P_{B-}(t))\) is larger than \(\alpha(t)\times P_{T}(t)\), \(\alpha(t)\) is decreased and \(\beta(t)\) is increased. If \(P_{T}(t)=0\) or \(P_{B+}(t)+P_{B-}(t)=0\), a corresponding weight \(\alpha(t)\) or \(\beta(t)\) is not updated. Then, the followings are used for adjusting the parameters:  
Figure 4: Example of the Or-opt. Black squares represent depots and gray circles represent ports. The red and blue lines represent each tour. (**a**) Tours before the Or-opt is executed. (**b**) A partial tour \(i-j\) is inserted into another tour \(k-l\) in the normal order (the orange rectangle). (**c**) The partial tour \(i-j\) is inserted into another tour \(k-l\) in the reverse order (the green rectangle). Figure 5: Example of the CROSS-exchange. Black squares represent depots and gray circles represent ports. Orange rectangles express normal order and green rectangles express a reverse order. Red and blue lines represent each tour. (**a**) Tours before the CROSS-exchange is executed. (**b**) Partial tours \(i-j\) and \(k-l\) are exchanged with the normal order. (**c**) Only the partial tour \(k-l\) is exchanged in the reverse order. (**d**) Only the partial tour \(i-j\) is exchanged in reverse order. (**e**) Partial tours \(i-j\) and \(k-l\) are exchanged in the reverse order. Figure 6: Example of the 2-opt method. Black squares represent depots and gray circles represent ports. (**a**) A tour before the 2-opt method is executed. (**b**) A tour after the 2-opt method is executed. if \(p_{\rm T}(t)=0\),  
\[\alpha(t+1) = \alpha(t), \tag{3}\]  
otherwise  
\[\alpha(t+1) = \left\{ \begin{array}{ll}\alpha(t)&\mbox{if }\,\,\alpha(t)p_{\rm T}(t)= \beta(t)(p_{\rm q+}(t)+p_{\rm q-}(t)),\\ \lambda\alpha(t),&\mbox{if }\,\,\alpha(t)p_{\rm T}(t)>\beta(t)(p_{\rm q+}(t)+p_{ \rm q-}(t)),\\ \mu\alpha(t),&\mbox{if }\,\,\alpha(t)p_{\rm T}(t)<\beta(t)(p_{\rm q+}(t)+p_{ \rm q-}(t)),\end{array} \right. \tag{4}\]  
and if \(p_{\rm q+}(t)+p_{\rm q-}(t)=0\),  
\[\beta(t+1)\,=\,\beta(t), \tag{5}\]  
otherwise  
\[\beta(t+1)\,=\,\left\{