Improved parrot algorithm-based thrust distribution method under failure mode of ship propeller

Abstract

The invention provides a thrust distribution method under a ship propeller fault mode based on an improved parrot algorithm, and belongs to the technical field of dynamic positioning ship thrust distribution. According to the technical scheme, the method comprises the following steps that 1, various parameters of a ship are initialized; step 2, constructing a thruster equality constraint term; 3, initializing parameters of the improved parrot optimization algorithm; 4, performing individual variable updating by applying an improved parrot optimization algorithm; 5, judging whether the maximum number of iterations is reached or not, if so, executing the step 6, otherwise, executing the step 4; and 6, outputting a thrust distribution result, updating ship state parameters, and skipping to the step 2 to carry out thrust distribution calculation of the next cycle. The method has the beneficial effects that the method has good algorithm optimization ability in both the failure mode and the non-failure mode of the ship propeller, the thrust optimal solution can be quickly and reliably obtained, and the ship thrust distribution performance is improved.

Description

Technical Field

[0001] The invention relates to the technical field of thrust distribution for dynamic positioning ships, and in particular to a thrust distribution method under a ship propeller failure mode based on an improved Parrot algorithm. Background Art

[0002] As a high-tech and advanced ship operating system, the dynamic positioning system has been widely used in various types of ships, such as dredgers, drilling platforms, shuttle tankers, scientific research ships, pipeline and cable laying ships, etc. This type of ship uses thrusters to resist the external environmental forces acting on the ship, so that the ship can maintain a set position on the horizontal plane. It is usually equipped with multiple thrusters. For the desired thrust and torque, there are multiple thrust distribution combinations. How to select the thrust distribution combination method is the key to improving the thrust distribution performance. In addition, faults may occur during the use of the thruster. If these faults are not resolved in time, they will affect the ship's positioning performance and even cause serious accidents. Therefore, considering the thruster failure situation when studying the thrust distribution optimization algorithm has important practical application significance and value.

[0003] Traditional thrust allocation methods, such as sequential quadratic programming and pseudo-inverse method, have the advantages of small computational complexity and good real-time performance, but they have problems such as poor flexibility, easy to fall into local extreme points, and weak adaptability under complex working conditions. The swarm intelligence optimization algorithm developed in recent decades provides another path for solving the thrust allocation optimization problem because it does not depend on the cost function and constraint form of the object problem. More and more researchers have begun to use this type of optimization method to solve the thrust allocation optimization problem. However, most of the current thrust allocation method studies based on swarm intelligence optimization algorithms do not consider the thruster failure mode and stable and reliable convergence under a small number of iterations, which affects the safety and reliability of ship operation and is not conducive to the further improvement of dynamic positioning performance. Summary of the invention

[0004] The purpose of the present invention is to provide a thrust distribution method under a ship propeller fault mode based on an improved Parrot algorithm. The present invention constructs different cost functions and equality constraints for the normal working mode and different fault modes of the propeller, respectively, selects the Parrot optimization algorithm with the characteristics of behavioral diversity, high performance, and strong parameter adaptability as the thrust distribution optimization algorithm, and performs multi-strategy improvements on the Parrot optimization algorithm, so that it can obtain a stable and reliable thrust distribution optimization solution under the condition of a smaller number of iterations, thereby improving the dynamic positioning performance and the safety and reliability of ship operation.

[0005] In order to achieve the above-mentioned invention object, the technical solution adopted by the present invention is specifically: a thrust distribution method under a ship propeller fault mode based on an improved Parrot algorithm, comprising the following steps:

[0006] Step 1, initialize various ship parameters: propeller position parameters, initial state values, physical limit values ​​(including propeller thrust and angle limit values, thrust change rate and angle change rate limit values), propeller power and thrust relationship coefficient;

[0007] Step 2: Determine whether the current thruster has a fault and the fault type (the faults mainly include thrust faults and rotation faults), construct thruster equation constraint items (the thrust of each thruster adopts the extended thrust mode, and when the thruster angle is fixed or the rotation fault is stuck, the x-axis thrust or y-axis thrust of the thruster is adopted), and calculate the specific solution and general solution of the thruster equation constraint items (the general solution coefficient is selected as an individual variable);

[0008] Step 3: According to the current working condition of the thruster, initialize the parameters of the improved Parrot optimization algorithm: population size N, maximum number of iterations t max , individual dimension, population upper and lower boundaries, time threshold t mf and t f (t mf <t f <t max ), initialize the population, calculate the population fitness value and sort it in ascending order, obtain the optimal individual fitness value and its position, etc.;

[0009] Step 4: Use the improved Parrot optimization algorithm to update individual variables. The specific implementation steps are as follows:

[0010] Step 4.1, select the individual variables to be updated in turn, and select the individual update method by randomly generating four integers 1, 2, 3 and 4. When the random number is 1, execute step 4.2, when the random number is 2, execute step 4.3, when the random number is 3, execute step 4.4, and when the random number is 4, execute step 4.5;

[0011] Step 4.2: Use the parrot's foraging behavior to update the individual position. The specific calculation is as follows:

[0012]

[0013] Among them, X(t) is the current position of the individual, X best represents the current best individual position, Levy (D) is used to represent the parrot's Levy flight behavior, D refers to the individual dimension, P mean indicates the position of the entire group, t is the number of iterations, rand(0,1) is a random number between 0 and 1, and X mean (t) is the average value of the current position of the entire group. After the individual is updated, it is determined whether the N individual variables are updated. If yes, execute step 4.6, otherwise execute step 4.1;

[0014] Step 4.3: Use the parrot's stay behavior to update the individual position. The specific calculation is as follows:

[0015] X(t+1)=X(t)+X best Levy(D)+rand(0,1)ones(1,D),

[0016] Among them, rand(0,1)ones(1,D) indicates the position of random stay on the owner, and ones(1,D) means that the D-dimensional data is generated as a 1 vector. After the individual is updated, it is determined whether the N individual variables are updated. If yes, execute step 4.6, otherwise execute step 4.1;

[0017] Step 4.4: Use the parrot's communication behavior to update the individual position. The specific calculation is as follows:

[0018] When the random generation probability p≤0.5, the update method is: X(t+1)=0.2rand(0,1)(1-t / t max )(X(t)-X mean (t)),

[0019] On the contrary, when p>0.5, the update method is: X(t+1)=0.2rand(0,1)e^(-t / (rand(0,1)t max )),

[0020] After the individual is updated, determine whether the N individual variables are updated. If yes, execute step 4.6; otherwise, execute step 4.1.

[0021] Step 4.5: Use the parrot's fear of strangers to update the individual position. The specific calculation is as follows:

[0022]

[0023] Step 4.6: Determine whether the current number of iterations t is greater than t mf , otherwise, execute step 4.7. If yes, select several individual variables with poor fitness values, and use the improved adaptive crossover operator strategy to update the selected individual variables. The specific update strategy is calculated as follows:

[0024]

[0025] Where: X w1 , X w2 is the position of the offspring individuals generated by the adaptive crossover operator, X w is the randomly selected individual position with poor fitness value, β is the adaptive adjustment coefficient, and by comparing the fitness values ​​of the two offspring, the offspring individual with higher fitness value enters the next iteration;

[0026] Step 4.7: In addition to the optimal individual variable, several individual variables are randomly selected and the improved Archimedean spiral optimization strategy is used for individual update. The specific update strategy is calculated as follows:

[0027] X A (t+1)=|X best -X A (t)|e bl cos(2πl)+|X worst -X A (t)|e bl sin(2πl),

[0028] Where X A (t) is a randomly selected individual variable, X worst is the position of the individual with the worst fitness value, b is the spiral shape constant, l is a random number between 0 and 1, and several l values ​​are randomly generated to guide the selected individuals to explore in l different spatial domains according to the update formula in step 4.7, and select the individual with the highest fitness value to participate in the next round of iteration;

[0029] Step 4.8: Determine whether the current number of iterations t is greater than t f , otherwise, execute step 5. If yes, select the optimal individual variable, and use the dual-scale variable step size "quasi-gradient method" to update the optimal individual variable. The basic idea of ​​this update strategy is: take the optimal individual position as the origin, the distance to the origin as the scale, set two different scales with large size differences, and form a large-scale exploration space and a small-scale development space. First, perform a multi-directional multi-step search in the small-scale development space (multi-directional multi-step search refers to searching along multiple set directions, obtaining the optimal search direction and performing multi-step searches along this direction). When there is a better individual in the small-scale space, select the best evolutionary individual from all the better individuals and conduct a multi-step search along this direction until the termination condition is met (the maximum number of exploration steps is reached or the best individual in the current direction is not further optimized). Then, repeat the above multi-directional multi-step search method with the new optimal individual position as the origin until the number of trials reaches the set upper limit. If a better individual is found in the small-scale development space, the evolution process ends. Otherwise, perform a multi-directional multi-step search in the large-scale exploration space.

[0030] Step 5: Determine whether the maximum number of iterations has been reached. If yes, execute step 6; otherwise, execute step 4.

[0031] Step 6: Output the thrust distribution result, update the ship state parameters, and jump to step 2 to perform the thrust distribution calculation for the next cycle.

[0032] Preferably, in step 2, the calculation method for constructing the propeller equation constraint term is as follows:

[0033] The constraint term of the propulsion equation based on the extended thrust is expressed as: e T e =τ, where B e =[B e1 ,B e2 ,…,B en ] is the thruster position construction matrix, T e = = [T e1 ,T e2 ,…,T en ] is the extended thrust vector, n is the total number of thrusters, τ is the expected force and torque, and B of the i-th thruster ei and T ei The specific calculation method is:

[0034] The thruster angle is not fixed: T ei =[T xi ,T yi ] T ,

[0035] The thruster angle is fixed or the rotation is stuck, and -45°≤α i ≤45°, or 135°≤α i ≤180°, or -180°<α i ≤-135°:

[0036] B ei =[1 tan(α i ) -l yi +l xi tan(α i ), T ei =T xi

[0037] The propeller angle is fixed or the rotation is stuck. For other angles: B ei =[cot(α i ) 1 -l yi cot(α i )+l xi ] T , T ei =T yi ,

[0038] Among them, l xi and l yi is the thruster hull coordinate, T xi and T yi is the thrust along the x-axis and y-axis, α i is the azimuth.

[0039] Preferably, in step 3 and step 4, the individual fitness value is calculated as follows:

[0040]

[0041] J is the fitness value, T i is the thrust of propeller i; Δα is the azimuth change rate vector; B(α) is the position construction matrix (where the expression of the i-th column is [cos(α i ) sin(α i ) -l yi cos(α i )+l xi sin(α i )] T ), S is the slack variable, Q and Ω are the weight coefficient matrices, K pti is the coefficient of relationship between the power and thrust of the ith propeller, δ is the singular coefficient, ε is the minimum value to prevent the denominator from being zero, and w i is the penalty coefficient of thruster i, A thri is the thrust fault flag of thruster i. When there is a thrust fault, the flag is 1, otherwise it is 0. The method for setting the parameters for calculating the individual fitness value of the thruster is as follows: when there is a thrust fault in the thruster, the thrust limit value of the thruster is set to the allowable extreme value in the fault mode, the thrust fault flag is set to 1, and other parameter values ​​are set according to the initial parameters in step 1; in other cases, all parameter values ​​are set according to the initial parameters in step 1, and all thruster thrust fault flags are set to 0.

[0042] Preferably, in step 4, the calculation method of the dual-scale variable step size "quasi-gradient method" is as follows:

[0043] Step 4.8.1. Initialize parameters, including two scale constants D small and D large , the number of search directions N d , the maximum number of exploration steps N S , the maximum number of attempts N t , individual dimension D, optimal individual position X best wait;

[0044] Step 4.8.2: Use the improved Tent chaos mapping method to generate a vector X with strong dispersion r , the calculation method of the improved Tent chaos mapping method is:

[0045] Where 1≤j≤D, 1≤i≤N d , c rj is a set value or a random number between 0 and 1. j When (i+1) is greater than 1, the value is 1;

[0046] Step 4.8.3: Calculate X in the small-scale development space best Along the N d Update the value in each direction, where the i-th direction, 1≤i≤N d , the update method for each direction is as follows: X ri The value is X r The i-th column vector of N d Update the fitness values ​​of individual variables in each direction, select the optimal fitness value and determine whether the optimal individual variable is optimized. If yes, execute step 4.8.5, otherwise execute step 4.8.4;

[0047] Step 4.8.4: Calculate X in large-scale exploration space best Along the N d Update the value in each direction, where the i-th direction, 1≤i≤N d , the update method for each direction is as follows: Calculate N d Update the fitness values ​​of individual variables in each direction, select the optimal fitness value and determine whether the optimal individual variable is optimized. If yes, execute step 4.8.5, otherwise execute step 4.8.7;

[0048] Step 4.8.5: Continue optimizing the optimal individual along the optimal direction. If the optimal individual is optimized and the maximum number of exploration steps N is not reached, S Then repeat step 4.8.5 to optimize the best individual, otherwise go to step 4.8.6;

[0049] Step 4.8.6: Has the maximum number of attempts N been reached? t , otherwise return to step 4.8.2, if yes, execute step 4.8.7;

[0050] Step 4.8.7: The optimization of the best individual is completed and the optimization result is output.

[0051] Compared with the prior art, the present invention has the following beneficial effects:

[0052] (1) Compared with traditional thrust allocation methods (such as sequential quadratic programming and pseudo-inverse method), the thrust allocation method proposed in the present invention can better improve the thrust allocation performance and achieve better thrust allocation effect under complex working conditions.

[0053] (2) Compared with many existing swarm intelligence optimization thrust allocation algorithms, the thrust allocation method proposed in the present invention obtains a more stable and reliable thrust allocation solution with a smaller number of iterations.

[0054] (3) The thrust distribution algorithm of the present invention takes into account the propeller failure situation, so that the propeller can obtain a better thrust distribution result even under certain failure conditions, thereby ensuring the safety and stability of the ship's operation and reducing the impact of the failure. BRIEF DESCRIPTION OF THE DRAWINGS

[0055] The accompanying drawings are used to provide further understanding of the present invention and constitute a part of the specification. They are used to explain the present invention together with the embodiments of the present invention and do not constitute a limitation of the present invention.

[0056] Figure 1 It is a schematic diagram of the process of the present invention.

[0057] Figure 2 4 is a graph showing the longitudinal force output in Example 1 of the present invention.

[0058] Figure 3 It is a lateral force output curve diagram in Example 1 of the present invention.

[0059] Figure 4 It is a curve diagram of the resultant torque output in Example 1 of the present invention.

[0060] Figure 5 4 is a graph showing an absolute value of an error in Embodiment 1 of the present invention.

[0061] Figure 6 It is a curve diagram of propeller power consumption in Example 1 of the present invention.

[0062] Figure 7 It is a longitudinal force output curve diagram in Example 2 of the present invention.

[0063] Figure 8 It is a lateral force output curve diagram in Example 2 of the present invention.

[0064] Fig. 9 It is a curve diagram of the resultant torque output in Example 2 of the present invention.

[0065] Fig.10 4 is a graph showing an absolute value of an error in Embodiment 2 of the present invention.

[0066] Fig.11 This is a curve diagram of propeller power consumption in Example 2 of the present invention. DETAILED DESCRIPTION

[0067] In order to make the purpose, technical solution and advantages of the present invention more clearly understood, the present invention is further described in detail below in conjunction with the accompanying drawings and embodiments. Of course, the specific embodiments described here are only used to explain the present invention and are not used to limit the present invention.

[0068] Example 1

[0069] In this embodiment, all propellers work normally in the first 99 sampling cycles, and in the last 21 sampling cycles, propellers 1 and 3 work normally, while propeller 2 has a rotation failure and is stuck.

[0070] See also Figures 2 to 6 , Table 2, This embodiment provides a technical solution, a method for distributing thrust in a ship propeller failure mode based on an improved Parrot algorithm, comprising the following steps:

[0071] Step 1. Initialize various ship parameters: propeller position parameters, initial state values, physical limit values ​​(including propeller thrust and angle limit values, thrust change rate and angle change rate limit values), propeller power and thrust relationship coefficient. Some parameters are shown in Table 1 below:

[0072] Table 1 Some parameters of the thruster

[0073]

[0074] Step 2: Determine whether the current thruster has a fault and the fault type (the faults mainly include thrust faults and rotation faults), construct thruster equation constraint items (the thrust of each thruster adopts the extended thrust mode, and when the thruster angle is fixed or the rotation fault is stuck, the x-axis thrust or y-axis thrust of the thruster is adopted), and calculate the specific solution and general solution of the thruster equation constraint items (the general solution coefficient is selected as an individual variable);

[0075] In this embodiment, when the propeller is not stuck due to rotation failure, the equality constraint term is constructed as follows:

[0076] τ is the desired force and torque, T xi and T yi is the thrust of the ith thruster along the x-axis and y-axis,

[0077] In this embodiment, when the propeller 2 is stuck due to a rotation failure in the 100th to 120th sampling period, if -45°≤α2≤45°, or 135°≤α2≤180°, or -180°≤α2≤-135°, the equality constraint is constructed as follows:

[0078]

[0079] In this embodiment, when the propeller 2 is stuck due to a rotation failure in the 100th to 120th sampling period, if 45°<α2<135°, or -135°<α2<-45°, the equality constraint term is constructed as follows:

[0080]

[0081] Step 3: According to the current working condition of the thruster, initialize the parameters of the improved Parrot optimization algorithm: population size N = 40, maximum number of iterations t max , population upper and lower boundaries, time threshold t mf and t f (This embodiment t mf =0.5tmax , t f =0.7t max ), individual dimension, initialize the population, calculate the population fitness value and sort it in ascending order, obtain the optimal individual fitness value and its position, etc.;

[0082] The calculation method of individual fitness value is as follows:

[0083]

[0084] J is the fitness value, T i is the thrust of propeller i; Δα is the azimuth change rate vector; B(α) is the position construction matrix (where the expression of the i-th column is [cos(α i ) sin(α i ) -l yi cos(α i )+l xi sin(α i )] T ); S is the slack variable; Q and Ω are weight coefficient matrices (in this example, they are 10000 and 2 diagonal matrices respectively); Kpt i is the power-thrust relationship coefficient of the i-th propeller, δ is the singular coefficient (the value is 0.1 in this embodiment), and ε is the minimum value to prevent the denominator from being zero (the value is 0.0000001 in this embodiment);

[0085] Step 4: Use the improved Parrot optimization algorithm to update individual variables. The specific implementation steps are as follows:

[0086] Step 4.1, select the individual variables to be updated in turn, and select the individual update method by randomly generating four integers 1, 2, 3 and 4. When the random number is 1, execute step 4.2, when the random number is 2, execute step 4.3, when the random number is 3, execute step 4.4, and when the random number is 4, execute step 4.5;

[0087] Step 4.2: Use the parrot's foraging behavior to update the individual position. The specific calculation is as follows:

[0088]

[0089] Among them, X(t) is the current position of the individual, X best represents the current best individual position, Levy (D) is used to represent the parrot's Levy flight behavior, D refers to the individual dimension, P mean indicates the position of the entire group, t is the number of iterations, rand(0,1) is a random number between 0 and 1, and X mean(t) is the average value of the current position of the entire group. After the individual is updated, it is determined whether the 40 individual variables are updated. If yes, execute step 4.6, otherwise execute step 4.1;

[0090] Step 4.3: Use the parrot's stay behavior to update the individual position. The specific calculation is as follows:

[0091] X(t+1)=X(t)+X best Levy(D)+rand(0,1)ones(1,D),

[0092] Among them, rand(0,1)ones(1,D) indicates the position of random stay on the owner, and ones(1,D) means that the D-dimensional data is generated as a 1 vector. After the individual is updated, it is determined whether the N individual variables are updated. If yes, execute step 4.6, otherwise execute step 4.1;

[0093] Step 4.4: Use the parrot's communication behavior to update the individual position. The specific calculation is as follows:

[0094] When the random generation probability p≤0.5, the update method is: X(t+1)=0.2rand(0,1)(1-t / t max )(X(t)-X mean (t)),

[0095] On the contrary, when p>0.5, the update method is: X(t+1)=0.2rand(0,1)e^(-t / (rand(0,1)t max )),

[0096] After the individual is updated, determine whether the 40 individual variables have been updated. If yes, execute step 4.6; otherwise, execute step 4.1.

[0097] Step 4.5: Use the parrot's fear of strangers to update the individual position. The specific calculation is as follows:

[0098]

[0099] Step 4.6: Determine whether the current number of iterations t is greater than t mf , otherwise, execute step 4.7. If yes, select several individual variables with poor fitness values, and use the improved adaptive crossover operator strategy to update the selected individual variables. The specific update strategy is calculated as follows:

[0100]

[0101] Where: X w1 , X w2 is the position of the offspring individuals generated by the adaptive crossover operator, X wis the randomly selected individual position with poor fitness value, β is the adaptive adjustment coefficient, and by comparing the fitness values ​​of the two offspring, the offspring individual with higher fitness value enters the next iteration;

[0102] Step 4.7: In addition to the optimal individual variable, several individual variables are randomly selected and the improved Archimedean spiral optimization strategy is used for individual update. The specific update strategy is calculated as follows:

[0103] X A (t+1)=|X best -X A (t)|e bl cos(2πl)+|X worst -X A (t)|e bl sin(2πl),

[0104] Where X A (t) is a randomly selected individual variable, X worst is the position of the individual with the worst fitness value. In this embodiment, b takes the value of 1, and l is a random number between 0 and 1. By randomly generating several l values, the selected individuals are guided to explore in l different spatial domains according to the update formula in step 4.7, and the individuals with the highest fitness value are selected to participate in the next round of iteration;

[0105] Step 4.8: Determine whether the current number of iterations t is greater than t f , otherwise, execute step 5. If yes, select the optimal individual variable and use the dual-scale variable step size "quasi-gradient method" to update the optimal individual variable. The specific implementation steps are as follows:

[0106] Step 4.8.1. Initialize parameters, including two scale constants D small =0.4 and D large =4, number of search directions N d =16, maximum number of exploration steps N S =8, maximum number of attempts N t =6, individual dimension D, optimal individual position X best wait;

[0107] Step 4.8.2: Use the improved Tent chaos mapping method to generate a vector X with strong dispersion r , the calculation method of the improved Tent chaos mapping method is:

[0108] Where 1≤j≤D, 1≤i≤N d , c rj is a set value or a random number between 0 and 1. j When (i+1) is greater than 1, the value is 1;

[0109] Step 4.8.3: Calculate X in the small-scale development space best Along the N d Update the value in each direction, where the i-th (1≤i≤N d ), the update method for each direction is as follows: X ri The value is X r The i-th column vector of N d Update the fitness values ​​of individual variables in each direction, select the optimal fitness value and determine whether the optimal individual variable is optimized. If yes, execute step 4.8.5, otherwise execute step 4.8.4;

[0110] Step 4.8.4: Calculate X in large-scale exploration space best Along the N d Update the value in each direction, where the i-th (1≤i≤N d ), the update method for each direction is as follows: Calculate N d Update the fitness values ​​of individual variables in each direction, select the optimal fitness value and determine whether the optimal individual variable is optimized. If yes, execute step 4.8.5, otherwise execute step 4.8.7;

[0111] Step 4.8.5: Continue optimizing the optimal individual along the optimal direction. If the optimal individual is optimized and the maximum number of exploration steps N is not reached, S Then repeat step 4.8.5 to optimize the best individual, otherwise go to step 4.8.6;

[0112] Step 4.8.6: Has the maximum number of attempts N been reached? t , otherwise return to step 4.8.2, if yes, execute step 4.8.7;

[0113] Step 4.8.7, the best individual optimization is completed, and the optimization result is output;

[0114] Step 5: Determine whether the maximum number of iterations has been reached. If yes, execute step 6; otherwise, execute step 4.

[0115] Step 6: Output the thrust distribution result, update the ship state parameters, and jump to step 3 to perform the thrust distribution calculation for the next cycle.

[0116] Table 2 Comparison of thrust allocation results of different swarm intelligent optimization algorithms (1)

[0117]

[0118] In Table 2, MSPO is the implementation method of this embodiment, the maximum number of iterations is 10, the maximum number of iterations of the SSA algorithm is 100, and the maximum number of iterations of other algorithms is 30. The cost functions of all algorithms are set the same, and each optimization algorithm runs continuously for 30 times. It can be seen from the data in Table 2 that only the implementation method of this embodiment has a small difference in the absolute value and sum of output deviations and the maximum and minimum values ​​of the total power at different running times, while the other algorithms have a large difference, indicating that the output thrust distribution result of the implementation method of this embodiment has the best stability. In addition, the total absolute value of the deviation is very small, indicating that the accuracy and reliability of the thrust distribution result of the implementation method of this embodiment are the best among all the compared algorithms. From the total power data, it can be seen that the average energy consumption of the implementation method of this embodiment is also the lowest. The above results show that the implementation method of this embodiment has the ability to quickly find the global optimization under the condition of a low number of iterations, the stability and reliability of the optimization result are high, and the thrust distribution accuracy is high and the energy consumption is low. Figures 2 to 6 The output result of the algorithm in this embodiment is compared with the traditional SQP algorithm and PI-PSO algorithm (the algorithm is a combination of pseudo-inverse method and particle swarm algorithm. If the pseudo-inverse method is used alone in this example, the expected value of thrust distribution cannot be tracked). Figures 2 to 5 It can be seen that the accuracy of the algorithm in this embodiment is higher than that of the other two optimization algorithms (the total absolute value of the deviation of the SQP algorithm is 73.788, the total absolute value of the deviation of the PI-PSO algorithm is 77.754, and the total absolute value of the deviation of the MSPO algorithm is 1.978), which shows that the traditional sequential quadratic programming method and pseudo-inverse method are prone to fall into local extreme points, and the thrust distribution accuracy is relatively poor. The implementation method of this embodiment has a strong global optimization capability and can obtain a higher-precision thrust distribution solution. Figure 6 From the total power consumption of data calculation (the total power consumed by the SQP algorithm is 30.155 kW, the total power consumed by the MSPO algorithm is 30.073 kW, and the total power consumed by the PI-PSO algorithm is 30.105 kW), it can be seen that the implementation method of this embodiment is more conducive to reducing energy consumption.

[0119] Example 2

[0120] In this embodiment, propeller 2 works normally, propeller 1 has a rotation fault and can only rotate within the range of -1 to 1 radian, and propeller 3 has a thrust fault and the maximum thrust can only reach 30% of normal operation.

[0121] See also Figures 7 to 11 , Table 3, the technical solution provided by this embodiment is, a thrust distribution method under a ship propeller fault mode based on an improved parrot algorithm, except for step 2 and step 3, the steps of this embodiment are arranged the same as those in embodiment 1, and step 2 and step 3 of this embodiment are arranged as follows:

[0122] Step 2: Determine whether the current thruster has a fault and the fault type (the faults mainly include thrust faults and rotation faults), construct thruster equation constraint items (the thrust of each thruster adopts the extended thrust mode, and when the thruster angle is fixed or the rotation fault is stuck, the x-axis thrust or y-axis thrust of the thruster is adopted), and calculate the specific solution and general solution of the thruster equation constraint items (the general solution coefficient is selected as an individual variable);

[0123] In this embodiment, the equality constraint term is constructed as follows:

[0124] τ is the desired force and torque, T xi and T yi is the thrust of the ith thruster along the x-axis and y-axis,

[0125] Step 3: According to the current working condition of the thruster, initialize the parameters of the improved Parrot optimization algorithm: population size N = 40, maximum number of iterations t max , population upper and lower boundaries, time threshold t mf and t f (This embodiment t mf =0.5t max , t f =0.7t max ), individual dimension, initialize the population, calculate the population fitness value and sort it in ascending order, obtain the optimal individual fitness value and its position, etc.;

[0126] The calculation method of individual fitness value is as follows:

[0127]

[0128] J is the fitness value, T i is the thrust of propeller i; Δα is the azimuth change rate vector; B(α) is the position construction matrix (where the expression of the i-th column is [cos(α i ) sin(α i ) -l yi cos(α i )+l xi sin(α i )] T ); S is the slack variable; Q and Ω are weight coefficient matrices (in this example, they are 10000 and 2 diagonal matrices respectively); K pti is the power-thrust relationship coefficient of the ith propeller, δ is the singular coefficient (the value is 0.1 in this embodiment), ε is the minimum value to prevent the denominator from being zero (the value is 0.0000001 in this embodiment), w i is the penalty coefficient of thruster i (set to 20 in this embodiment), A thri is the thrust fault indicator of thruster i (in this embodiment, Athr3 =1, others are 0), the thrust limit value of thruster thrust 3 is reset to 30% of the normal working condition;

[0129] Table 3 Comparison of thrust allocation results of different swarm intelligent optimization algorithms (2)

[0130]

[0131]

[0132] In this embodiment, the maximum number of iterations, the cost function setting and the number of runs are the same as those in Embodiment 1. From the data in Table 3, it can be seen that the implementation method of this embodiment is still the best in terms of stability, reliability, accuracy and energy consumption compared with other algorithms. It further shows that under different types of fault modes, the implementation method of this embodiment has the ability to quickly find the optimal solution with a low number of iterations, improves the thrust distribution accuracy, reduces energy consumption, and ensures the reliability and safety of ship operation. Figures 7 to 10 It can be seen that the other two optimization algorithms have large deviations. The total absolute value of the deviations of the SQP algorithm is 85.087, and the total absolute value of the deviations of the PI-PSO algorithm is 67.078), while the total absolute value of the deviations of the MSPO algorithm in this embodiment is only 0.886, which further illustrates that the traditional sequential quadratic programming method and pseudo-inverse method are prone to fall into local extreme points, and the thrust distribution accuracy is relatively poor under some working conditions. The implementation method of this embodiment has a stronger global optimization capability and can obtain a high-precision thrust distribution solution. Fig.11 The total power consumption of data calculation (the total power consumed by the SQP algorithm is 15.967kW, the total power consumed by the MSPO algorithm is 15.505kW, and the total power consumed by the PI-PSO algorithm is 16.036kW. It can be seen that the implementation method of this embodiment has lower energy consumption.

[0133] From the analysis of the above two embodiments, it can be seen that the MSPO thrust distribution method proposed in this embodiment can improve the thrust distribution accuracy and reduce energy consumption. The thrust distribution result is stable and reliable. Even under certain thruster failure conditions, a better thrust distribution solution can be obtained to ensure the safety and stability of the ship's operation. In addition, the thrust distribution algorithm of this embodiment can modify each index coefficient according to actual needs to adapt to the thrust distribution optimization of different working conditions.

[0134] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, improvements, etc. made within the spirit and principle of the present invention should be included in the protection scope of the present invention.

Claims

1. A thrust distribution method under ship propeller failure mode based on improved Parrot algorithm, characterized in that: The steps include: Step 1: Initialize various ship parameters: propeller position parameters, initial state values, physical limit values, propeller power and thrust relationship coefficients; The physical limit values ​​include thrust and angle limit values, thrust change rate and angle change rate limit values ​​of the propeller; Step 2: Determine whether the current thruster has a fault and the fault type. The fault types include thrust fault and rotation fault. Construct the thruster equation constraint item. The thrust of each thruster adopts the extended thrust mode. When the thruster angle is fixed or the rotation fault is stuck, the x-axis thrust or y-axis thrust of the thruster is used. Calculate the specific solution and general solution of the thruster equation constraint item. The general solution coefficient is selected as the individual variable. Step 3: According to the current working condition of the thruster, initialize the parameters of the improved Parrot optimization algorithm: population size N, maximum number of iterations t max , individual dimension, population upper and lower boundaries, time threshold t mf and t f , t mf <t f <t max , initialize the population, calculate the fitness value of the population and sort it in ascending order, and obtain the optimal individual fitness value and its position; Step 4: Use the improved Parrot optimization algorithm to update individual variables. The steps are as follows: Step 4.1, select the individual variables to be updated in turn, and select the individual update method by randomly generating four integers 1, 2, 3 and 4. When the random number is 1, execute step 4.2, when the random number is 2, execute step 4.3, when the random number is 3, execute step 4.4, and when the random number is 4, execute step 4.5; Step 4.2: Use the parrot's foraging behavior to update the individual position. The calculation is as follows: Among them, X(t) is the current position of the individual, X best represents the current best individual position, Levy (D) is used to represent the parrot's Levy flight behavior, D refers to the individual dimension, P mean indicates the position of the entire group, t is the number of iterations, rand(0,1) is a random number between 0 and 1, and X mean (t) is the average value of the current position of the entire group. After the individual is updated, it is determined whether the N individual variables are updated. If yes, execute step 4.6, otherwise execute step 4.1; Step 4.3: Use the parrot's stay behavior to update the individual position. The calculation is as follows: X(t+1)=X(t)+X best Levy(D)+rand(0,1)ones(1,D), Among them, rand(0,1)ones(1,D) indicates the position of random stay on the owner, and ones(1,D) means that the D-dimensional data is generated as a 1 vector. After the individual is updated, it is determined whether the N individual variables are updated. If yes, execute step 4.6, otherwise execute step 4.1; Step 4.4: Use the parrot's communication behavior to update the individual position. The calculation is as follows: When the random generation probability p≤0.5, the update method is: X(t+1)=0.2rand(0,1)(1-t / t max )(X(t)-X mean (t)), On the contrary, when p>0.5, the update method is: X(t+1)=0.2rand(0,1)e^(-t / (rand(0,1)t max )), After the individual is updated, determine whether the N individual variables are updated. If yes, execute step 4.6; otherwise, execute step 4.

1. Step 4.5: Use the parrot's fear of strangers to update the individual position. The calculation is as follows: Step 4.6: Determine whether the current number of iterations t is greater than t mf , otherwise, execute step 4.

7. If yes, select several individual variables with poor fitness values, and use the improved adaptive crossover operator strategy to update the selected individual variables. The update strategy is calculated as follows: Where: X w1 , X w2 is the position of the offspring individuals generated by the adaptive crossover operator, X w is the randomly selected individual position with poor fitness value, β is the adaptive adjustment coefficient, and by comparing the fitness values ​​of the two offspring, the offspring individual with higher fitness value enters the next iteration; Step 4.7: In addition to the optimal individual variable, several individual variables are randomly selected and the improved Archimedean spiral optimization strategy is used to update the individual variables. The update strategy is calculated as follows: X A (t+)=|X best -X A (t)|e bl cos(2πl)+|X worst -X A (t)|e bt sin(2πl), Where X A (t) is a randomly selected individual variable, X worst is the position of the individual with the worst fitness value, b is the spiral shape constant, l is a random number between 0 and 1, and several l values ​​are randomly generated to guide the selected individuals to explore in l different spatial domains according to the update formula in step 4.7, and select the individual with the highest fitness value to participate in the next round of iteration; Step 4.8: Determine whether the current number of iterations t is greater than t f , otherwise, execute step 5. If yes, select the optimal individual variable, and use the dual-scale variable step size "quasi-gradient method" to update the optimal individual variable. The basic idea of ​​this update strategy is: take the optimal individual position as the origin, the distance to the origin as the scale, set two different scales with large size differences, constitute a large-scale exploration space and a small-scale development space, first perform a multi-directional multi-step search in the small-scale development space (multi-directional multi-step search refers to searching along multiple set directions, obtaining the optimal search direction and performing multi-step searches along this direction), when there is a better individual in the small-scale space, select the best evolutionary individual from all the better individuals and conduct a multi-step search along this direction until the termination condition is met, the maximum number of exploration steps is reached or the best individual in the current direction is not further optimized, then repeat the above multi-directional multi-step search method with the new optimal individual position as the origin until the number of trials reaches the set upper limit, if a better individual is found in the small-scale development space, the evolution process ends, otherwise a multi-directional multi-step search is performed in the large-scale exploration space; Step 5: Determine whether the maximum number of iterations has been reached. If yes, execute step 6; otherwise, execute step 4. Step 6: Output the thrust distribution result, update the ship state parameters, and jump to step 2 to perform the thrust distribution calculation for the next cycle.

2. The thrust distribution method under the ship propeller failure mode based on the improved Parrot algorithm according to claim 1 is characterized in that: In step 2, the calculation method for constructing the thruster equation constraint term is as follows: The constraint term of the propulsion equation based on the extended thrust is expressed as: e T e =τ, where B e =[B e1 ,B e2 ,…,B en ] is the thruster position construction matrix, T e = = [T e1 ,T e2 ,…,T en ] is the extended thrust vector, n is the total number of thrusters, τ is the expected force and torque, and B of the i-th thruster ei and T ei The calculation method is: The thruster angle is not fixed: T ei =[T xi ,T yi ] T , The thruster angle is fixed or the rotation is stuck, and -45°≤α i ≤45°, or 135°≤α i ≤180°, or -180°<α i ≤-135°: B ei \[1 to(α i ) -l yi +l xi tan(α i )] T ,T ei =T xi , The propeller angle is fixed or the rotation is stuck. For other angles: B ei =[cot(α i ) 1 -l yi cot(α i )+l xi ] T , T ei =T yi , Among them, l xi and l yi is the thruster hull coordinate, T xi and T yi is the thrust along the x-axis and y-axis, α i is the azimuth.

3. The thrust distribution method under the ship propeller failure mode based on the improved Parrot algorithm according to claim 1 is characterized in that: In step 3 and step 4, the calculation method of individual fitness value is as follows: J is the fitness value, T i is the thrust of propeller i; Δα is the azimuth change rate vector; B(α) is the position construction matrix (where the expression of the i-th column is [cos(α i ) sin(α i ) -l yi cos(α i )+l xi sin(α i )] T ), S is the slack variable, Q and Ω are the weight coefficient matrices, K pti is the coefficient of relationship between the power and thrust of the ith propeller, δ is the singular coefficient, ε is the minimum value to prevent the denominator from being zero, and w i is the penalty coefficient of thruster i, A thri is the thrust fault flag of thruster i. When there is a thrust fault, the flag is 1, otherwise it is 0. The method for setting the parameters for calculating the individual fitness value of the thruster is as follows: when there is a thrust fault in the thruster, the thrust limit value of the thruster is set to the allowable extreme value in the fault mode, the thrust fault flag is set to 1, and other parameter values ​​are set according to the initial parameters in step 1; in other cases, all parameter values ​​are set according to the initial parameters in step 1, and all thruster thrust fault flags are set to 0.

4. The method for distributing thrust under a ship propeller failure mode based on the improved Parrot algorithm according to claim 1, characterized in that: In step 4, the calculation method of the dual-scale variable step size "quasi-gradient method" is as follows: Step 4.8.

1. Initialize parameters, including two scale constants D small and D large , the number of search directions N d , the maximum number of exploration steps N S , the maximum number of attempts N t , individual dimension D, optimal individual position X best ; Step 4.8.2: Use the improved Tent chaos mapping method to generate a vector X with strong dispersion r , the calculation method of the improved Tent chaos mapping method is: Where 1≤j≤D, 1≤i≤N d , c rj is a set value or a random number between 0 and 1. j When (i+1) is greater than 1, the value is 1; Step 4.8.3: Calculate X in the small-scale development space best Along the N d Update the value in each direction, where the i-th direction, 1≤i≤N d , the update method for each direction is as follows: X ri The value is X r The i-th column vector of N d Update the fitness values ​​of individual variables in each direction, select the optimal fitness value and determine whether the optimal individual variable is optimized. If yes, execute step 4.8.5, otherwise execute step 4.8.4; Step 4.8.4: Calculate X in large-scale exploration space best Along the N d Update the value in each direction, where the i-th direction, 1≤i≤N d , the update method for each direction is as follows: Calculate N d Update the fitness values ​​of individual variables in each direction, select the optimal fitness value and determine whether the optimal individual variable is optimized. If yes, execute step 4.8.5, otherwise execute step 4.8.7; Step 4.8.5: Continue optimizing the optimal individual along the optimal direction. If the optimal individual is optimized and the maximum number of exploration steps N is not reached, S Then repeat step 4.8.5 to optimize the best individual, otherwise go to step 4.8.6; Step 4.8.6: Has the maximum number of attempts N been reached? t , otherwise return to step 4.8.2, if yes, execute step 4.8.7; Step 4.8.7: The optimization of the best individual is completed and the optimization result is output.

Images