## 1 Introduction  
Over the past five years, the global logistics industry has undergone a revolutionary digitalization transformation. According to Statista (2024) [1], the global parcel delivery quantity exceeded 150 billion in 2024, with an annual growth rate of 18.7%, and is projected to surge to 350 billion by 2030. In this trend, the optimization goals of logistics systems have shifted from traditional large-scale transportation to high-precision resource coordination [2]. The logistics industry now needs to improve the loading efficiency of individual delivery vehicles while optimizing delivery routes, as the synergy between these two factors directly impacts cost-saving potential. The parcels delivery problem is a common problem of logistics task in factories, hospitals, universities, communities, etc. Raw materials, medicines, and commodities need to be duly and correctly delivered to their destination by a delivery vehicle. Since the capacity and quantity of the delivery vehicles in a delivery center is limited, the delivery scheme should be properly planned to achieve a high delivery efficiency, especially when the parcels are numerous. Generally, as the shape, size, and destination of parcels are different, it is difficult to fully utilize the volume of the delivery vehicle while making sure the vehicle is driving along the shortest delivery route. Therefore, the capacitated vehicle routing problem with three-dimensional loading constraints (3L-CVRP) is an NP-hard problem [3]. The 3L-CVRP research dates back to 2006 when [PERSON] et al. [4] proposed a two-stage tabu search algorithm to introduce the rule of last-in-first-out (LIFO) and stability of parcels into route planning. The core value of 3L-CVRP research is achieving Pareto optimality in resource efficiency with limited space and time, and many researchers focused on this point from the aspects of three-dimensional (3D) loading and vehicle route planning. [PERSON] et al. [5] developed two improved relocation constraints to achieve the balance between transportation cost and operational complexity in 3L-CVRP, and the results demonstrate the significance of relocation allowances in reducing overall costs and improving the capacity utilization of the delivery vehicle. [PERSON] et al. [6] utilized a deep reinforcement learning (DRL) approach to solve the 3D bin packing problem, meaning the parcels loading problem can be solved by sequence, orientation, and position, independently. Similarly, [PERSON] et al. [7] proposed a DRL approach to solve the 3D packing problem, and presented a novel state representation of a packing environment. [PERSON] et al. [8] proposed an approach to program constraints in 3L-CVRP, and improved the results of 36 of 93 benchmark problems. [PERSON] et al. [9] developed a software toolbox to assist in solving the 3D bin packing problem in some specific industries. [PERSON] et al. [10,11] deemed that the routing problem is superior to the packing problem, and proposed an adaptive genetic algorithm with elastic strategy to solve the 3L-CVRP for many dynamic demands. The 3L-CVRP research is significant in the logistics industry, and many constraints can be considered, such as time windows [12; 13; 14], resource sharing [15], the stability of parcels [16], the heterogeneity of parcels [17], etc. These constraints affect the loading scheme and route planning. Therefore, to simplify the modeling and computing process, many researchers adopt a step-by-step optimization strategy of route planning followed by parcels loading [5] or parcels loading followed by route planning [18]. However, the traditional step-by-step optimization strategy exhibits low efficiency in both steps and low global suboptimality in high-density delivery scenarios (e.g., campus logistics), as the loading sequence and the route topology are nonlinearly correlated [19].