Container stowage method and system based on multi-stage optimization algorithm

Abstract

The application discloses a container loading method and system based on a multi-stage optimization algorithm and belongs to the technical field of container loading. The method comprises the following steps: obtaining cabin layout data of a target ship and physical attribute data, destination data and loading and unloading time window data of containers to be loaded; generating an initial loading scheme through a multi-stage optimization algorithm; in the first stage, performing mode matching and pre-grouping on the containers to be loaded based on a historical loading case library to obtain a candidate grouping set; in the second stage, performing cabin allocation on the candidate grouping set by combining ship stability constraints and strength constraints and adopting a tabu search algorithm to obtain a locally optimal loading scheme; and in the third stage, taking the locally optimal loading scheme as an initial population and performing multi-objective global optimization by using a genetic algorithm based on non-dominated sorting to output a final loading scheme. The method can efficiently solve a loading scheme that takes into account stability, strength and operation efficiency under complex constraints.

Description

Technical Field

[0001] This invention relates to the field of container loading and stowage technology, specifically to a container loading and stowage method and system based on a multi-stage optimization algorithm. Background Technology

[0002] With the trend towards larger container ships, the container capacity per vessel continues to increase, and the layout of cargo holds is becoming increasingly complex. When developing a stowage plan, it is necessary to comprehensively consider multiple attributes such as container weight, size, destination port, and loading / unloading time windows, while simultaneously meeting constraints related to ship stability, hull strength, and port operational efficiency. Existing stowage methods largely rely on manual experience or single heuristic rules, allocating cargo holds sequentially through a pre-set priority strategy. This approach is somewhat practical for handling simple stowage tasks with limited scale.

[0003] However, loading methods relying on human experience and single heuristic rules have significant drawbacks. When dealing with thousands of cargo holds and containers, the computational complexity of manual loading increases dramatically, making it impossible to arrive at an optimal solution that satisfies multiple constraints within a reasonable timeframe. Single heuristic algorithms are prone to getting trapped in local optima, lacking the ability to escape local extremes during the search process. This results in a significant discrepancy between the final loading solution and the theoretical optimal solution in terms of overall layout balance, yard turnover rate, and ship stability. While some existing loading methods based on evolutionary algorithms possess global search capabilities, their initial solution construction remains relatively crude, resulting in numerous infeasible solutions in the early stages of the search, low convergence efficiency, and an inability to meet the timeliness requirements of actual terminal operations for solution generation.

[0004] In container loading and stowage problems, the actual operations in the yard change frequently. Real-time factors such as quay crane scheduling and truck turnaround can cause deviations between the planned stowage sequence and the actual execution sequence. Static stowage schemes are insufficient to adapt to such disturbances. Meanwhile, each port terminal accumulates a large amount of historical stowage data in actual operation. This data implicitly contains proven container combination patterns and space allocation experience. However, conventional methods lack an effective mechanism to systematically integrate successful grouping experiences from historical cases into the initial solution construction process. This causes the starting point of the optimization search to deviate from the high-quality solution region, requiring subsequent iterations to consume significant computational resources to approach the feasible space. Summary of the Invention

[0005] This paper proposes a container loading and stowage method that can integrate historical loading experience and approach the global optimal solution in a hierarchical and progressive manner. Through a multi-stage coupled optimization strategy, it can construct a high-quality initial solution while maintaining an effective search of complex constraint space, thereby improving the overall performance of the loading scheme in terms of ship stability, structural strength and on-site operational adaptability.

[0006] To achieve the above objectives, this invention provides the following technical solution: This invention provides a container loading and stowage method based on a multi-stage optimization algorithm, comprising: acquiring the target vessel's hold layout data and the physical attribute data, destination data, and loading / unloading time window data of the containers to be stowed; and generating an initial stowage plan based on the above data using a multi-stage optimization algorithm. This multi-stage optimization algorithm decomposes the complex stowage problem into three interconnected stages in a progressive manner, balancing solution efficiency and solution quality.

[0007] In the first stage, pattern matching and pre-grouping of containers to be loaded are performed based on a historical loading case library. Specifically, container combination features are extracted from historical loading records to construct feature vectors for the current batch of containers, and the cosine similarity between the two is calculated. Historical grouping methods with similarity higher than a preset threshold are included in the candidate grouping set. This process fully utilizes historical successful loading experience, quickly forms a high-quality initial grouping framework, significantly compresses the search space for subsequent optimization, and avoids the computational overhead of blindly combining from scratch.

[0008] The second stage, based on the candidate group set and incorporating ship stability and strength constraints, employs a tabu search algorithm to refine container allocation within each group, obtaining a locally optimal loading scheme. Specifically, by defining neighborhood operations such as swapping the positions of different containers within two groups or moving a single container to an empty space, the algorithm aims to minimize the maximum stacking height and the number of container turnovers. It iterative optimization is performed while satisfying strict safety constraints such as heel angle, trim angle, initial metacentric height, deck allowable load, and hull bending moment. A tabu table stores summary information of visited solutions to prevent the search process from getting stuck in loops, thus converging to a locally optimal loading structure with good stacking characteristics and operational convenience within a finite number of iterations.

[0009] The third stage uses the locally optimal loading scheme obtained in the second stage as the elite individuals in the initial population. A genetic algorithm based on non-dominated sorting is used for multi-objective global optimization to output the final loading scheme. The loading scheme is encoded as a chromosome reflecting the correspondence between cargo space and containers. After generating an initial population containing this locally optimal solution, the quality of individuals is evaluated using non-dominated sorting and crowding calculation. The population is evolved through selection, crossover, and mutation operations. Crossover swaps the container allocation of the corresponding cargo space columns of parent chromosomes, and mutation randomly changes the containers corresponding to cargo spaces to legal replacement values, thus maintaining the diversity and feasibility of the population. This evolutionary process is repeated until the termination condition is met. Finally, the fitness of the Pareto optimal solution set is comprehensively evaluated based on preset weights to select the overall optimal final loading scheme. This stage overcomes the limitations of single-objective optimization, achieving an effective balance among multiple mutually constraining objectives such as stability, intensity, and operational efficiency, resulting in a more applicable and practical scheme.

[0010] Preferably, during the execution of the multi-stage optimization algorithm, real-time port operation data such as the progress of quay crane operations and the real-time location of containers in the yard are dynamically accessed. When the deviation between the actual loading and unloading sequence and the planned sequence in the final stowage scheme exceeds a preset deviation threshold, the third-stage genetic algorithm is triggered to re-optimize the current scheme. Through this dynamic adjustment mechanism, the present invention can effectively cope with uncertainties in the operation process, ensure that the stowage scheme is always consistent with the actual situation on site, and avoid yard congestion and ship delays caused by deviations from the plan.

[0011] As a further improvement to this invention, adaptive crossover and mutation probabilities are introduced into the genetic algorithm. Operator parameters are dynamically adjusted based on the current population's non-dominant rank distribution and individual crowding distance, enabling the algorithm to maintain strong exploratory capabilities in the early stages of evolution and focus on fine-grained local searches in the later stages, thereby accelerating global convergence and obtaining a more uniformly distributed Pareto front. In the tabu search, the tabu table uses a hash table to store the hash values ​​of the cargo space allocation matrix and the objective function value. The tabu object is a complete loading scheme, and the tabu length is dynamically set according to the total number of containers to be loaded. This mechanism reduces storage overhead and enhances the algorithm's ability to escape local optima.

[0012] This invention also provides a container loading and stowage system based on a multi-stage optimization algorithm, including a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it implements the steps of the container loading and stowage method based on the multi-stage optimization algorithm described above. This system deeply integrates case-based reasoning pre-grouping, tabu search that strengthens local search, and genetic algorithms that consider multiple objective trade-offs. It can quickly generate stowage schemes that strictly satisfy stability and strength constraints while taking into account yard turnover rate and loading and unloading efficiency, and has the ability to respond in real time and re-optimize the scheme under operational disturbances.

[0013] The technical effects and advantages provided by the present invention in the above technical solution are as follows: The first stage utilizes a historical loading case database to perform pattern matching and pre-grouping of containers to be loaded, constructing feature vectors for the current batch of containers. Cosine similarity is then calculated between these vectors and container combination features from historical loading records, such as weight distribution, destination combinations, and loading / unloading sequences, to filter out candidate groups with matching scores exceeding a preset threshold. This historical case embedding mechanism changes the conventional algorithm's search mode, which starts from random or simple rule-based solutions. The optimization search directly unfolds within the neighborhood of validated good solutions, reducing the overhead of generating and evaluating a large number of infeasible solutions in the early stages. This ensures the initial solution has good physical rationality and accelerates the convergence of the subsequent fine-tuning process.

[0014] The third stage uses the locally optimal loading scheme obtained from tabu search as the initial population members of a genetic algorithm based on non-dominated sorting. These individuals, along with randomly generated individuals, participate in multi-objective global optimization. Real-time port operation data is dynamically integrated during execution. When the deviation between the actual loading / unloading sequence and the planned sequence exceeds a preset threshold, the genetic algorithm is triggered to re-optimize. This coupling design of local optimization results and global optimization overcomes the limitations of a single tabu search (limited scope) and the initial blindness of a single genetic algorithm. It also endows the loading scheme with responsiveness to the progress of quay crane operations and real-time changes in the yard, exhibiting self-adjusting characteristics to operational disturbances. This results in a more balanced allocation of multiple objectives, including ship stability constraints, strength constraints, and loading / unloading efficiency, adapting to the dynamically changing pace of terminal production. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of this application or the prior art, the drawings used in the embodiments will be briefly introduced below. Obviously, the drawings described below are only some embodiments recorded in this invention. For those skilled in the art, other drawings can be obtained based on these drawings.

[0016] Figure 1 This is a flowchart of a container loading and stowage method based on a multi-stage optimization algorithm; Figure 2 This is a flowchart for generating container grouping schemes based on historical loading case matching; Figure 3 This is a flowchart of a container loading optimization method based on a non-dominated sorting genetic algorithm; Figure 4 These are probability density curves of weight distribution for different types of containers; Figure 5 It is a curve showing the actual loading and unloading sequence deviation as a function of the number of containers already loaded onto the ship. Detailed Implementation

[0017] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0018] See Figure 1This invention provides a container loading and stowage method based on a multi-stage optimization algorithm, including acquiring the target vessel's hold layout data and the physical attribute data, destination data, and loading / unloading time window data of the containers to be loaded. Based on the acquired physical attribute data, destination data, and loading / unloading time window data, an initial stowage plan is generated through a multi-stage optimization algorithm. The multi-stage optimization algorithm consists of three stages: First, based on a historical stowage case library, pattern matching and pre-grouping of the containers to be loaded are performed to obtain a candidate group set with a matching degree higher than a preset threshold. Second, based on the obtained candidate group set, combined with ship stability constraints and strength constraints, a tabu search algorithm is used to allocate hold space to the containers within the group to obtain a locally optimal stowage plan. Third, the locally optimal stowage plan is used as part of the initial population, and a multi-objective global optimization is performed using a genetic algorithm based on non-dominated sorting to output the final stowage plan. The hold layout data includes the column number, layer number, and corresponding load-bearing capacity information of each compartment of the vessel. The physical attribute data includes the container weight, size, and type. The destination data includes the unloading port sequence of each container. The loading and unloading time window data includes the estimated arrival and departure times. Ship stability constraints include limits on heel angle, trim angle, and lower limit on initial metacentric height; strength constraints include allowable deck loads, allowable bilge loads, and hull bending moment limits. The final loading scheme obtained through a multi-stage optimization algorithm satisfies the above constraints and optimizes multiple objectives.

[0019] Example 1: In specific implementation, please refer to Figure 2 A historical stowage case database is established. This database contains multiple historical stowage records, each corresponding to a completed container loading and stowage task. Each historical stowage record includes container combination characteristics, including weight distribution, destination combination, and loading / unloading sequence. Weight distribution refers to the percentage of containers in each weight range within the historical stowage record; destination combination refers to the percentage of containers in each port of discharge; and loading / unloading sequence refers to the order in which the containers were loaded or unloaded within the historical stowage record.

[0020] Extract container combination features for each historical loading record from the historical loading case database. For each historical loading record, extract the weight distribution to construct a weight distribution vector, and the vector's _i_th _j ... Each component is , indicating that it falls into the . The proportion of containers in each weight range to the total number of containers in the historical loading record; extract destination combinations to construct a destination combination vector, the vector's first... Each component is This indicates that the port of destination is the first The proportion of containers at each port to the total number of containers in the historical loading record; extract the loading and unloading sequence to construct a loading and unloading sequence vector, the vector's first... Each component is , indicating the first The normalized sequence number values ​​of containers loaded and unloaded at each sequential position in the historical stowage record are used. The weight distribution vector, destination combination vector, and loading / unloading sequence vector are concatenated to form the container combination feature vector of the historical stowage record.

[0021] Using the physical attribute data, destination data, and loading / unloading time window data of the containers to be loaded, a container feature vector for the current batch is constructed. The physical attribute data includes container weight, the destination data includes the order of unloading ports, and the loading / unloading time window data includes the estimated arrival and departure times. Based on the container weights of the containers to be loaded, the proportion of containers in each weight range is calculated, generating a weight distribution vector for the current batch. The vector's _i_th_i ... Each component is This indicates that the container in the current batch awaiting loading has fallen into the first category. The proportion of containers in each weight range to the total number of containers awaiting loading. Based on the destination data of the containers awaiting loading, the proportion of containers at each port of discharge is calculated, generating a destination combination vector for the current batch. The vector's _i_th__... Each component is This indicates that the destination port of the containers awaiting loading in the current batch is the [number]. The proportion of containers at each port to the total number of containers awaiting loading. Based on the loading and unloading time window data of the containers awaiting loading, the estimated arrival time of each container is extracted, and the containers are sorted according to time to generate a loading and unloading order vector. The first position of the vector is... Each component is , indicating the first The normalized sequence number of the containers at each sequential position. The weight distribution vector, destination combination vector, and loading / unloading sequence vector are concatenated to form the container feature vector of the current batch.

[0022] Calculate the cosine similarity between the feature vector of the current batch of containers and the feature vector of the container combination in each historical loading record. The formula for calculating the cosine similarity is as follows:

[0023] in, This represents the cosine similarity between the feature vector of the current batch of containers and the feature vector of a container combination in a historical loading record. Indicates the total number of weight ranges; This represents the weight distribution vector of the current batch. The value of each component; The first element in the weight distribution vector of historical load records The value of each component; Indicates the total number of ports of discharge; This represents the first element in the current batch destination combination vector. The value of each component; This represents the first term in the destination combination vector of historical load records. The value of each component; The dimension of the loading / unloading sequence vector is numerically equal to the total number of containers to be loaded; This represents the first element in the current batch loading / unloading sequence vector. The value of each component; This represents the first element in the historical loading and unloading sequence vector. The value of each component.

[0024] The container grouping methods corresponding to historical loading records with a cosine similarity greater than a preset similarity threshold are used as candidate grouping schemes. The preset similarity threshold is set based on the statistical distribution of historical loading records in the historical loading case library. The setting method is as follows: for all historical loading records in the historical loading case library, calculate the cosine similarity between the container combination feature vectors of every two different historical loading records to obtain a similarity set; take the 90th percentile of the similarity set as the preset similarity threshold.

[0025] All obtained candidate grouping schemes are deduplicated. Specifically, each candidate grouping scheme is represented as a group identifier sequence. The group identifier sequences of multiple candidate grouping schemes are compared. If two candidate grouping schemes have the same group identifier sequence, one candidate grouping scheme is retained, and the other is deleted. Schemes with the same grouping in the deduplicated candidate grouping schemes are merged. Specifically, multiple candidate grouping schemes with the same group identifier are merged into one grouping scheme. The groups in this merged scheme include all containers from the merged schemes. After merging, a set of candidate groups with a matching degree higher than a preset threshold is generated. The matching degree is determined by the cosine similarity of the candidate grouping schemes, and the cosine similarity value is the matching degree value.

[0026] Example 2: In practice, an initial cabin allocation matrix is ​​established based on the cabin area corresponding to each group in the candidate grouping set. A group corresponds to a continuous cabin area on the ship, and the cabin area consists of multiple cabins, each with a cabin column number and a cabin deck number. The rows of the initial cabin allocation matrix represent the cabin column numbers, and the matrix's rows represent the cabin column numbers. The corresponding cabin class and column number is All cabin classes; the columns of the initial cabin allocation matrix represent cabin class numbers, and the matrix's first column represents the cabin class number. The corresponding cabin deck number is All cabin classes. Elements of the initial cabin allocation matrix. Indicates that it is located at the th Line number The container number assigned to the slot in the column, when When the corresponding cargo space is not allocated any containers, It is assigned a vacant slot identifier, which takes the value -1.

[0027] Ship stability constraints and strength constraints are treated as hard constraints. Ship stability constraints include limits on heel angle, trim angle, and the lower limit of the initial metacentric height. Strength constraints include allowable deck loads, allowable bilge loads, and hull moment limits. After each change to the loading scheme, the corresponding heel angle, trim angle, initial metacentric height, deck load, bilge load, and hull moment are calculated and compared with the limits specified in the ship stability and strength constraints. If any of the following occurs: heel angle exceeds the heel angle limit; trim angle exceeds the trim angle limit; initial metacentric height is below the lower limit of the initial metacentric height; deck load exceeds the allowable deck load; bilge load exceeds the allowable bilge load; or hull moment exceeds the hull moment limit, the loading scheme is deemed to violate the hard constraints and is discarded as a candidate solution.

[0028] The neighborhood operation in tabu search is defined as two operations. The first neighborhood operation is to swap the slot positions of two different containers located in the same group. The first container and the second container are selected from the same group; the first container is currently assigned to slot number [missing information]. The second container is currently allocated to the slot number [missing information]. The first container will be allocated to the designated space. The second container was allocated to the hold. The second type of neighborhood operation involves moving a single container to an available slot within the same group. A third container is selected from the same group; this third container is currently assigned to slot number [missing information]. Select an available cabin class within the same group. The cabin class number for the available cabin class is: The third container was allocated to the hold. and cabin Set to empty cabin.

[0029] The optimization objectives are to minimize the maximum stacking height and the number of container repositioning operations. The maximum stacking height is defined as the maximum number of container stacking layers across all cargo slots. The number of container repositioning operations is defined as the total number of times other containers need to be moved during unloading because the target container is pressed down by other containers. An objective function is established, and its value is expressed as:

[0030] in, Represents the objective function value; This indicates the maximum stacking height of the current loading scheme; This indicates the number of times the container needs to be flipped in the current loading scheme; The weighting coefficient represents the maximum stacking height value, and its value ranges from (0,1). The weighting coefficient represents the number of times the box was searched. Weighting coefficients and The method for obtaining this information is as follows: it is set according to the relative importance that the ship operator places on the stacking height control target and the container turnover number control target. When the ship operator places more emphasis on controlling the maximum stacking height, A value greater than 0.5 is used when the ship operator places greater emphasis on controlling the number of times containers are re-tucked. Take a value less than 0.5.

[0031] Iteratively perform neighborhood operations. Starting from the current loading scheme, generate all possible neighborhood solutions. Perform hard constraint verification on each neighborhood solution. If a neighborhood solution violates the hard constraints, remove it from the candidate set. Calculate the objective function value of each neighborhood solution that passes the hard constraint verification, and select the neighborhood solution with the smallest objective function value as the candidate optimal neighborhood solution. Compare this candidate optimal neighborhood solution with the summary information of visited solutions stored in the tabu list.

[0032] A tabu list is set up to store summary information of solutions visited in the most recent iterations. The summary information includes the hash value of the space allocation matrix and the objective function value. The hash value of the space allocation matrix is ​​obtained by operating a secure hash algorithm on the element sequence of the space allocation matrix, and is used to uniquely identify a loading scheme. The tabu object is a complete loading scheme. When the hash value of the space allocation matrix of a candidate optimal neighborhood solution is equal to any hash value stored in the tabu list, and the objective function value of the candidate optimal neighborhood solution is the same as the objective function value of the corresponding record's hash value, that candidate optimal neighborhood solution is determined to be tabu. The tabu length is set to a dynamic value positively correlated with the total number of containers to be loaded. The formula for calculating the tabu length is... ,in Forbidden length value, Forbidden length coefficient, The value was set to 0.2 based on empirical statistics. This represents the total number of containers awaiting loading. Indicates rounding up. Forbidden length value. It increases with the total number of containers waiting to be loaded.

[0033] In each iteration, if the candidate optimal neighborhood solution is in a tabu state, the neighborhood solution with the smallest objective function value is selected from the remaining non-tabulated neighborhood solutions as the selected solution for the current iteration; if all neighborhood solutions are in a tabu state, the solution that was first added to the tabu list in the tabu state is selected as the selected solution for the current iteration. The selected solution is updated to the current loading scheme, and the summary information of the selected solution is added to the tabu list. If the number of summary information stored in the tabu list exceeds the tabu length value, the earliest added summary information is removed. The selected solution and its update number are recorded after each iteration.

[0034] The search stops when the preset maximum number of iterations is reached, or when the objective function value shows no improvement after several consecutive iterations. The determination of no improvement after several consecutive iterations is as follows: count the number of iterations within a monitoring window (set to 100 iterations). If the minimum objective function value remains unchanged within the most recent 100 iterations, then the objective function value is considered to have no improvement. After stopping the search, the current optimal solution is output as the locally optimal loading scheme. The current optimal solution is the one with the smallest objective function value among all recorded selected solutions.

[0035] Example 3: In specific implementation, please refer to Figure 3 The loading information corresponding to the locally optimal loading scheme obtained in the second stage is encoded as a chromosome. The locally optimal loading scheme consists of the container allocation results for each hold. During encoding, a gene segment position is assigned to each container-carrying hold on the ship, and the order of the gene segments is consistent with the order of the hold numbers. A chromosome consists of multiple gene segments, each gene segment corresponding to a hold, and the value stored in the gene segment is the container number allocated to that hold. If a hold is not allocated a container, the corresponding gene segment stores an empty hold identifier, which is set to -1.

[0036] Multiple individuals, including those representing locally optimal loading schemes, are randomly generated to form an initial population. The initial population is generated as follows: the chromosome encoded by the locally optimal loading scheme is used as one individual in the initial population. The remaining individuals are generated by iterating through the gene segments corresponding to all cargo slots for each individual to be generated, and randomly assigning a valid container number to each gene segment. The random selection range includes containers within the current group that have not yet been assigned to other cargo slots for this individual, ensuring that the same container is not repeatedly assigned to multiple cargo slots and maintaining the validity of the loading scheme. The initial population size is set to a fixed value, determined based on the scale of the loading problem: when the total number of containers to be loaded is less than or equal to 500, the initial population size is set to 100; when the total number of containers to be loaded is greater than 500, the initial population size is set to 200.

[0037] The initial population is subjected to a non-dominated sort. During the non-dominated sort, multiple optimization objective values ​​are calculated for each individual in the population. The optimization objectives include minimizing the maximum stack height and minimizing the number of box flips. Based on the multiple optimization objective values, Pareto dominance is compared between individuals. If individual A is not inferior to individual B in all optimization objectives, and is strictly superior to individual B in at least one of the optimization objectives, then individual A dominates individual B. All individuals in the population are divided into multiple non-dominated levels. The set of individuals with non-dominated level 1 is the set of individuals in the population that are not dominated by any other individual. After removing individuals from the first non-dominated level of the population, individuals that are not dominated by any remaining individuals are identified as non-dominated level 2. This process is repeated until all individuals are assigned to the corresponding non-dominated level. Calculate the crowding distance for each individual within its respective non-dominated level. The crowding distance is calculated as follows: for each optimization objective, sort the individuals within that level in ascending order according to the optimization objective value, and set the crowding distance between the first and last individuals in the sort to a maximum value; for the middle individuals, the crowding distance is the sum of the normalized differences between the two adjacent individuals before and after that individual at each optimization objective.

[0038] The tournament selection process involves choosing parent individuals from the population. The specific steps are as follows: Each time, two individuals are randomly selected from the population. Their non-dominance levels are compared. If the non-dominance levels are different, the individual with the lower non-dominance level is selected. If the non-dominance levels are the same, their crowding distance is compared, and the individual with the larger crowding distance is selected. This selected individual is then used as the parent individual. This process is repeated until the number of selected parent individuals meets the requirements for the crossover operation.

[0039] A single-point crossover operation is used to swap the container allocations in corresponding slot columns of two parent chromosomes. During single-point crossover, a crossover point is randomly selected, defined between two gene segments. The gene segment sequence before the crossover point on the first parent chromosome is concatenated with the gene segment sequence after the crossover point on the second parent chromosome to generate the first offspring chromosome; the gene segment sequence before the crossover point on the second parent chromosome is concatenated with the gene segment sequence after the crossover point on the first parent chromosome to generate the second offspring chromosome. After the crossover operation, the offspring chromosomes undergo validity correction: the system checks for duplicate or missing container numbers. If duplicate allocations are found, the redundant allocations in the duplicate slots are removed, and the missing container slots are randomly assigned valid container numbers that have not yet appeared.

[0040] The site mutation operation randomly changes the container number of a specific cargo hold to a valid replacement value. During site mutation, for each gene segment in a chromosome, an adaptive mutation probability is used to determine whether a mutation has occurred. If a mutation is determined, the cargo hold corresponding to that gene segment is selected, and a container number is randomly chosen from the container numbers in the current group that do not appear in other cargo holds on that chromosome as a replacement value to replace the original value of the gene segment, thus ensuring the validity of the mutated loading scheme. All offspring individuals generated by the crossover and mutation operations constitute the next generation population.

[0041] Repeat the steps of non-dominated sorting, crowding distance calculation, tournament selection, single-point crossover, and site mutation until the maximum number of generations is reached. The maximum number of generations is set based on the initial population size and the total number of containers to be loaded: when the total number of containers to be loaded is less than or equal to 500, the maximum number of generations is set to 300; when the total number of containers to be loaded is greater than 500, the maximum number of generations is set to 500.

[0042] The Pareto optimal solution set is extracted from the last generation of the population. This set consists of all individuals with a non-dominant level of 1. The objective function values ​​of each solution in the Pareto optimal solution set are weighted and summed according to a predefined weight vector to calculate the overall fitness. The formula for calculating the overall fitness is:

[0043] in, This represents the overall fitness value of a particular solution; This indicates the maximum stacking height value corresponding to the solution; This indicates the number of times the box was searched for corresponding to the solution; This represents the weighting coefficient of the maximum stack height value in the overall fitness score, with a value ranging from 0 to 1; This represents the weighting coefficient of the number of times the box was searched in the overall fitness score. Weighting coefficients and Based on the ship operator's preferences, when controlling the maximum stacking height takes precedence, The value should be greater than 0.5; when the control over the number of times the box is turned over takes priority, The value is less than 0.5. The solution with the smallest overall fitness value is selected as the final load arrangement scheme.

[0044] Adaptive crossover and mutation probabilities are used in crossover and mutation operations. The adaptive crossover probability is dynamically adjusted based on the current distribution of non-dominant ranks in the population. The adjustment method is as follows: count the number of individuals in the current population at the first non-dominant rank. Let the total number of individuals in the population be... Then adaptive crossover probability The calculation is expressed as ,in To maximize the crossover probability, the value is set to 0.9, following the standard genetic algorithm setting. The minimum crossover probability is set to 0.6. When the proportion of individuals with the first non-dominant rank in the population is high, the crossover probability decreases to preserve the genes of high-quality individuals; when the proportion of individuals with the first non-dominant rank is low, the crossover probability increases to enhance the population's exploration ability.

[0045] The adaptive mutation probability is dynamically adjusted based on the coefficient of variation of the crowding distance between individuals. The adjustment method is as follows: calculate the average crowding distance of all individuals in the current population. and standard deviation Coefficient of variation of crowding distance Adaptive mutation probability From the basic mutation probability Determined by both the coefficient of variation and the variance coefficient, the calculation method is as follows: Basic mutation probability The value is set to 0.05, determined based on the commonly used mutation rate range in genetic algorithms. When the coefficient of variation of the crowding distance increases, it indicates that the distribution differences between individuals are becoming larger. In this case, the mutation probability is increased to maintain diversity. When the coefficient of variation decreases, the mutation probability is reduced to minimize the destruction of good solutions.

[0046] Example 4: In practice, physical attribute data includes container weight, container dimensions, and container type. Container weight is expressed in kilograms and obtained from the terminal weighing system or shipper declaration data. Container dimensions include 20-foot containers, 40-foot containers, and corresponding height information. The external length of a 20-foot container is approximately 6058 mm, and the external length of a 40-foot container is approximately 12192 mm. Height information includes two specifications: standard height 2591 mm and high cube height 2896 mm. Container types include general dry cargo containers, refrigerated containers, and dangerous goods containers. General dry cargo containers are used to load goods without special temperature control requirements, refrigerated containers are used to load goods requiring a specific temperature environment, and dangerous goods containers are used to load dangerous goods classified according to the International Maritime Dangerous Goods Code. Container type data is obtained from the booking system and stored in data records as type codes.

[0047] The destination data includes the port-of-call sequence for each container, which is the sequential numbering of the ports the vessel calls at in order of port call. This sequence is determined by the shipping company's pre-published voyage plan, which lists all ports the vessel calls at sequentially from its port of origin. Each container's destination data records its port-of-call sequence number, starting from 1 and incrementing. 1 represents the first port of call, and higher numbers indicate later ports of call. The port-of-call sequence is extracted from the booking information and linked to the corresponding containers.

[0048] The loading and unloading time window data includes the estimated arrival and departure times for each container. The estimated arrival time refers to the expected arrival time of the vessel carrying the container at the port of discharge, and the estimated departure time refers to the expected departure time of the vessel from the port of discharge. These estimated arrival and departure times are obtained from the port information system, which generates the time window data based on the vessel's voyage plan, speed, and port berth arrangements. Each container's loading and unloading time window data corresponds one-to-one with the order of the ports of discharge. If a container has multiple time windows at multiple ports of discharge, only the set of time window data corresponding to the port of discharge is recorded.

[0049] Ship stability constraints include limits on heel angle, pitch angle, and the lower limit on initial metacentric height. The heel angle limit is determined based on the maximum permissible heel angle specified in the ship's stability manual, and the pitch angle limit is also determined based on the maximum permissible pitch angle specified in the ship's stability manual, which is provided by the ship's design firm or classification society. The lower limit on initial metacentric height is set to be no less than 0.15 meters, determined according to the minimum requirements for initial metacentric height in the International Maritime Organization's Complete Stability Code 2008. The initial metacentric height is calculated using the following formula during the evaluation of loading schemes:

[0050] in, This represents the initial stable height value, in meters. This indicates the height of the ship's transverse metacenter from the baseline, in meters, and is obtained from the ship's hydrostatic parameters table based on the current displacement. This represents the height of the ship's center of gravity, in meters, and is calculated from the weight of each container, the height of each container's center of gravity, and the ship's own weight and center of gravity position. The calculated initial metacentric height must be greater than or equal to the lower limit of the initial metacentric height; otherwise, the corresponding loading scheme will be deemed to violate ship stability constraints.

[0051] Strength constraints include allowable deck loads, allowable bilge loads, and hull bending moment limits. Allowable deck loads refer to the maximum weight per unit area that each deck level can withstand, measured in tons per square meter (t / m²), and are obtained from the ship's structural drawings. Allowable bilge loads refer to the maximum weight per unit area that each compartment's bottom can withstand, measured in tons per square meter (t / m²), and are also obtained from the ship's structural drawings. Hull bending moment limits include the maximum allowable sagging moment and the maximum allowable camber moment, measured in kilonewton-meters (kN·m), and are obtained from the ship's loading manual. During the evaluation of the stowage plan, the total weight of the containers loaded in each compartment is divided by the compartment's floor area to obtain the bilge load value, which is then compared with the allowable bilge load. Similarly, the total weight of the containers loaded on each deck is divided by the deck area to obtain the deck load value, which is also compared with the allowable deck load. Simultaneously, the hull bending moment value generated by the overall ship weight distribution is calculated and compared with the hull bending moment limits.

[0052] During the second phase of tabu search, after each neighborhood solution is generated, the corresponding heel angle, trim angle, initial metacentric height, deck load, bilge load, and hull bending moment are calculated. Each calculated value is then compared with the limits specified in the ship's stability and strength constraints. Neighborhood solutions that violate any constraint are marked as infeasible and are not included in the subsequent calculation and comparison of the objective function value, nor are they incorporated into the candidate solution set.

[0053] In the third stage of the genetic algorithm execution based on non-dominated sorting, after decoding the loading scheme for each individual in the population, the values ​​of heel angle, pitch angle, initial metacentric height, deck load, bilge load, and hull bending moment are calculated and compared with the limits specified in the ship's stability and strength constraints. Individuals that violate any constraint are assigned a penalty term, a value much larger than the normal objective function value, causing this infeasible individual to be ranked lower in the non-dominated sorting and thus naturally eliminated during the evolutionary process. Handling infeasible solutions through a penalty function ensures that the optimization search is conducted within the principal feasible region.

[0054] See Figure 4 The figure shows the weight probability density distribution based on container type in Example 4. The horizontal axis represents the container weight in kilograms, ranging from 5,000 kg to 40,000 kg, and the vertical axis represents the probability density value of the corresponding weight. The legend distinguishes three container types: general dry cargo containers (blue curve), refrigerated containers (orange curve), and dangerous goods containers (red curve).

[0055] The probability density curves show that refrigerated containers have a higher probability density in the lighter weight range, peaking between approximately 13,000 kg and 15,000 kg, indicating that their weight is concentrated in this range. Ordinary dry goods containers have a more even weight distribution, peaking between approximately 22,000 kg and 24,000 kg, with a lower probability density peak than refrigerated containers, showing a heavier and wider distribution. Dangerous goods containers, on the other hand, have a heavier weight distribution, peaking between approximately 28,000 kg and 30,000 kg, with a flatter probability density curve, indicating that dangerous goods containers are heavier and have a wider distribution range.

[0056] Example 5: In practice, the multi-stage optimization algorithm dynamically incorporates real-time port operation data during execution. This real-time port operation data includes the progress of quay crane operations and the real-time location of containers in the yard. Quay crane operation progress data is obtained from the terminal operating system, which records the container numbers currently being loaded / unloaded by each quay crane, a list of completed container numbers, and the percentage of the planned quay crane operation. Real-time container location data in the yard is obtained from the yard management system, which maintains the container's bay number, column number, layer number, and actual movement record within the yard.

[0057] Real-time port operation data is compared with the final stowage plan. The final stowage plan includes the planned loading and unloading sequence for each container, which specifies the order in which containers are hoisted onto the ship during the loading process. The actual loading and unloading sequence is extracted from the real-time port operation data; this actual sequence is the order in which the quay cranes actually complete the container hoisting operations onto the ship. The deviation between the actual and planned loading and unloading sequences is calculated using the following formula:

[0058] in, This indicates the magnitude of the deviation between the actual loading and unloading sequence and the planned loading and unloading sequence. This indicates the total number of containers that have actually been loaded onto the ship; The index number indicates the container that has actually been loaded onto the ship, and the index numbers are arranged in chronological order of the actual loading completion time. Indicates the index number is The priority value of the container in the actual loading and unloading sequence, with the priority value starting from 1 and incrementing sequentially; Indicates the index number is The priority value of the container in the planned loading and unloading sequence; Indicates the index number is The absolute value of the difference between the actual and planned loading / unloading sequence of containers. Deviation value. It reflects the overall degree of deviation between the actual work sequence and the planned sequence.

[0059] The deviation value is compared with the preset deviation threshold. The preset deviation threshold is set to 3.0. The basis for this value is that when the deviation value reaches 3.0, it means that the average sequence deviation of the loaded containers has reached more than 3 positions. At this time, the expected number of container turnings in the yard bay where the containers waiting to be loaded are located may differ significantly from the actual number of container turnings, and the loading plan needs to be re-optimized to correct it.

[0060] When the detected deviation exceeds a preset deviation threshold, the third stage of the genetic algorithm based on non-dominated sorting is triggered. Upon triggering, optimization is re-started from the current actual state. The current actual state includes the current hold information of containers that have been loaded onto the ship, the real-time location information of unloaded containers in the yard, and the current operation progress of the quay crane. The hold allocation results of containers that have been loaded onto the ship are fixed and no longer participate in the gene segment mutation operation in the genetic algorithm. The unloaded containers are treated as a set of containers to be reloaded, generating a new initial population. Each individual in the new initial population has the hold allocation of loaded containers fixed in its gene segment, and only the gene segment of unloaded containers is randomly initialized. The non-dominated sorting, crowding distance calculation, tournament selection, single-point crossover, and site mutation operations in the third stage are executed until the maximum number of generations is reached, generating the updated final loading scheme.

[0061] In the third stage of the genetic algorithm, the adaptive crossover probability is calculated based on the distribution ratio of individuals at each non-dominant level in the current population. Specifically, the calculation method is as follows: count the number of individuals in the current population at non-dominant level 1, denoted as... Count the number of individuals in the current population that are at the non-dominant level 2, and denote them as follows: The total number of individuals in the population is Adaptive crossover probability The calculation method is as follows ,in To maximize the crossover probability, the value is set to 0.9. The minimum crossover probability is set to 0.6. The crossover probability decreases when the proportion of individuals with a non-dominance level of 1 in the population excluding individuals with a non-dominance level of 2 increases; the crossover probability increases when the proportion of individuals with a non-dominance level of 1 decreases.

[0062] In the third stage of the genetic algorithm, the adaptive mutation probability is dynamically adjusted based on the coefficient of variation of the crowding distance between individuals. Specifically, the adjustment method involves calculating the standard deviation of the crowding distance between all individuals in the current population, denoted as... Calculate the average crowding distance of all individuals in the current population, denoted as . Coefficient of variation of crowding distance Adaptive mutation probability The calculation method is as follows ,in The base mutation probability is set to 0.05, a conservative value within the commonly used range for mutation probabilities in genetic algorithms, to reduce solution quality fluctuations caused by excessive mutation. As the mutation coefficient increases, the crowding distribution among individuals becomes more differentiated, leading to an increase in the adaptive mutation probability to enhance population diversity; conversely, as the mutation coefficient decreases, the adaptive mutation probability is reduced to maintain a good solution structure.

[0063] In the second stage of tabu search, the tabu table uses a hash table to store the summary information of visited solutions. The summary information includes the hash value of the space allocation matrix and the objective function value. The hash value of the space allocation matrix is ​​obtained by operating a secure hash algorithm on the element sequence of the space allocation matrix; the hash value is used to uniquely identify a loading scheme. After generating a neighborhood solution, the elements of the space allocation matrix corresponding to the neighborhood solution are arranged in row-major order into a one-dimensional array. This one-dimensional array is input into the secure hash algorithm, which outputs a hash value of fixed length. The hash value of the neighborhood solution is compared one by one with the hash values ​​already stored in the tabu table. If the hash value of the neighborhood solution is the same as the hash value of a record in the tabu table, and the objective function value of the neighborhood solution is the same as the objective function value of that record in the tabu table, then the neighborhood solution is determined to be tabu. The tabu object is a complete loading scheme, and the tabu length is set according to the total number of containers to be loaded. The tabu length is calculated as follows: ,in Forbidden length value, Forbidden length coefficient, The system is segmented based on the total number of containers to be loaded: when the total number of containers to be loaded is less than or equal to 200... The value is 0.15; when the total number of containers to be loaded is greater than 200 and less than or equal to 500, The value is 0.20; when the total number of containers to be loaded is greater than 500, The value is 0.25. This represents the total number of containers awaiting loading. This indicates rounding up. The larger the total number of containers to be loaded, the greater the forbidden length value, thus effectively avoiding cyclic searches within a larger search space.

[0064] The container loading and stowage system based on a multi-stage optimization algorithm includes a memory, a processor, and a computer program stored in the memory and running on the processor. When the processor executes the computer program, it performs the following steps: acquiring the target ship's hold layout data and the physical attribute data, destination data, and loading / unloading time window data of the containers to be stowed; implementing the first stage of pattern matching and pre-grouping, which involves extracting container combination features from a historical stowage case library, constructing the feature vector of the current batch of containers, calculating cosine similarity, and using the container grouping methods corresponding to historical stowage records with cosine similarity greater than a preset similarity threshold as candidate grouping schemes, and generating a candidate group set through deduplication and merging; and implementing the second stage of tabu search hold allocation, which involves establishing an initial hold allocation matrix, defining neighborhood operations as exchanging the hold positions of different containers within the same group or moving a single container to an empty hold within the same group, using ship stability constraints and strength constraints as hard constraints, and minimizing the maximum stacking height and the number of container flipping operations. To optimize the objective, a hash table is used to store the summary information of visited solutions. Tabu search iterations are performed to output locally optimal loading schemes. The third stage implements a multi-objective global optimization step based on a non-dominated sorting genetic algorithm. This involves encoding the locally optimal loading scheme as a chromosome, randomly generating an initial population, performing non-dominated sorting and crowding distance calculations, generating the next generation population through tournament selection, single-point crossover, and site mutation, and using adaptive crossover and adaptive mutation probabilities to repeat the evolution until the maximum number of generations is reached. The final loading scheme is selected from the Pareto optimal solution set based on the weight vector. The system dynamically accesses the quay crane operation progress and real-time container location data in the yard, calculates the deviation value and compares it with a preset deviation threshold. When the deviation value exceeds the preset deviation threshold, the third stage genetic algorithm is triggered to re-optimize and output the updated final loading scheme.

[0065] See Figure 5 In the graph, the horizontal axis represents the number of containers loaded onto the ship, in TEUs, ranging from 0 to 600. The vertical axis represents the deviation (DEV) between the actual loading / unloading sequence and the planned loading / unloading sequence, ranging from 0 to 3.2. The solid blue line curve represents the trend of the actual deviation as the number of containers loaded onto the ship changes, while the dashed red line represents the preset deviation threshold, which is fixed at 3.0.

[0066] The curve shows that the actual deviation gradually increases with the number of loaded containers. When the number of loaded containers is approximately 0 to 150, the deviation remains between 0.5 and 0.8, indicating that the actual loading and unloading sequence deviates little from the planned sequence, and the loading operation is relatively stable. Subsequently, the deviation increases slowly, and within the range of 150 to 400 TEU, the fluctuation range gradually increases, reaching approximately 1.5, indicating that the loading and unloading sequence begins to show some deviation but is still within a controllable range.

[0067] Once the number of loaded containers exceeds 400 TEU, the rate of deviation increases significantly, and the curve tends to approach the preset deviation threshold of 3.0. In the range of approximately 550 to 600 TEU, the actual deviation repeatedly approaches but does not significantly exceed the threshold. This stage reflects that the deviation between the loading / unloading sequence and the planned sequence is close to the preset threshold, posing a potential risk to the number of container re-handling operations and operational efficiency.

[0068] According to the deviation amount and preset deviation threshold settings in Example 5, when the deviation amount reaches or exceeds 3.0, a re-optimization based on the non-dominated sorting genetic algorithm needs to be triggered to adjust the loading scheme. The figure shows that during this operation, although the actual deviation amount gradually increases, it has not yet exceeded the threshold, indicating that although the current operation has a deviation, it has not yet reached the critical point to trigger re-optimization and is still within the range of dynamic monitoring and adjustment.

[0069] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A container loading and stowage method based on a multi-stage optimization algorithm, characterized in that, Includes the following steps: Obtain the target vessel's hold layout data and the physical attribute data, destination data, and loading / unloading time window data of the containers to be loaded; Based on the physical attribute data, destination data, and loading / unloading time window data, an initial loading scheme is generated through a multi-stage optimization algorithm. The multi-stage optimization algorithm includes: In the first stage, the containers to be loaded are pattern matched and pre-grouped based on the historical loading case library to obtain a set of candidate groups with a matching degree higher than a preset threshold. In the second stage, based on the candidate group set, the tabu search algorithm is used to allocate space to containers within the group, taking into account ship stability constraints and strength constraints, to obtain a locally optimal loading scheme. In the third stage, the locally optimal loading scheme is used as the initial population, and a multi-objective global optimization is performed using a genetic algorithm based on non-dominated sorting to output the final loading scheme.

2. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 1, characterized in that, The first stage involves pattern matching and pre-grouping of the containers to be loaded based on a historical loading case database, including: Extract the container combination characteristics of each historical loading record from the historical loading case database, including weight distribution, destination combination, and loading and unloading sequence; Using the physical attribute data, destination data, and loading / unloading time window data of the containers to be loaded, a feature vector of the current batch of containers is constructed. Calculate the cosine similarity between the container feature vector of the current batch and the container combination features of each historical loading record, and use the container grouping method corresponding to the historical loading record with a cosine similarity greater than a preset similarity threshold as a candidate grouping scheme. The candidate grouping schemes are deduplicated and merged to generate a set of candidate groups with a matching degree higher than a preset threshold.

3. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 2, characterized in that, The second stage employs a tabu search algorithm to allocate space to containers within a group, including: Based on the cargo area corresponding to each group in the candidate group set, an initial cargo space allocation matrix is ​​established. The rows of the initial cargo space allocation matrix represent the cargo space column number, the columns represent the cargo space layer number, and the elements represent the container numbers of the corresponding cargo space allocation. The neighborhood operation is defined as swapping the slot positions of different containers within two groups, or moving a single container to an available slot within the same group. Using ship stability and strength constraints as hard constraints, and minimizing the maximum stacking height and the number of crate flipping operations as optimization objectives, the neighborhood operation is performed iteratively, and the solution and its update count are recorded after each iteration. Set up a tabu table to store solutions that have been visited in the most recent iterations to avoid circular searches; When the preset maximum number of iterations is reached or the target value is not improved after several consecutive iterations, the current optimal solution is output as the local optimal loading scheme.

4. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 3, characterized in that, In the third stage, the locally optimal loading scheme is used as the initial population, and a multi-objective global optimization is performed using a genetic algorithm based on non-dominated sorting, including: The chromosome corresponding to the local optimal loading scheme is encoded as a loading sequence. Each chromosome consists of multiple gene segments, and each gene segment corresponds to a cargo space and its assigned container number. Multiple initial individuals, including the locally optimal loading scheme, are randomly generated to form an initial population. The initial population is sorted non-dominated, and the crowding distance of each individual is calculated. The next generation population is generated through tournament selection, single-point crossover, and site mutation. The crossover operation swaps the container allocations of corresponding compartment columns in the two parent chromosomes, and the mutation operation randomly changes the container number of a certain compartment to a valid replacement value. Repeat the non-dominated sorting, selection, crossover, and mutation steps until the maximum number of generations is reached; From the Pareto optimal solution set of the last generation, the comprehensive fitness of each solution is calculated according to the preset weight vector, and the solution with the highest comprehensive fitness is selected as the final loading scheme.

5. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 1, characterized in that, The physical attribute data includes container weight, size, and type; the destination data includes the order of unloading ports; and the loading and unloading time window data includes the estimated arrival time and estimated departure time.

6. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 1, characterized in that, The ship stability constraints include heel angle limits, longitudinal angle limits, and lower limit of initial metacentric height, while the strength constraints include allowable deck loads, allowable bilge loads, and hull bending moment limits.

7. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 1, characterized in that, The multi-stage optimization algorithm dynamically incorporates real-time port operation data during execution, including the progress of quay crane operations and the real-time location of containers in the yard. When the actual loading and unloading sequence deviates from the planned sequence in the final loading scheme by more than a preset deviation threshold, the genetic algorithm in the third stage is triggered to re-optimize and update the final loading scheme.

8. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 1, characterized in that, In the genetic algorithm based on non-dominated sorting, adaptive crossover probability and adaptive mutation probability are used. The adaptive crossover probability is dynamically adjusted according to the non-dominated level distribution of the current population, and the adaptive mutation probability is dynamically adjusted according to the mutation coefficient of the crowding distance of individuals.

9. The container loading and stowage method based on a multi-stage optimization algorithm as described in claim 1, characterized in that, The tabu search algorithm uses a hash table to store the summary information of visited solutions. The summary information includes the hash value of the space allocation matrix and the objective function value. The tabu object is a complete loading scheme, and the tabu length is set to a dynamic value that is positively correlated with the total number of containers to be loaded.

10. A container loading and stowage system based on a multi-stage optimization algorithm, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the steps of the container loading and stowage method based on a multi-stage optimization algorithm as described in any one of claims 1 to 9.

Images