Unmanned vehicle yaw motion optimal control method based on symplectic pseudo-spectral method

Abstract

According to the unmanned vehicle yaw motion optimal control method based on the symplectic pseudo-spectral method, the optimal control problem of an original continuous time system is approximated by executing the pseudo-spectral method, and the designed line search and convex optimization technology can effectively guarantee the iteration speed and the numerical stability; the problem caused by improper initial grid division can be improved by combining model precision test and a grid adaptive updating strategy, and a smooth and continuous control law can be generated finally, so that more accurate vehicle system yaw motion control is realized.

Description

Technical Field

[0001] The present invention belongs to the technical field of yaw motion control of vehicle trajectory tracking, and in particular relates to a yaw motion optimal control method using a symplectic pseudospectral method in an indirect method of optimal control problems as an optimization solver. Background Art

[0002] The vehicle yaw control problem is a complex dynamic optimization problem, which is particularly important for ensuring the vehicle driving stability and safety under complex conditions such as high-speed driving and emergency obstacle avoidance in autonomous driving. The commonly used optimization methods for vehicle yaw control problems in the prior art mainly include precise control algorithms based on dynamic models, intelligent control algorithms, and integrated optimization algorithms. Among them, precise control algorithms based on dynamic models, such as model predictive control and linear quadratic regulators, can achieve high control accuracy, but the huge amount of calculation makes it difficult to apply them in real-time systems; intelligent control algorithms such as fuzzy logic control and neural network control have good adaptability and robustness, but the design is complex and may have overfitting problems; integrated optimization algorithms such as methods that combine model predictive control with intelligent control can combine the advantages of the former two to achieve efficient and precise control, but there are still disadvantages of complex design and implementation processes. Summary of the invention

[0003] In view of this, in order to solve the above technical problems in the art, the present invention provides an optimal control method for yaw motion of an unmanned vehicle based on symplectic pseudospectral method, which specifically includes the following steps:

[0004] Step 1: Determine the optimal control problem of the yaw motion of the unmanned vehicle and its corresponding model, including the dynamic model related to the vehicle's center of mass sideslip angle, yaw rate, and the classic magic tire formula (MF98):

[0005]

[0006] where x = [β, ω, δ] T ,

[0007] ψ(x)=(ω-ω des ) 2

[0008]

[0009] Where β, ω, and δ are the vehicle's center of mass sideslip angle, yaw rate, and wheel angle, respectively; mass is the vehicle's mass, inertial z is the angular inertia of the vehicle around the z-axis center, F y,r 、F y,fare the lateral forces provided by the rear wheel and the front wheel respectively, t is the time, t s ,t f are the start and end time respectively, ω des is the target yaw rate, u is the control variable, x is the state variable, is its first-order derivative;

[0010] A virtual weight σ is introduced into the lateral force to describe the transition of the model from approximately linear to nonlinear:

[0011] F y =F y,linear σ-F y,nonlinear (σ-1)

[0012] F y,linear =F z DCBs y

[0013] F y,nonlinear =F z Dsin(C atan(B s y +E(s y B-atan(s y B))))

[0014] Among them, F z is the vertical load of the tire, B, C, D, E are the parameters of the classic magic tire formula;

[0015] Step 2: Hamiltonian function established for the optimal control problem According to the optimality condition, it is rewritten into a new Hamiltonian function form H = H (x, λ) containing only the state quantity x and the co-state quantity λ;

[0016] Step 3: In the time domain [t s ,t f ] Given an initial grid solution, [t s ,t f ] is divided into Q sub-intervals (the distances between each interval may be unequal), and in any j-th sub-interval, a generating function S is constructed based on the Hamiltonian function H. (j) and use N (j) The pseudo-spectrum of the order is used to approximate each generating function; based on the generating function S in the total time domain and the action and last action principle of the generating function in each sub-interval, a nonlinear equation group F(r)=0 of the first-order optimality condition is established, and then the Jacobian matrix K of the nonlinear equation group with respect to the combined state quantity r is calculated to obtain a basic iteration KΔr+F=0, and Δr is solved by iteration; wherein,

[0017]

[0018] x0 and x Q are the initial state and mode after pseudo-spectral approximation, Δr = -K -1 (r)F(r); the superscript “-” indicates that the pseudo-spectral approximation of the state quantity or co-state quantity is obtained;

[0019] Step 4: For the solution Δr obtained in the kth iteration k , regard it as the search direction of this iteration and define a new value function Φ, use the Armijo line search method to solve an optimal subproblem so that the increment of r decreases with respect to Φ;

[0020] Step 5: After the current iteration meets the predetermined conditions, check whether the virtual quantity weight σ needs to be updated. If it needs to be updated, return to step 2; if it does not need to be updated, proceed to the next step;

[0021] Step 6: Check the model accuracy at the pseudo-spectrum matching points and other points outside the given matching points. If the model accuracy does not meet the requirements, a new sub-interval division is obtained through the grid adaptive update strategy, and return to step 2; otherwise, the iteration ends and proceeds to the next step;

[0022] Step 7: Use the sampling controller to act on any moment in the time domain, restore the state quantity and corresponding control quantity at the current moment through the point information and pseudo-spectrum information, and use them as the output at the current moment.

[0023] Furthermore, in step 2, the following Hamiltonian function is first established:

[0024]

[0025] According to the first-order optimality condition Solve and get u=u(x,λ) and substitute it into A new Hamiltonian function form H = H (x, λ) is obtained which only contains the state quantity x and the co-state quantity λ.

[0026] Furthermore, in step 3, specifically for the time domain [t s ,t f The jth subinterval among the Q subintervals divided by ] is defined as:

[0027]

[0028] have Continuous function x in subinterval (j) (t) and λ (j) (t)N (j) The LGL-type pseudospectral approximations of the orders are:

[0029]

[0030] The generating function in the subinterval is:

[0031]

[0032] Use N (j) The LGL-type pseudo-spectrum of order is approximated to obtain:

[0033]

[0034] in, is H(x,λ) at τ k The value at

[0035] For S (j) Find the variation and apply the first-order optimality condition:

[0036]

[0037] a and b are the variational lower and upper bounds, respectively;

[0038] get:

[0039]

[0040] For x0 and x after pseudospectral approximation Q have:

[0041]

[0042] And for the LGL type pseudo-spectrum, we have:

[0043]

[0044] Furthermore, the first-order optimality condition for establishing the nonlinear equation system F(r)=0 in step 3 is specifically:

[0045]

[0046] When establishing the nonlinear equation system F(r)=0, first remember:

[0047]

[0048] Where m = 0, 1, 2, ..., N (j) , according to the above first-order optimality condition, we can derive:

[0049]

[0050] Taking partial derivative of f, we get:

[0051]

[0052] remember:

[0053]

[0054] The nonlinear equation system F(r) is obtained:

[0055]

[0056] remember:

[0057]

[0058] Then the Jacobian matrix K of equation F with respect to r is:

[0059]

[0060] According to the characteristics of the optimal control problem that the initial state is given and the terminal is free, we get:

[0061]

[0062] Therefore, for a given point r k At , the obtained F(r) and K(r) are used to perform Newton's method iteration calculation:

[0063] Δr k =-K -1 (r k )F(r k )

[0064] r k+1 =r k +Δr k .

[0065] Furthermore, in step 4, the value function Φ of the following form is specifically defined:

[0066] Φ(r)=||F(r)||2

[0067] When using the Armijo line search method to solve the optimal subproblem, the k-th result Δr k As the search direction, find a suitable α that satisfies:

[0068]

[0069] Furthermore, in step 5, the following predetermined conditions are specifically set:

[0070] Δr k ≤erro max or iter ≥ iter max

[0071] Among them, erro max is the maximum allowable error, iter, iter maxare the current number of iterations and the maximum number of iterations allowed respectively. The virtual quantity weight σ is updated once each time the above conditions are met.

[0072] Furthermore, in step six, the model accuracy is checked based on the following formula:

[0073]

[0074] Among them, e (j) is the model accuracy of the j-th subinterval, is the model error at any point selected within the subinterval, which is defined as follows:

[0075]

[0076] in, Calculated through pseudo-spectral transformation information;

[0077] The model accuracy must meet the following requirements:

[0078] If it is not satisfied, the grid adaptive update strategy is executed to obtain a new sub-interval division:

[0079] First, define the model curvature as:

[0080]

[0081] For all points in the current subinterval And a given parameter R, perform the following update steps:

[0082] ①If The subscripts max and ave represent the maximum value and the average value, respectively, dividing the subinterval into A new uniform subinterval;

[0083] ②If Then the number of pseudo-spectral points in the subinterval is increased to

[0084] Furthermore, in step 7, the process of restoring the state quantity at the current moment by using the point matching information and the pseudo spectrum information includes:

[0085] First, find the subinterval t∈[t j ,t j-1 ], and perform pseudo-spectral transformation to obtain τ:

[0086]

[0087] For the pseudo-spectral transformation parameters L and D corresponding to the sub-interval, define the sampling vector L τ , Dτ , and the following calculation formula is obtained:

[0088] L τ (k)=polyval([L(k,1);L(k,2);…;L(k,N (j) )],τ),k=0,1,2,…,N (j)

[0089] D τ (k)=polyval([D(k,1);D(k,2);…;D(k,N (j) )],τ),k=0,1,2,…,N (j)

[0090] Among them, polyval(p,τ) means calculating the value of polynomial p at τ, p is a vector of length n+1, and the elements are the coefficients of the n-degree polynomial in descending order;

[0091] According to the properties of the pseudo-spectrum, we can obtain:

[0092]

[0093]

[0094] For the case where higher-order derivatives need to be calculated, taking the second order as an example, the first-order derivative at the collocation point can be calculated using the formula Then calculate the second-order derivative at a given time t:

[0095]

[0096] The optimal control method for the yaw motion of an unmanned vehicle based on the symplectic pseudospectral method provided by the present invention approximates the optimal control problem of the original continuous-time system by executing the pseudospectral method. The designed line search and convex optimization techniques can effectively ensure the iteration speed and numerical stability. The combination of model accuracy verification and grid adaptive update strategy can improve the problems caused by improper initial grid division. The present invention can ultimately generate a smooth and continuous control law, thereby achieving more accurate yaw motion control of the vehicle system. BRIEF DESCRIPTION OF THE DRAWINGS

[0097] Figure 1 A framework diagram of the method provided by the present invention;

[0098] Figure 2 Based on the single optimization result in the example of the present invention;

[0099] Figure 3 The results are compared with the examples of the present invention and other prior arts. DETAILED DESCRIPTION

[0100] The technical solution of the present invention will be described clearly and completely below in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of the present invention.

[0101] The optimal control method for yaw motion of an unmanned vehicle based on the symplectic pseudospectral method provided by the present invention specifically comprises the following steps:

[0102] Step 1: Determine the optimal control problem of the yaw motion of the unmanned vehicle and its corresponding model, including the dynamic model related to the vehicle's center of mass sideslip angle, yaw rate, and the classic magic tire formula (MF98):

[0103]

[0104] where x = [β, ω, δ] T ,

[0105] ψ(x)=(ω-ω des ) 2

[0106]

[0107] Where β, ω, and δ are the vehicle's center of mass sideslip angle, yaw rate, and wheel angle, respectively; mass is the vehicle's mass, inertial z is the angular inertia of the vehicle around the z-axis center, F y,r 、F y,f are the lateral forces provided by the rear wheel and the front wheel respectively, t is the time, t s ,t f are the start and end time respectively, ω des is the target yaw rate, u is the control variable, x is the state variable, is its first-order derivative;

[0108] A virtual weight σ is introduced into the lateral force to describe the transition of the model from approximately linear to nonlinear:

[0109] F y =F y,linear σ-F y,nonlinear (σ-1)

[0110] F y,linear =F z DCBs y

[0111] F y,nonlinear =F zDsin(C atan(B s y +E(s y B-atan(s y B))))

[0112] Among them, F z is the vertical load of the tire, B, C, D, E are the parameters of the classic magic tire formula;

[0113] The virtual quantity weight σ decreases gradually with the iteration process. In the example of the present invention, the specific values ​​of σ selected with the iteration are [1, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0];

[0114] Step 2: Hamiltonian function established for the optimal control problem According to the optimality condition, it is rewritten into a new Hamiltonian function form H = H (x, λ) containing only the state quantity x and the co-state quantity λ;

[0115] Step 3: In the time domain [t s ,t f ] Given an initial grid solution, [t s ,t f ] is divided into Q sub-intervals (the distances between each interval may be unequal), and in any j-th sub-interval, a generating function S is constructed based on the Hamiltonian function H. (j) and use N (j) The pseudo-spectrum of the order is used to approximate each generating function; based on the generating function S in the total time domain and the action and last action principle of the generating function in each sub-interval, a nonlinear equation group F(r)=0 of the first-order optimality condition is established, and then the Jacobian matrix K of the nonlinear equation group with respect to the combined state quantity r is calculated to obtain a basic iteration KΔr+F=0, and Δr is solved by iteration; wherein,

[0116]

[0117] x0 and x Q are the initial state and mode after pseudo-spectral approximation, Δr = -K -1 (r)F(r); the superscript “-” indicates that the pseudo-spectral approximation of the state quantity or co-state quantity is obtained;

[0118] Step 4: For the solution Δr obtained in the kth iteration k , regard it as the search direction of this iteration and define a new value function Φ, use the Armijo line search method to solve an optimal subproblem so that the increment of r decreases with respect to Φ;

[0119] Step 5: After the current iteration meets the predetermined conditions, check whether the virtual quantity weight σ needs to be updated. If it needs to be updated, return to step 2; if it does not need to be updated, proceed to the next step;

[0120] Step 6: Check the model accuracy at the pseudo-spectrum matching points and other points outside the given matching points. If the model accuracy does not meet the requirements, a new sub-interval division is obtained through the grid adaptive update strategy, and return to step 2; otherwise, the iteration ends and proceeds to the next step;

[0121] Step 7: Use the sampling controller to act on any moment in the time domain, restore the state quantity and corresponding control quantity at the current moment through the point information and pseudo-spectrum information, and use them as the output at the current moment.

[0122] In a preferred embodiment of the present invention, in step 2, the following Hamiltonian function is first established:

[0123]

[0124] According to the first-order optimality condition Solve and get u=u(x,λ) and substitute it into A new Hamiltonian function form H = H (x, λ) is obtained which only contains the state quantity x and the co-state quantity λ.

[0125] Furthermore, in step 3, specifically for the time domain [t s ,t f The jth subinterval among the Q subintervals divided by ] is defined as:

[0126]

[0127] have

[0128] Continuous function x in subinterval (j) (t) and λ (j) (t)N (j) The LGL-type pseudospectral approximations of the orders are:

[0129]

[0130] The generating function in the subinterval is:

[0131]

[0132] Use N (j) The LGL-type pseudo-spectrum of order is approximated to obtain:

[0133]

[0134] in, is H(x,λ) at τ k The value at

[0135] For S (j) Find the variation and apply the first-order optimality condition:

[0136]

[0137] a and b are the variational lower and upper bounds, respectively;

[0138] get:

[0139]

[0140] For x0 and x after pseudospectral approximation Q have:

[0141]

[0142] And for the LGL type pseudo-spectrum, we have:

[0143]

[0144] In a preferred embodiment of the present invention, the first-order optimality condition for establishing the nonlinear equation system F(r)=0 in step 3 is specifically:

[0145]

[0146] When establishing the nonlinear equation system F(r)=0, first remember:

[0147]

[0148] Where m = 0, 1, 2, ..., N (j) , according to the above first-order optimality condition, we can derive:

[0149]

[0150] Taking partial derivative of f, we get:

[0151]

[0152] remember:

[0153]

[0154] The nonlinear equation system F(r) is obtained:

[0155]

[0156] remember:

[0157]

[0158] Then the Jacobian matrix K of equation F with respect to r is:

[0159]

[0160] According to the characteristics of the optimal control problem that the initial state is given and the terminal is free, we get:

[0161]

[0162] Therefore, for a given point r k At , the obtained F(r) and K(r) are used to perform Newton's method iteration calculation:

[0163] Δr k =-K -1 (r k )F(r k )

[0164] r k+1 =r k +Δr k .

[0165] In a preferred embodiment of the present invention, the value function Φ of the following form is specifically defined in step 4:

[0166] Φ(r)=||F(r)||2

[0167] When using the Armijo line search method to solve the optimal subproblem, the k-th result Δr k As the search direction, find a suitable α that satisfies:

[0168]

[0169] In a preferred embodiment of the present invention, the following predetermined conditions are specifically set in step 5:

[0170] Δr k ≤erro max or iter ≥ iter max

[0171] Among them, erro max is the maximum allowable error, iter, iter max are the current number of iterations and the maximum number of iterations allowed respectively. The virtual quantity weight σ is updated once each time the above conditions are met.

[0172] In a preferred embodiment of the present invention, the model accuracy is checked in step 6 based on the following formula:

[0173]

[0174] Among them, e (j) is the model accuracy of the j-th subinterval, is the model error at any point selected within the subinterval, which is defined as follows:

[0175]

[0176] in, Calculated through pseudo-spectral transformation information;

[0177] The model accuracy must meet the following requirements:

[0178] If it is not satisfied, the grid adaptive update strategy is executed to obtain a new sub-interval division:

[0179] First, define the model curvature as:

[0180]

[0181] For all points in the current subinterval And a given parameter R, perform the following update steps:

[0182] ①If The subscripts max and ave represent the maximum value and the average value, respectively, dividing the subinterval into A new uniform subinterval;

[0183] ②If Then the number of pseudo-spectral points in the subinterval is increased to

[0184] In a preferred embodiment of the present invention, the process of restoring the state quantity at the current moment by using the point matching information and the pseudo spectrum information in step seven includes:

[0185] First, find the subinterval t∈[t j ,t j-1 ], and perform pseudo-spectral transformation to obtain τ:

[0186]

[0187] For the pseudo-spectral transformation parameters L and D corresponding to the sub-interval, define the sampling vector L τ , D τ , and the following calculation formula is obtained:

[0188] L τ (k)=polyval([L(k,1);L(k,2);…;L(k,N (j) )],τ),k=0,1,2,…,N (j)

[0189] D τ (k)=polyval([D(k,1);D(k,2);…;D(k,N (j) )],τ),k=0,1,2,…,N (j)

[0190] Among them, polyval(p,τ) means calculating the value of polynomial p at τ, p is a vector of length n+1, and the elements are the coefficients of the n-degree polynomial in descending order;

[0191] According to the properties of the pseudo-spectrum, we can obtain:

[0192]

[0193] For the case where higher-order derivatives need to be calculated, taking the second order as an example, the first-order derivative at the collocation point can be calculated using the formula Then calculate the second-order derivative at a given time t:

[0194]

[0195] Figure 1 The system optional framework of the present invention when applied to the real vehicle control of an unmanned vehicle is shown, which is divided into two layers of controllers. The upper layer controller executes in sequence:

[0196] ① Set the allowable error, maximum number of iterations, vehicle system model parameters, pseudo-spectral transformation parameters and grid parameters.

[0197] ② The upper controller starts the loop, obtains vehicle motion information when the next loop starts, and initializes the iteration parameters.

[0198] ③Give an initial guess for the optimal solution and start iteration.

[0199] ④ For the solution Δr obtained in the kth iteration k , regard it as the search direction of this iteration, and use the Armijo line search technique to solve an optimal subproblem so that the increment of r decreases with respect to Φ.

[0200] ⑤ For the current iteration, check the virtual quantity weight. If the virtual quantity weight needs to be updated, return to step 2; if it does not need to be updated, proceed to the next step.

[0201] ⑥ Check the model accuracy. If the model accuracy does not meet the requirements, apply the grid adaptive update strategy to obtain a new sub-interval and return to step 3. Otherwise, end the iteration and output the iteration result, and wait for the next cycle.

[0202] The lower controller calculates the control quantity at the current moment and outputs it in each cycle.

[0203] Figure 2 The results of a single optimization in the example of the present invention are shown. Figure 3 By comparing the optimization results of this example with those obtained by applying the common linear model predictive controller LMPC, it can be seen that the present invention has a significant accuracy advantage over the prior art.

[0204] It should be understood that the size of the serial numbers of the steps in the embodiment of the present invention does not mean the order of execution. The execution order of each process should be determined by its function and internal logic, and should not constitute any limitation on the implementation process of the embodiment of the present invention.

[0205] Although embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions and variations may be made to the embodiments without departing from the principles and spirit of the present invention, and that the scope of the present invention is defined by the appended claims and their equivalents.

Claims

1. An optimal control method for yaw motion of an unmanned vehicle based on the symplectic pseudospectral method, characterized by: The specific steps include: Step 1: Determine the optimal control problem and corresponding model of the yaw motion of the unmanned vehicle, including the dynamic model related to the vehicle's center of mass sideslip angle, yaw rate, and the classic magic tire formula: x(t s )=x s Where x=[β,ω,δ] T , ψ(x)=(ω-ω des ) 2 Among them, β, ω, δ are the vehicle's center of mass sideslip angle, yaw rate and wheel angle, mass is the vehicle's mass, and intertial z is the angular inertia of the vehicle around the z-axis center, F y,r 、F y,f are the lateral forces provided by the rear wheel and the front wheel respectively, t is the time, t s ,t f are the start and end time respectively, ω des is the target yaw rate, u is the control quantity, x is the state quantity, is its first-order derivative; A virtual weight σ is introduced into the lateral force to describe the transition of the model from approximately linear to nonlinear: F y =F y,linear σ-F y,nonlinear (p-1) F y,linear =F z DCBs y F y,nonlinear =F z Dsin(Catan(Bs y +E(s y B-atan(s y B)))) Among them, F z is the vertical load of the tire, B, C, D, E are the parameters of the classic magic tire formula, s y is the slip rate; Step 2: Hamiltonian function established for the optimal control problem According to the optimality condition, it is rewritten into a new Hamiltonian function form H = H (x, λ) containing only the state quantity x and the co-state quantity λ; Step 3: In the time domain [t s ,t f ] Given an initial grid solution, [t s ,t f ] is divided into Q subintervals, and in any j-th subinterval, a generating function S is constructed based on the Hamiltonian function H. (j) and use N (j) The pseudo-spectrum of the order is used to approximate each generating function; based on the generating function S in the total time domain and the action and last action principle of the generating function in each sub-interval, a nonlinear equation group F(r)=0 of the first-order optimality condition is established, and then the Jacobian matrix K of the nonlinear equation group with respect to the combined state quantity r is calculated to obtain a basic iteration KΔr+F=0, and Δr is solved by iteration; wherein, x0 and x Q are the initial state and mode after pseudo-spectral approximation, Δr = -K -1 (t)F(r); the superscript "-" indicates that the pseudo-spectral approximation of the state quantity or co-state quantity is obtained; Step 4: For the solution Δr obtained in the kth iteration k , regard it as the search direction of this iteration and define a new value function Φ, use the Armijo line search method to solve an optimal subproblem so that the increment of r decreases with respect to Φ; Step 5: After the current iteration meets the predetermined conditions, check whether the virtual quantity weight σ needs to be updated. If it needs to be updated, return to step 2; if it does not need to be updated, proceed to the next step; Step 6: Check the model accuracy at the pseudo-spectrum matching points and other points outside the given matching points. If the model accuracy does not meet the requirements, a new sub-interval division is obtained through the grid adaptive update strategy, and return to step 2; otherwise, the iteration ends and proceeds to the next step; Step 7: Use the sampling controller to act on any moment in the time domain, restore the state quantity and corresponding control quantity at the current moment through the point information and pseudo-spectrum information, and use them as the output at the current moment.

2. The method according to claim 1, characterized in that: In step 2, we first establish the following Hamiltonian function: According to the first-order optimality condition Solve and get u=u(x,λ) and substitute it into A new Hamiltonian function form H = H (x, λ) is obtained which only contains the state quantity x and the co-state quantity λ.

3. The method according to claim 2, characterized in that: In step 3, specifically for the time domain [t s ,t f The jth subinterval among the Q subintervals divided by ] is defined as: have Continuous function x in subinterval (j) (t) and λ (j) (t)N (j) The LGL-type pseudospectral approximations of the orders are: The generating function in the subinterval is: Use N (j) The LGL-type pseudo-spectrum of order is approximated to obtain: in, is H(x,λ) at τ k The value at is the coefficient, is the differential matrix; For S (j) Find the variation and apply the first-order optimality condition: a and b are the variational lower and upper bounds, respectively; get: For x0 and x after pseudospectral approximation Q have: And for the LGL type pseudo-spectrum, we have:

4. The method according to claim 3, characterized in that: The first-order optimality condition for establishing the nonlinear equation system F(r)=0 in step 3 is specifically: When establishing the nonlinear equation system F(r)=0, first remember: Where m = 0, 1, 2, ..., N (j) , according to the above first-order optimality condition, we can derive: Taking partial derivative of f, we get: I is the identity matrix; remember: The nonlinear equation system F(r) is obtained: remember: Then the Jacobian matrix K of equation F with respect to r is: According to the characteristics of the optimal control problem that the initial state is given and the terminal is free, we get: Therefore, for a given point r k At , the obtained F(r) and K(r) are used to perform Newton's method iteration calculation: Δr k =-K -1 (r k )F(r k ) r k+1 =r k +Δr k 。 5. The method according to claim 4, characterized in that: In step 4, the value function Φ is specifically defined in the following form: Φ(r)=||F(r)||2 When using the Armijo line search method to solve the optimal subproblem, the k-th result Δr k As the search direction, find a suitable α that satisfies:

6. The method according to claim 5, characterized in that: In step 5, the following pre-set conditions are specifically set: Δr k ≤erro max or iter ≥ iter max Among them, erro max is the maximum allowable error, iter, iter max are the current number of iterations and the maximum number of iterations allowed respectively. The virtual quantity weight σ is updated once each time the above conditions are met.

7. The method according to claim 6, characterized in that: In step 6, the model accuracy is checked based on the following formula: Among them, e (j) is the model accuracy of the j-th subinterval, is the model error at any point selected within the subinterval, which is defined as follows: in, Calculated through pseudo-spectral transformation information; The model accuracy must meet the following requirements: e (j) ≤erro max ; If it is not satisfied, the grid adaptive update strategy is executed to obtain a new sub-interval division: First, define the model curvature as: For all points in the current subinterval And a given parameter R, perform the following update steps: ①If The subscripts max and ave represent the maximum value and the average value, respectively, dividing the subinterval into A new uniform subinterval; ②If Then the number of pseudo-spectral matching points in the subinterval is increased to 8. The method according to claim 7, characterized in that: The process of restoring the state quantity at the current moment by using the point information and pseudo-spectrum information in step 7 includes: First, find the subinterval t∈[t j ,t j-1 ], and perform pseudo-spectral transformation to obtain τ: For the pseudo-spectral transformation parameters L and D corresponding to the sub-interval, define the sampling vector L τ , D τ , and the following calculation formula is obtained: L τ (k)=polyval([L(k,1);L(k,2);…;L(k,N (j) )],τ),k=0,1,2,…,N (j) D τ (k)=polyval([D(k,1);D(k,2);…;D(k,N (j) )],τ),k=0,1,2,…,N (j) Among them, polyval(p,τ) means calculating the value of polynomial p at τ, p is a vector of length n+1, and the elements are the coefficients of the n-degree polynomial in descending order; According to the properties of the pseudo-spectrum, we can obtain:

Images