A split alternating optimization thrust distribution method for a full-rotation ship

Abstract

The application discloses a kind of full-rotation ship decomposed type alternate optimization thrust distribution method, belongs to ship propeller control technical field.This method will full-rotation ship thrust distribution problem be divided into thrust amplitude and azimuth alternate optimization sub-problems, by constructing the augmented Lagrangian function with box constraint, in combination with gradient descent and backtracking line search strategy, give the explicit solution of optimal thrust and azimuth, realize propeller thrust and azimuth collaborative optimization.This method can unify the processing singularity of geometric matrix, thrust and azimuth size, change rate, energy consumption, accuracy and other constraint requirements, with low computational complexity, high efficiency, strong real-time and numerical stability, avoid the mutation phenomenon of traditional algorithm under boundary constraint, suitable for the high-frequency control demand of intelligent tugboat and other full-rotation propelling ships, with good engineering adaptability and popularization value.

Description

Technical Field

[0001] This invention belongs to the field of ship propulsion thrust distribution technology, and relates to a thrust distribution method for azimuth-rotating ships, specifically a decomposed alternating optimized thrust distribution method for azimuth-rotating ships. Background Technology

[0002] With the rapid development of the marine industry and port shipping, the number of vessels equipped with azimuth thrusters has increased rapidly. Typical examples include intelligent port tugboats, marine engineering support vessels, large offshore platform auxiliary vessels, and engineering barges. These vessels, operating in complex environments such as berthing and unberthing assistance, towing and escorting large vessels, precise positioning operations, and emergency rescue, need to perform tasks such as low-speed precise maneuvering, lateral berthing, large-angle turns, and attitude fine-tuning. This places higher demands on the coordinated control of the propulsion system. In these multi-thrust azimuth thrusters, the thrust distribution unit's task is to rationally distribute the desired resultant force and torque of the three degrees of freedom given by the upper controller to each azimuth thruster. This requires satisfying physical constraints such as thrust amplitude, azimuth angle, and rate of change, while also considering energy consumption, thruster wear, and maneuverability. It is a crucial link in achieving safe and efficient operation of intelligent tugboats and other azimuth thrusters.

[0003] Currently, thrust allocation methods mainly include pseudo-inverse algorithms, grouping decomposition methods, swarm intelligence algorithms, and sequential quadratic programming (SQP) algorithms.

[0004] Regarding the aforementioned pseudo-inverse algorithm, Chinese patent application number 2023103086709 discloses a control and thrust distribution method for a straight-wing propeller. This technology mainly targets the constant-speed pitch control method for straight-wing propellers, establishing a control command modulation method based on a cubic fitting model of thrust coefficient and eccentricity. At the thrust distribution level, a nonlinear optimization problem with minimum power consumption as the objective function is constructed, and a solution method based on KT conditions is adopted. Furthermore, for the quadratic programming subproblem in the solution process, an improved SQP algorithm using a pseudo-inverse method is proposed to determine the initial optimization value, solving the problem of efficient control of straight-wing propellers in dynamic positioning systems, achieving the effects of improving thrust output efficiency, reducing energy consumption, and reducing wear and noise. However, this technology has the following drawbacks: the method is highly dependent on the thrust coefficient-eccentricity model of a specific propeller, resulting in poor versatility; although the initial value provided by the pseudo-inverse method improves the convergence of SQP, it is still prone to getting trapped in local extrema under complex nonlinear constraints; the algorithm as a whole relies on multiple matrix decompositions and active set switching, resulting in high computational complexity and difficulty in meeting the requirements of high-frequency real-time control.

[0005] Regarding the aforementioned grouping and decomposition method, Chinese patent application number 2012101774759 discloses a thrust allocation method for dynamically positioned vessels based on thrust allocation management. This technology mainly employs a thrust allocation management module, which groups the thrusters, sets offset values, and linearizes the maximum thrust constraints of the omnidirectional thrusters. It then applies quadratic programming, particularly the smooth Newton method, to optimize the thrust allocation solution, thus solving the thrust allocation problem for multi-thruster systems. This achieves the effects of improving the thruster system's execution capability, ensuring the controller's thrust requirements, making thruster use more rational, and simplifying operation. However, this technology has the following drawbacks: linearization approximation introduces modeling errors, affecting thrust allocation accuracy; the iterative process of the smooth Newton method relies on multiple Jacobian matrix decompositions and updates, resulting in high computational overhead and poor real-time performance; furthermore, this method fails to uniformly consider thruster azimuth rate constraints and geometric singularities, making it prone to solution instability or failure under complex operating conditions.

[0006] Regarding the aforementioned swarm intelligence algorithm, Chinese patent application No. 2023100833175 discloses a ship thrust allocation method based on particle swarm optimization (PSO). This technology mainly employs PSO combined with an optimal azimuth angle zone delineation mechanism and a switching flag to dynamically select either PSO or SQP for thrust allocation, thereby optimizing indicators such as energy consumption, deviation, and wear. This solves the problems of thrust allocation accuracy and global optimization, achieving improved thrust allocation accuracy and energy conservation and emission reduction. However, this technology has the following drawbacks: the PSO algorithm has a slow convergence speed and strong randomness, and the stability of the results is significantly affected by parameter adjustments; the switching logic between PSO and SQP is complex, and the algorithm structure is not fixed, leading to uncertain periodic responses in the control system; furthermore, swarm intelligence algorithms require extensive iterative searches, resulting in long computation times and making real-time operation difficult in embedded systems.

[0007] Regarding the aforementioned algorithm combining SQP and swarm intelligence, Chinese patent application number 2018106861356 discloses a thrust allocation method based on a combination of particle swarm optimization and sequential quadratic programming. This technique primarily uses PSO for global optimization, then uses the result as the initial value for the SQP algorithm for local fine-tuning to solve the thrust allocation problem in dynamic positioning systems, achieving the effects of improving the algorithm's global optimization performance and reducing thruster wear and energy consumption. However, this technique has the following drawbacks: while combining the two algorithms improves optimization performance, it significantly increases computational complexity and runtime; the convergence characteristics of the PSO part depend on random factors, and the output results may fluctuate significantly under different operating conditions; furthermore, the local optimization of SQP still requires repeated matrix factorization, making it difficult to guarantee the real-time performance and stability of the overall method.

[0008] To address the aforementioned swarm intelligence algorithms, Chinese patent application No. 2020114780930 discloses a ship thrust allocation method based on a mixed swarm optimization algorithm. This technology primarily employs an innovative "maximum pursuit" mixed swarm optimization algorithm, combining a bidirectional optimization strategy, performing local optimization on both the best and worst individuals, and introducing a sector combination selection mechanism based on dwell time and lag switching technology. This solves the global optimization problem in thrust allocation, achieving the effect of reducing ship system energy consumption while effectively reducing propeller wear and avoiding thrust singular structures. However, this technology has the following drawbacks: the search logic of the mixed swarm optimization algorithm is complex, the population evolution and parameter adjustment processes involve a large amount of computation, making it difficult to guarantee the real-time response of the algorithm in high-frequency control; at the same time, the multi-threshold switching mechanism is sensitive to the dynamic performance of the system, and the engineering implementation and debugging process is complex, limiting its real-time deployment and hardware portability.

[0009] In summary, while the grouping decomposition method and pseudo-inverse algorithm are computationally fast, they are difficult to comprehensively consider the actual constraints of the thruster, such as mechanical limiting, rate constraints, and power limits. The SQP algorithm can handle constrained problems, but its optimization process relies on multiple matrix decompositions and active constraint switching, making the computational steps cumbersome and prone to numerical instability. Swarm intelligence algorithms, such as particle swarm optimization (PSO), genetic algorithms, and dung beetle algorithms, have global optimization capabilities, but their convergence speed is slow, parameters are sensitive, and computational structures are not fixed, making it difficult to meet the computational cycle requirements of real-time control.

[0010] Currently, with the widespread application of azimuth thrusters on intelligent tugboats and other vessels, the nonlinearity, coupling, and multi-constraint characteristics of thrust distribution problems are becoming increasingly prominent. Traditional methods suffer from a series of common problems, including high computational complexity, unstable structure, sensitivity to parameters and initial values, poor real-time performance, and inconsistent constraint handling. These methods struggle to simultaneously meet the requirements of computational speed, numerical stability, and solution feasibility, and cannot satisfy the stringent requirements of high-frequency real-time control and thrust smoothness for multi-thrust azimuth ships under complex operating conditions. Therefore, it is necessary to design a new thrust distribution method for azimuth ships. Summary of the Invention

[0011] The purpose of this invention is to address a series of problems existing in the thrust allocation methods for multi-thrust azimuth ships in practical engineering, such as high computational complexity, insufficient real-time performance, scattered constraint processing, and unstable solutions. This invention proposes a decompositional alternating optimization thrust allocation method for azimuth ships. By decomposing and co-optimizing thrust and azimuth, this method can simultaneously meet multiple constraints such as thrust, azimuth and rate of change, energy consumption, and accuracy. It significantly reduces computational complexity, improves algorithm real-time performance, ensures the smoothness and convergence stability of azimuth adjustment, and avoids abrupt changes at constraint boundaries that occur in traditional methods. This method is suitable for real-time, stable, and engineering-feasible optimal thrust allocation for azimuth ships such as intelligent tugboats.

[0012] This application provides a decomposed alternating optimized thrust distribution method for azimuth-revolving ships, which adopts the following technical solution:

[0013] A decomposed alternating optimized thrust distribution method for azimuth-based ships, characterized by the following steps:

[0014] S1: Establish a thrust distribution model for a multi-propeller azimuth ship. The model includes the thrust of each propeller, the resultant force and resultant moment vector, the propeller geometry matrix, and the thrust amplitude vector of each propeller.

[0015] S2: Under the condition of simultaneously considering the singularity of the geometric matrix, the magnitude and rate of change of the thruster thrust and azimuth angle, energy consumption, accuracy and other actual physical constraints, construct an objective function that minimizes the thrust distribution error.

[0016] S3: Solve the thrust distribution problem using a decompositional alternating optimization method, including:

[0017] S3-1: Under the condition of fixed thruster azimuth angle, solve the thrust amplitude subproblem. The subproblem is a convex quadratic programming form with box constraints. It is solved by constructing an augmented Lagrangian function and performing an alternating optimization update step.

[0018] S3-2: With a fixed thrust amplitude, optimize the thruster azimuth angle by combining gradient descent and backtracking search methods to update the azimuth angle and perform constrained projection.

[0019] By adopting the above technical solution, a thrust distribution model for a multi-thrust, azimuth-rotor ship is established, an optimization objective function is constructed, and a decompositional alternating optimization method is applied to solve it. This method solves for the thrust amplitude when the thruster azimuth is fixed, and optimizes the thruster azimuth when the thrust amplitude is fixed. This solves the problems of high computational complexity, insufficient real-time performance, scattered constraint handling, and unstable solutions in existing thrust distribution methods, making it difficult to meet the requirements of high-frequency real-time control. It improves the computational efficiency of thrust distribution, meeting the requirements of high-frequency real-time control; it can uniformly handle multiple constraint conditions, ensuring the global feasibility and stability of the solution; and it reduces computational complexity, shortens computation time, and improves the real-time performance and reliability of the system.

[0020] Furthermore, the mathematical description of the thrust distribution model is as follows:

[0021] ;

[0022] in, The resultant force / moment vector of the ship; It is a geometric matrix; Here are the azimuth vectors for each thruster. This is the thrust amplitude vector.

[0023] By adopting the above technical solution, a mathematical description of the thrust distribution model is provided, clarifying the relationship between the ship's total resultant force / moment vector, geometric matrix, and thrust amplitude vector. This solves the problem that existing models lack a precise mathematical description of the thrust distribution problem, leading to a complex solution process and difficulty in optimization. It facilitates the understanding and implementation of the thrust distribution algorithm, providing a foundation for the subsequent construction of the objective function and constraints, ensuring the accuracy and operability of the model.

[0024] Furthermore, the optimization objective function described in step S2 is:

[0025] ;

[0026] in, For the desired control signal, This is the weight matrix. For regularization parameters, It is the identity matrix. This is the singularity penalty coefficient. It is a constant.

[0027] By adopting the above technical solution, a weight matrix, a singularity penalty coefficient, and a regularization term are introduced into the objective function to minimize thrust allocation error and consider geometric matrix singularity constraints. This solves the problem that existing optimization objective functions fail to fully consider the impact of geometric matrix singularity on thrust allocation, leading to instability in solutions under complex operating conditions. It improves the accuracy and stability of thrust allocation under geometric matrix singularity conditions, and the weight matrix helps balance the optimization objective, avoiding overfitting and underfitting problems.

[0028] Furthermore, the constraints described in step S2 include:

[0029] Thrust amplitude constraint: ,in , These represent the minimum and maximum thrust, respectively.

[0030] Azimuth size constraints: ,in , These are the minimum and maximum limits of the azimuth angle, respectively.

[0031] Thrust amplitude change rate constraint: ,in , These represent the minimum and maximum rates of thrust change, respectively. The magnitude of the thrust at the previous moment. The sampling period;

[0032] Azimuth rate of change constraint: ,in , These represent the minimum and maximum rates of change of azimuth, respectively. This represents the azimuth angle at the previous moment.

[0033] By adopting the above technical solution, the constraints in thrust distribution are clarified, including upper and lower limits for thrust amplitude, azimuth angle, and rate of change. This solves the problem in existing methods where constraints are scattered and treated inconsistently, leading to complex solutions and the tendency for sudden thrust changes or solution failures under boundary conditions. It can uniformly handle multiple constraints, ensuring the global feasibility and stability of the solution, avoiding sudden thrust changes and solution failures, and improving the numerical stability and control stability of the system.

[0034] Furthermore, the solution steps of the alternating optimization method described in step S3 include:

[0035] (1) Constructing the augmented Lagrange function

[0036] ;

[0037] in, , , , , For box constraint indicator functions, The penalty coefficient is... As an auxiliary variable, To scale the dual variable.

[0038] (2) Perform alternating optimization and update steps

[0039] (2-1) Main Variable renew: subscript Indicates the number of program iterations. This represents the total number of steps.

[0040] (2-2) Auxiliary variables renew: ,in For projection operators;

[0041] (2-3) Dual variables renew: .

[0042] (3) After executing the above three-step alternating optimization procedure, a stopping criterion is set. When all three conditions are met:

[0043] ;

[0044] ;

[0045] Terminate program iteration, where, Represents the original residual, used to measure the degree to which the consistency constraints of variables are satisfied; This represents the dual residual, used to measure the stability of auxiliary variable updates; and These are the convergence thresholds for the original and dual residuals, respectively. At this point, the optimal solution for the thrust amplitude is obtained. .

[0046] If the above convergence conditions are not met, the above three-step alternating optimization procedure will continue to be executed, and updates will be performed. Continue until the stopping criterion is met, and then output the optimal solution.

[0047] By adopting the above technical solution, the solution steps of the alternating optimization method are described in detail, including the construction of the augmented Lagrangian function, the alternating update steps, and the stopping criterion. This solves the problems of high computational complexity and poor real-time performance of existing solution methods, making it difficult to meet the requirements of high-frequency real-time control. It improves solution efficiency, reduces computation time, and meets the needs of high-frequency real-time control. It also enhances the numerical stability and convergence of the algorithm, ensuring reliability under complex operating conditions.

[0048] Furthermore, the optimization step of the thruster azimuth angle in step S3-2 includes:

[0049] (1) Calculate the gradient of the objective function

[0050] According to the objective function right Take the partial derivative to obtain the azimuth gradient. , is used to characterize the fastest direction of increase of the objective function;

[0051] (2) The azimuth angle is updated by combining gradient descent and backtracking search.

[0052] Set initial step size Retrospective factor and the decreasing control constant In each iteration, the trial azimuth angle is calculated first:

[0053] ;

[0054] in, , , This is the projection operator.

[0055] If the objective function value of the trial solution satisfies the descent condition:

[0056] ;

[0057] The program will then terminate and output the optimal thruster azimuth angle. ;

[0058] Otherwise, let Repeat the above trial until the descent condition is met, and output the optimal thruster azimuth angle.

[0059] By employing the aforementioned technical solution, the optimization steps for the thruster azimuth angle are described in detail, including the gradient calculation of the objective function, the update strategy combining gradient descent and backtracking search, and the application of the azimuth projection operator. This addresses the problem of existing methods lacking effective update strategies and constraint handling mechanisms when optimizing the thruster azimuth angle, leading to poor optimization results. It improves the accuracy and smoothness of azimuth angle optimization, avoids the thrust abrupt change problem in traditional methods, and ensures the stability and reliability of the optimization process through backtracking search and the azimuth projection operator.

[0060] Re-execute step S3-1 to solve for the thrust amplitude, and output the optimal thrust and azimuth results at the current sampling time: update the optimal azimuth from step S3-2. Substitute the thrust calculation process described in step S3-1, recalculate the thrust amplitude, and obtain the optimal thrust allocation result that matches the current azimuth angle.

[0061] Finally, the thrust amplitude and azimuth angle commands of each thruster are output as the execution signal of the current sampling period and sent to the propulsion system. They are also used as the initial values ​​for the next sampling period to participate in the subsequent decompositional alternating optimization thrust distribution calculation, so as to achieve continuous and smooth real-time control.

[0062] In summary, the present invention has the following beneficial technical effects:

[0063] (1) Structured decomposition to improve convergence stability. This invention proposes a decompositional alternating optimization algorithm structure of "thrust-azimuth angle decomposition solution" for thrust allocation in azimuth-rotating ships. Through structured decomposition, the original nonlinear coupled optimization problem is divided into a thrust amplitude subproblem and an azimuth angle subproblem. The thrust subproblem is solved using the alternating direction multiplier method, realizing convex quadratic programming optimization with a fixed matrix structure. The azimuth angle subproblem is solved using a strategy combining gradient descent and backtracking search, and an azimuth angle projection constraint and barrier term smoothing adjustment mechanism are introduced. This decompositional alternating optimization structure can ensure that each iteration is physically feasible and numerically stable, and converges quickly within a finite number of steps, making it particularly suitable for high-dimensional coupled scenarios of azimuth-rotating thrusters.

[0064] (2) The computational efficiency is greatly improved, meeting the requirements of high-frequency real-time control. This invention decomposes the problem into matrix multiplication and projection operations with a fixed structure through an alternating optimization framework, and the core matrix can be decomposed in advance. This avoids complex calculations in online iteration, realizes a fixed and low-complexity solution process, greatly improves the computational speed, and can stably meet the high-frequency real-time control requirements of 50-100Hz.

[0065] (3) The output is smooth and stable, which significantly improves the thrust change problem during constraint switching. In the process of solving the fixed thrust amplitude and azimuth angle alternately, the present invention introduces a dual variable smoothing adjustment mechanism. This mechanism enables the thruster to approach the physical boundary smoothly rather than suddenly hitting the boundary, avoiding the thrust change phenomenon caused by the traditional method during constraint switching, and improving the numerical stability and control stability of the system.

[0066] (4) Enhanced constraint handling capability, ensuring global feasibility and robustness of the solution. This invention employs a box-constrained projection operator to uniformly handle various physical constraints such as thrust amplitude, azimuth angle, and rotational speed. This method requires no manual intervention or judgment of the active set, automatically and rigorously constraining all solutions within the physically feasible range, achieving a 100% constraint satisfaction rate, and greatly enhancing the reliability and robustness of the algorithm under various working conditions.

[0067] (5) It has strong engineering practicality and is easy to deploy on embedded platforms. The computational steps required by this invention are constant and the memory usage is controllable, making it extremely suitable for implementation on embedded hardware with limited computing resources. This not only reduces system power consumption and hardware costs, but also provides technical support for the widespread application of high-performance dynamic positioning control algorithms in real ship control systems.

[0068] (6) Wide range of applications, suitable for various types of azimuth-rotor vessels such as intelligent tugboats. This invention does not depend on specific ship types or specific propeller model parameters, and can be widely used in various vessels equipped with azimuth-rotor propellers, such as intelligent port tugboats, engineering support vessels, and offshore platform operation vessels. Attached Figure Description

[0069] Figure 1 This is a schematic diagram of the method flow of the present invention.

[0070] Figure 2 This is a schematic diagram of the thruster arrangement in an embodiment of the present invention.

[0071] Figure 3 This is a schematic diagram of the thrust variation curve of the thruster 1 in working condition one of the embodiments of the present invention.

[0072] Figure 4 This is a schematic diagram of the thrust variation curve of the thruster 2 in working condition one of the embodiments of the present invention.

[0073] Figure 5 This is a schematic diagram of the azimuth angle change curve of the thruster 1 in working condition one of the embodiments of the present invention.

[0074] Figure 6 This is a schematic diagram of the azimuth angle change curve of the thruster 2 in working condition one of the embodiments of the present invention.

[0075] Figure 7 This is a schematic diagram of the combined thrust and torque variation curves in working condition one of the embodiments of the present invention.

[0076] Figure 8 This is a schematic diagram of the thrust change rate curve of the thruster 1 in working condition one of the embodiments of the present invention.

[0077] Figure 9 This is a schematic diagram of the thrust change rate curve of the thruster 2 in working condition one of the embodiments of the present invention.

[0078] Figure 10 This is a schematic diagram of the azimuth angle change rate curve of the thruster 1 in working condition one of the embodiments of the present invention.

[0079] Figure 11 This is a schematic diagram of the azimuth angle change rate curve of the thruster 2 in working condition one of the embodiments of the present invention.

[0080] Figure 12 This is a schematic diagram comparing the smoothness of three algorithm instructions in working condition one of the embodiments of the present invention.

[0081] Figure 13 This is a schematic diagram comparing the real-time performance of three algorithm instructions in working condition one of the embodiments of the present invention.

[0082] Figure 14 This is a schematic diagram of the thrust variation curve of the thruster 1 in working condition 2 of the present invention.

[0083] Figure 15 This is a schematic diagram of the thrust variation curve of the thruster 2 in working condition two of the present invention.

[0084] Figure 16 This is a schematic diagram of the azimuth angle change curve of thruster 1 in working condition 2 of the present invention.

[0085] Figure 17 This is a schematic diagram of the azimuth angle change curve of the thruster 2 in working condition two of the present invention.

[0086] Figure 18 This is a schematic diagram of the combined thrust and torque variation curves in working condition two of the present invention.

[0087] Figure 19 This is a schematic diagram of the tugboat trajectory change curve in working condition two of the present invention.

[0088] Figure 20 This is a schematic diagram of the thrust change rate curve of thruster 1 in working condition 2 of the present invention.

[0089] Figure 21 This is a schematic diagram of the thrust change rate curve of the thruster 2 in working condition 2 of the present invention.

[0090] Figure 22 This is a schematic diagram of the azimuth angle change rate curve of the thruster 1 in working condition 2 of the present invention.

[0091] Figure 23 This is a schematic diagram of the azimuth angle change rate curve of the thruster 2 in working condition two of the present invention.

[0092] Figure 24 This is a schematic diagram comparing the smoothness of three algorithm instructions in working condition two of the present invention.

[0093] Figure 25 This is a schematic diagram comparing the real-time performance of three algorithm instructions in working condition two of the present invention. Detailed Implementation

[0094] The present invention will be further described below with reference to the accompanying drawings and embodiments:

[0095] A decomposed alternating optimized thrust distribution method for azimuth-based ships, characterized by the following specific steps:

[0096] like Figure 2 As shown, establish the attached coordinate system. ,in For the ship's center of gravity, The axis points towards the bow. Pointing to the starboard side of the ship. In order to produce the desired motion behavior, the ship controller will generate the desired control signal. ,in / They are longitudinal / lateral thrust, This is the turning torque.

[0097] A thrust distribution model for a multi-propeller, azimuth-rotor ship is established, and its mathematical description is as follows:

[0098] Ship equipped with Taiwan full-rotation thruster, the first Taiwan thruster The position below is The magnitude of the thrust generated is azimuth angle is Then each thruster in The resultant force / torque vector generated below can be expressed as:

[0099] ;

[0100] in: For geometric matrices, Here are the azimuth vectors for each thruster. This represents the thrust amplitude vector for each thruster.

[0101] Considering physical constraints such as geometric matrix singularity, thruster thrust and azimuth angle magnitudes and rates of change, and energy consumption, an optimization objective function is constructed to minimize thrust distribution error:

[0102] ;

[0103] in, This is the weight matrix. For regularization parameters, It is the identity matrix. This is the singularity penalty coefficient. It is a constant.

[0104] Physical constraints such as the magnitude and rate of change of thruster thrust and azimuth angle include:

[0105]

[0106] in, , These represent the minimum and maximum thrust of the thruster, respectively. These are the thruster azimuth limits, , , and These are the limits for the rate of change of thrust and azimuth angle, respectively. and The magnitudes of the thrust and azimuth at the previous moment. The sampling period.

[0107] The thrust distribution problem is solved using a decompositional alternating optimization method. First, with the thruster azimuth angle fixed, the thrust amplitude subproblem is solved. This subproblem can be expressed as a convex quadratic programming problem with box constraints:

[0108] ;

[0109] in, , , , .

[0110] Solving this convex quadratic programming problem involves the following steps:

[0111] First, to facilitate handling box constraints, auxiliary variables are introduced. and impose consistency constraints Then the problem is equivalent to

[0112] ;

[0113] in, Let be the box constraint indicator function, To facilitate numerical implementation, we introduce a scaling dual variable for the unscaled Lagrange multipliers. , Then the augmented Lagrange function can be written as:

[0114] ;

[0115] in, This is the penalty coefficient.

[0116] Then, based on the augmented Lagrangian function, a three-step alternating optimization procedure is executed to update sequentially. .

[0117] The first step is to fix , Minimize ,right Taking the partial derivative and setting it to zero, we get

[0118] ;

[0119] main variables Updated to:

[0120] ;

[0121] Among them, subscript Indicates the number of program iterations. This represents the total number of iterations. This step is equivalent to solving a set of linear equations, and its coefficient matrix... The matrix remains constant during system operation. During the initialization phase, a Cholesky decomposition can be performed on the matrix and its inverse value can be stored. In subsequent solution processes, the explicit thrust solution can be obtained by substitution, which greatly reduces the amount of computation and significantly improves the real-time performance of the computation.

[0122] The second step is to fix it. and In this case, update the auxiliary variable. :

[0123] ;

[0124] in, The projection operator restricts the optimization results to the feasible region.

[0125] The third step is to update the dual variables. :

[0126] ;

[0127] This step gradually strengthens consistency constraints, making and They converge to the same optimal point within the feasible region.

[0128] After performing the above three steps, a stopping criterion is set. The criterion is set when all of the following conditions are met:

[0129] ;

[0130] ;

[0131] Terminate program iteration, where, Represents the original residual, used to measure the degree to which the consistency constraints of variables are satisfied; This represents the dual residual, used to measure the stability of auxiliary variable updates; and These are the convergence thresholds for the original and dual residuals, respectively. At convergence, we have... Therefore, the optimal solution for the thrust amplitude is obtained. .

[0132] If the above convergence conditions are not met, the above three-step alternating optimization procedure will continue to be executed, and updates will be performed. Continue until the stopping criterion is met, and then output the optimal solution.

[0133] After terminating thrust amplitude optimization, the thruster azimuth angle is optimized. Maintain... Without changing the initial position, an azimuth angle update strategy combining gradient descent and backtracking search is adopted to solve for and constrain the projection of the azimuth angles of each thruster. Specifically, the following steps are included:

[0134] First, calculate the azimuth gradient: based on the objective function. right Take the partial derivative to obtain the azimuth gradient. , is used to characterize the fastest direction of increase of the objective function, and its specific expression is:

[0135] ;

[0136] The gradient of the singular term can be approximated using numerical difference.

[0137] Then, perform a backtracking search and update the azimuth angle: set the initial step size. Retrospective factor and the decreasing control constant In each iteration, the trial azimuth angle is first calculated:

[0138] ;

[0139] in, , , This is the projection operator.

[0140] If the objective function value of the trial solution satisfies the descent condition:

[0141] ;

[0142] Then, terminate the program and output the optimal thruster azimuth angle. .

[0143] Otherwise, let Repeat the above trial until the descent condition is met, and output the optimal thruster azimuth angle.

[0144] Re-execute step S3-1 to solve for the thrust amplitude, and output the optimal thrust and azimuth results at the current sampling time: update the optimal azimuth from step S3-2. Substitute the thrust calculation process described in step S3-1, recalculate the thrust amplitude, and obtain the optimal thrust allocation result that matches the current azimuth angle.

[0145] Finally, the thrust amplitude and azimuth angle commands of each thruster are output as the execution signal of the current sampling period and sent to the propulsion system. They are also used as the initial values ​​for the next sampling period to participate in the subsequent decompositional alternating optimization thrust distribution calculation, so as to achieve continuous and smooth real-time control.

[0146] Example

[0147] To verify the effectiveness and superiority of the thrust distribution method for azimuth-rotor vessels based on the decompositional alternating optimization algorithm, a scaled-down model of the "Yungang Electric Tug No. 1" intelligent tugboat equipped with two azimuth thrusters was used as the research object to conduct a comprehensive simulation study.

[0148] In this embodiment, the tugboat is equipped with two azimuth thrusters at the stern, such as... Figure 2 As shown. The position coordinates of the two thrusters are (-1.2, 0.8) and (-1.2, -0.8) respectively; the initial state of the thrusters. , The parameters of the objective function are: , , Constraints , , , , , Sampling period Initial values ​​in thrust iteration , , Regularization parameters Penalty coefficient Stop condition threshold , Initial step size in azimuth update related parameters Retrospective factor .

[0149] To verify the effectiveness of the method, two working conditions were designed to verify the performance of the thrust allocation algorithm.

[0150] Operating Condition 1: In practical applications, the generalized force instruction is desired. Typically generated in real-time by upper-level navigation, path planning, and attitude control modules, its values ​​vary with the task, environment, and control objectives, and can be considered as bounded time-varying inputs. To verify the applicability and robustness of the method of this invention under such general inputs, a set of continuous time-varying commands containing different frequencies and phases is constructed as representative samples. The input signal, under the condition of simultaneous excitation of three channels, can form various amplitude combinations and variation trends, thus enabling its use in simulation. The arbitrary time-varying values ​​that may occur in engineering projects place high demands on the allocatability, constraint handling capability, and numerical stability of the thrust allocation algorithm.

[0151] Working Condition 2: To verify the effectiveness of the algorithm in practical engineering applications, the desired control input designed in reference (DIO No. 10.1016 / j.oceaneng.2017.09.062) was used. .

[0152] In addition, to demonstrate the superiority of the decompositional alternating optimization thrust allocation algorithm of the present invention, the method of the present invention was compared with the particle swarm optimization (PSO) algorithm and the sequential quadratic programming (SQP) algorithm, respectively, and comprehensive simulation performance tests were carried out.

[0153] Table 1. Simulation data for three thrust distribution algorithms under operating condition 1

[0154] algorithm RMSE mean Mean IAE Total time elapsed (SECs(s)) PSO 0.169 1.740 6.355 SQP 0.125 1.357 3.398 Algorithm of this invention 0.092 0.410 0.320

[0155] Table 2. Simulation data for three thrust distribution algorithms under working condition 2.

[0156] algorithm RMSE mean Mean IAE Total time elapsed (SECs(s)) PSO 0.193 0.729 5.487 SQP 0.197 1.164 2.730 Algorithm of this invention 0.192 0.510 0.350

[0157] Simulation results are shown in Tables 1 and 2. Figures 3 to 25As shown. The command output of the dual full-spinning thrusters is given by the thrust curve and azimuth curve in the two operating conditions, respectively. The thrust and azimuth time histories for operating conditions one and two are shown below. Figures 3-6 and Figures 14-17 As shown in the figures, both sets of results indicate that PSO and SQP exhibit more pronounced abrupt changes and chattering during certain periods, and the command sequences are more prone to high-frequency adjustments. In contrast, the thrust and azimuth sequences of the algorithm presented in this invention are more continuous and smooth overall, which is beneficial for improving the stability of tugboat attitude control and the lifespan of the actuator.

[0158] The resultant force and moment tracking results for the two sets of working conditions are as follows: Figure 7 and Figure 18 As shown, the three algorithms can generally track both longitudinal and lateral force channels, but there are significant differences in transition moments and regions where allotropy deteriorates. PSO and SQP are more prone to local spikes or short-term deviations, with particularly pronounced fluctuations in the yaw moment channel. In contrast, the curve of the decompositional alternating optimization algorithm is closer to the reference overall and has smaller fluctuations, indicating that it has more stable tracking performance under non-convex constraint switching and geometric condition deterioration. The statistical results of RMSE and IAE in Tables 1 and 2 further quantify this advantage. The root mean square error and absolute integral error of the decompositional alternating optimization algorithm are the smallest in both sets of working conditions, indicating that the algorithm of this invention has higher resultant force tracking accuracy under typical working conditions and the most unfavorable rotating load conditions.

[0159] To verify that the rate of change satisfies the physical constraints of the thruster, the simulation results of the rate of change for conditions one and two are shown below. Figures 8-11 and Figures 20-23 The black dashed lines in the figure represent the upper and lower limits of the rate of change. Under both operating conditions, none of the three algorithms exhibited rate of change exceeding the limits. However, PSO and SQP approached the limits more frequently, showing edge-running characteristics in some sections. The algorithm of this invention, on the other hand, leaves a larger margin overall, indicating that it can more effectively suppress unnecessary rapid adjustments while meeting tracking accuracy requirements, thereby improving instruction smoothness and actuator friendliness. Furthermore, the instruction smoothness statistics for operating conditions one and two are as follows: Figure 12 and Figure 24 As shown, the decompositional alternating optimization algorithm addresses the total thrust variation. and total azimuth variation Both metrics are smaller, which corroborates the phenomenon of "more continuous and less chattering" in the thrust and azimuth time history.

[0160] At the motion response level, the planar trajectory results for condition two are as follows: Figure 19 As shown, the trajectory generated by the decompositional alternating optimization algorithm has a higher degree of fit and smaller deviation from the expected trajectory compared to the SQP algorithm. This indicates that its thrust distribution output not only has a lower error at the resultant force and moment level, but also brings a more stable trajectory response performance at the kinematic level.

[0161] Regarding real-time performance, as shown in Tables 1 and 2, the total time consumption of the decompositional alternating optimization algorithm in the two sets of working conditions is 0.320 s and 0.350 s, respectively, significantly lower than that of PSO (6.355 s and 5.487 s) and SQP (3.398 s and 2.730 s). In working condition one, the total time consumption of the algorithm of this invention is reduced by approximately 94.96% and 90.58% compared to PSO and SQP, respectively, and the calculation speed is increased by approximately 19.8 times and 10.6 times, respectively. In working condition two, the total time consumption of the decompositional alternating optimization algorithm is reduced by approximately 93.62% and 87.17% compared to PSO and SQP, respectively, and the calculation speed is increased by approximately 15.7 times and 7.8 times, respectively. Further conversion based on the average calculation time per step, the decompositional alternating optimization algorithm maintains sub-millisecond levels in both sets of working conditions, while PSO and SQP are in the millisecond level, indicating that the algorithm of this invention can run stably at higher refresh rates. Combining the RMSE and IAE statistics in Tables 1 and 2, as well as the thrust, azimuth, and rate of change curves in the figure, it can be seen that the algorithm of this invention can significantly reduce the computational cost while still maintaining a smaller or equivalent tracking error and outputting a smoother thruster command sequence.

[0162] In summary, under the same amplitude and rate of change constraints, the decompositional alternating optimization thrust allocation algorithm of this invention maintains a small resultant force tracking error in both sets of operating conditions and outputs a more continuous and smoother thrust and azimuth command sequence. Compared with PSO and SQP, the algorithm of this invention significantly reduces computation time while ensuring tracking accuracy and smoothness, demonstrating stronger real-time online solution capabilities. Therefore, the algorithm of this invention achieves a better balance between accuracy, smoothness, and computational efficiency, and can meet the real-time thrust allocation requirements of intelligent tugboats and other azimuth-driven vessels under complex operating conditions.

Claims

1. A decomposed alternating optimized thrust distribution method for azimuth-based ships, characterized in that, Includes the following steps: S1: Establish a thrust distribution model for a multi-propeller azimuth ship. The model includes the thrust of each propeller, the resultant force and resultant moment vector, the propeller geometry matrix, and the thrust amplitude vector of each propeller. S2: Under the condition of simultaneously considering the singularity of the geometric matrix, the magnitude and rate of change of the thruster thrust and azimuth angle, energy consumption, accuracy and other actual physical constraints, construct an objective function that minimizes the thrust distribution error. S3: Solve the thrust distribution problem using a decompositional alternating optimization method, including: S3-1: Under the condition of fixed thruster azimuth angle, solve the thrust amplitude subproblem. The subproblem is a convex quadratic programming form with box constraints. It is solved by constructing an augmented Lagrangian function and performing an alternating optimization update step. S3-2: With a fixed thrust amplitude, optimize the thruster azimuth angle by combining gradient descent and backtracking search methods to update the azimuth angle and perform constrained projection.

2. The method for decomposed alternating optimized thrust distribution of a fully azimuth-rotating ship according to claim 1, characterized in that, The mathematical description of the thrust distribution model in step S1 is as follows: ； in, The resultant force / moment vector of the ship; It is a geometric matrix; Here are the azimuth vectors for each thruster. This represents the thrust amplitude vector for each thruster.

3. The method for decomposed alternating optimized thrust distribution for a fully azimuth-rotating ship according to claim 1, characterized in that, The optimization objective function mentioned in step S2 is: ； in, For the desired control signal, This is the weight matrix. For regularization parameters, It is the identity matrix. This is the singularity penalty coefficient. It is a constant.

4. The method for decomposed alternating optimized thrust distribution for a fully azimuth-rotating ship according to claim 1, characterized in that, The constraints mentioned in step S2 include: Thrust amplitude constraint: ,in , These represent the minimum and maximum thrust, respectively. Azimuth size constraints: ,in , These are the minimum and maximum limits of the azimuth angle, respectively. Thrust amplitude change rate constraint: ,in , These represent the minimum and maximum rates of thrust change, respectively. The magnitude of the thrust at the previous moment. The sampling period; Azimuth rate of change constraint: ,in , These represent the minimum and maximum rates of change of azimuth, respectively. This represents the azimuth angle at the previous moment.

5. The method for decomposed alternating optimized thrust distribution for a fully azimuth-rotating ship according to claim 1, characterized in that, The solution steps of the decompositional alternating optimization method described in step S3 include: (1) Constructing the augmented Lagrange function ； in, , , , , For box constraint indicator functions, The penalty coefficient is... As an auxiliary variable, To scale the dual variable; (2) Perform alternating optimization and update steps (2-1) Main Variable renew: subscript Indicates the number of program iterations. This represents the total number of steps. (2-2) Auxiliary variables renew: ,in For projection operators; (2-3) Dual variables renew: ; (3) After executing the above three-step alternating optimization procedure, a stopping criterion is set. When all three conditions are met: ； ； Terminate program iteration, where, Represents the original residual, used to measure the degree to which the consistency constraints of variables are satisfied; This represents the dual residual, used to measure the stability of auxiliary variable updates; and These are the convergence thresholds for the original and dual residuals, respectively; at this point, the optimal solution for the thrust amplitude is obtained. ; If the above convergence conditions are not met, the above three-step alternating optimization procedure will continue to be executed, and updates will be performed. Continue until the stopping criterion is met, and then output the optimal solution.

6. The method for decomposed alternating optimized thrust distribution of a fully azimuth-rotating ship according to claim 1, characterized in that, The optimization steps for the thruster azimuth angle in step S3-2 include: (1) Calculate the gradient of the objective function According to the objective function right Take the partial derivative to obtain the azimuth gradient. , is used to characterize the fastest direction of increase of the objective function; (2) The azimuth angle is updated by combining gradient descent and backtracking search. Set initial step size Retrospective factor and the decreasing control constant In each iteration, the trial azimuth angle is calculated first: ； in, , , For projection operators; If the objective function value of the trial solution satisfies the descent condition: ； The program will then terminate and output the optimal thruster azimuth angle. ; Otherwise, let Repeat the above trial until the descent condition is met, and output the optimal thruster azimuth angle.

7. The method for decomposed alternating optimized thrust distribution for a fully azimuth-rotating ship according to claim 1, characterized in that, After updating in step S3-2, the optimal azimuth angle will be... Substitute the thrust calculation process described in step S3-1, recalculate the thrust amplitude, and obtain the optimal thrust allocation result that matches the current azimuth angle.

Images