An Adaptive Switching Method for Rigid / Non-rigid Algorithms in Solving System Simulation Models

Abstract

This invention discloses an adaptive switching method for rigid / non-rigid algorithms in solving system simulation models. During the system simulation model solving process, the rigidity of the ODE-IVP equations is monitored in real time, and an appropriate solution method is automatically selected based on the current rigidity value. To maximize solution efficiency, this invention does not calculate the rigidity of the current equation at every solution step. Instead, it dynamically "samples" the model according to a certain probability, ensuring that information on model rigidity changes is not lost while maintaining computational efficiency.

Description

Technical Field

[0001] This invention belongs to the field of system simulation technology and relates to an adaptive switching method for rigid / non-rigid algorithms for solving system simulation models. Background Technology

[0002] System simulation models can all be expressed as initial value problems of an N-dimensional ordinary differential equation in the real number field, i.e., ODE-IVP (Ordinary Differential Equation-Initial Value Problem):

[0003]

[0004] For ODE-IVP, the stiffness of its equations is an important indicator. If the solution components of ODE-IVP simultaneously contain solution components that change very quickly (with a small time scale) and those that change very slowly (with a large time scale), this property is mathematically called stiffness, and the ODE-IVP equations describing this process are called stiff equations.

[0005] The stiffness of the equations does not depend on the numerical method used to solve the ODE-IVP equations. However, this very property makes numerical solutions to some large stiffness equations very difficult, necessitating the use of implicit schemes. But implicit scheme methods require numerical iteration, resulting in low computational efficiency, making them unsuitable for applications with high real-time requirements.

[0006] In the field of system modeling and simulation, engineers cannot intuitively obtain the rigidity of the model after constructing it based on mathematical and physical principles, making it difficult to select a suitable solution algorithm. Currently, solving system simulation models typically involves selecting a solution method before computation begins, and then using that method throughout the entire computation process. If the model itself is not very rigid, but the engineer chooses an implicit scheme algorithm capable of solving rigid equations, this results in unnecessary wasted efficiency; if the model itself is very rigid, but the engineer chooses an explicit scheme algorithm that cannot solve rigid equations, the equations are difficult to converge, often leading to unsuccessful solutions.

[0007] Furthermore, due to the complexity of physical processes, the stiffness of the same ODE-IVP equation often changes with various factors such as time, external excitation, internal state, and parameter selection. A solution method that is relatively suitable under a certain set of parameters may not be suitable when a different set of parameters is input. Summary of the Invention

[0008] The problem solved by this invention is to propose an adaptive switching method for rigid / non-rigid models applicable to the field of system simulation. During the solution process of the system simulation model, the rigidity of the ODE-IVP equation can be monitored in real time and an appropriate solution method can be automatically selected based on the current rigidity value, thereby improving the solution efficiency of the model.

[0009] This invention is achieved through the following technical solution:

[0010] An adaptive switching method for rigid / non-rigid algorithms in solving system simulation models includes the following operations:

[0011] 1) Begin solving the simulation model in the system simulation and obtain the corresponding Jacobian matrix based on its ODE-IVP equation;

[0012] 2) Random numbers are generated based on probability density. The rigidity of the ODE-IVP equation is initially determined based on the random number sampling rate. If no random number is selected, it is considered a non-rigid equation. If a random number is selected, the rigidity is dynamically monitored based on the cyclic memory window.

[0013] Construct a circular memory window to record the stiffness values ​​corresponding to the current solution step and the previous few solution steps, and calculate the average stiffness value of the solution steps within the circular memory window;

[0014] If the average rigidity of the cyclic memory window corresponding to multiple consecutive solution steps is less than the rigidity judgment threshold M, then the ODE-IVP equation of the solution step is considered to be a non-rigid equation; otherwise, it is determined whether the Jacobian matrix J of the current solution step has more than two eigenvalues ​​with negative real parts.

[0015] If there is no equation, it is considered a non-rigid equation; if there is an equation, it is a rigid equation. Then, calculate the Jacobian matrix J at the current time and calculate the rigidity value of the ODE-IVP equation.

[0016] 3) Record the stiffness value corresponding to the current solution step. Non-rigid equations are all recorded as M-2, while rigid equations are recorded as the calculated stiffness value. Refresh the corresponding values ​​in the loop memory window.

[0017] 4) Based on the current stiffness value, adaptively select the solution method and perform a single-step solution:

[0018] If the current stiffness value Stiff > M, then the implicit scheme method is used for solving;

[0019] If the initial stiffness value Stiff ≤ M, then the explicit solution method is used.

[0020] 5) Step-by-step simulation until the entire simulation process is completed.

[0021] The method for initially determining the rigidity of the ODE-IVP equation based on the random number sampling rate is as follows:

[0022] Generate a random number x in the range [0,1] using the average distribution as the probability density function;

[0023] If x ≤ 0.99, then no one was selected; otherwise, they were selected.

[0024] The circular memory window is constructed as follows:

[0025] The circular memory window records the stiffness values ​​corresponding to the current solution step and the previous four solution steps;

[0026] If the average stiffness value corresponding to 5 consecutive solution steps is less than M, then the stiffness value is recorded as M-2; otherwise, the stiffness of the current solution step is determined:

[0027] Determine whether the Jacobian matrix J of the current solution step has more than two eigenvalues ​​with negative real parts. If so, calculate the stiffness value of the current solution step; otherwise, denote the stiffness value as M-2.

[0028] Furthermore, if the current stiffness value Stiff > M, then the DISPRK22 method is used for solving;

[0029] If the initial stiffness value Stiff ≤ M, then the ODE45 method is used to solve it.

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

[0031] The proposed method for adaptive switching between rigid and non-rigid models in system simulation can monitor the rigidity of the ODE-IVP equations in real time and automatically select the appropriate solution method based on the current rigidity value during the solution process of the system simulation model, which significantly improves the solution efficiency of the system simulation model.

[0032] In order to maximize the solution efficiency, this invention does not calculate the stiffness of the current equation in every solution step. Instead, it dynamically "samples" the model according to a certain probability, which can ensure that the stiffness change information of the model is not lost while also taking into account the computational efficiency.

[0033] This invention is applicable to solving models in multidisciplinary simulations where multiple disciplines are interconnected and coupled, especially simulation scenarios in the same simulation project that simultaneously include large-scale models (hydraulic, heat transfer, fuel, pneumatic, etc.) and small-scale models (motor, electronic control, etc.). Attached Figure Description

[0034] Figure 1 This is a schematic diagram of the process of the present invention;

[0035] Figure 2A schematic diagram illustrating the application of a cyclic memory window to record the stiffness values ​​of the ODE-IVP equation over five consecutive solution steps;

[0036] Figure 3 This is a schematic diagram illustrating how the stiffness of the equation varies with the model.

[0037] Figure 4 This is a schematic diagram illustrating the adaptive switching of the solution method as the equation becomes rigid. Detailed Implementation

[0038] To address the shortcomings of traditional system simulation model solution methods, this invention focuses on dynamically and in real-time monitoring the stiffness of the ODE-IVP equations and automatically selecting the solution method, proposing an adaptive switching method for rigid / non-rigid models applicable to the field of system simulation.

[0039] An adaptive switching method for rigid / non-rigid algorithms in solving system simulation models includes the following operations:

[0040] 1) Begin solving the simulation model in the system simulation and obtain the corresponding Jacobian matrix based on its ODE-IVP equation;

[0041] 2) Random numbers are generated based on probability density. The rigidity of the ODE-IVP equation is initially determined based on the random number sampling rate. If no random number is selected, it is considered a non-rigid equation. If a random number is selected, the rigidity is dynamically monitored based on the cyclic memory window.

[0042] Construct a circular memory window to record the stiffness values ​​corresponding to the current solution step and the previous few solution steps, and calculate the average stiffness value of the solution steps within the circular memory window;

[0043] If the average rigidity of the cyclic memory window corresponding to multiple consecutive solution steps is less than the rigidity judgment threshold M, then the ODE-IVP equation of the solution step is considered to be a non-rigid equation; otherwise, it is determined whether the Jacobian matrix J of the current solution step has more than two eigenvalues ​​with negative real parts.

[0044] If there is no equation, it is considered a non-rigid equation; if there is an equation, it is a rigid equation. Then, calculate the Jacobian matrix J at the current time and calculate the rigidity value of the ODE-IVP equation.

[0045] 3) Record the stiffness value corresponding to the current solution step. Non-rigid equations are all recorded as M-2, while rigid equations are recorded as the calculated stiffness value. Refresh the corresponding values ​​in the loop memory window.

[0046] 4) Based on the current stiffness value, adaptively select the solution method and perform a single-step solution:

[0047] If the current stiffness value Stiff > M, then the implicit scheme method is used for solving;

[0048] If the initial stiffness value Stiff ≤ M, then the explicit solution method is used.

[0049] 5) Step-by-step simulation until the entire simulation process is completed.

[0050] The improvements made in this invention include:

[0051] I. This invention proposes a "spotting method" based on probability density sampling to dynamically monitor the rigidity of the ODE-IVP equation without affecting the solution efficiency as much as possible;

[0052] The stiffness of the ODE-IVP equations can be obtained by calculating the ratio of the negative real parts of the largest and smallest eigenvalues ​​of the state matrix (obtained by linearizing the ODE-IVP equations at a specific time step). While performing the calculation of the state matrix and its eigenvalues ​​at every solution step can achieve the goal of dynamically monitoring the stiffness of the ODE-IVP equations, this approach incurs a significant computational burden, impacting overall solution efficiency. To balance overall solution efficiency with minimizing the omission of any significant changes in the stiffness of the ODE-IVP equations, this invention proposes a dynamic stiffness monitoring method based on probability-based point marking.

[0053] II. Based on rigid dynamic monitoring, this invention proposes an adaptive solution method selection.

[0054] ODE45 (Dormand-Prince) and DISPRK22 (Diagonally Implicit Symplectic P-Runge-Kutta with 2) were selected. nd stage and 2 nd The algorithm is based on the second-order, second-level, diagonal implicit symplectic block Runge-Kutta method. It achieves automatic switching between explicit and implicit scheme algorithms by identifying the rigidity of the ODE-IVP equation and the relevant states in the solution process.

[0055] The following specific embodiments illustrate each part.

[0056] like Figure 1 As shown, the adaptive switching method for rigid / non-rigid algorithms for solving system simulation models proposed in this invention includes the following operations:

[0057] Step 1: Dynamically monitor the current stiffness of the ODE-IVP equation

[0058] Step 1.1: Start solving the simulation model in the system simulation and obtain the corresponding Jacobian matrix based on its ODE-IVP equation;

[0059] For example, when solving a system simulation model, the Jacobian matrix of f(t,y) in y'=f(t,y) is calculated.

[0060] The corresponding Jacobian matrix can be obtained from its ODE-IVP equation [J]. dim×dim ,Right now:

[0061]

[0062] Where y is the solution vector of the current ODE-IVP equation, y' is the differential of the solution vector of the current ODE-IVP equation, t is the current computation time, and dim is the dimension of y.

[0063] For example: y = [y1, y2, ... y dim ],y'=[y1',y'2,…y' dim ];

[0064] Step 1.2: Conduct rigidity monitoring based on the "point-marking method"

[0065] Monitoring and recording the stiffness of the ODE-IVP equations at every integration step would lead to the following two problems:

[0066] The first problem is the increased computational load; since obtaining the specific rigidity also takes time, frequent calculations will inevitably slow down the overall solution efficiency.

[0067] The second problem is that numerical fluctuations can lead to drastic fluctuations in stiffness, making it difficult to select the solution algorithm. Due to factors such as numerical errors, the monitored stiffness may be discontinuous and fluctuate drastically. In order to obtain relatively stable stiffness values ​​and avoid introducing inaccurate stiffness values ​​that could cause misjudgments in algorithm switching, it is necessary to perform a certain "filtering" process on the stiffness obtained from real-time monitoring to remove interference "noise".

[0068] This invention proposes a "dot method" to solve the above two problems. The specific operation is as follows:

[0069] ① Generate a random number x in [0,1] using the average distribution as the probability density function;

[0070] The random number x is unrelated to the ODE-IVP equation and the Jacobian matrix. The purpose of setting the random number x is to prioritize the explicit scheme solution (ODE45 algorithm) at the beginning of the calculation, regardless of the rigidity of the equation. The explicit scheme calculation is generally more efficient. Even if it is not suitable, this algorithm is generally more stable. The implicit scheme solution algorithm is prone to non-convergence and error (occasionally).

[0071] ② If x ≤ 0.99, jump directly to ⑥; if x > 0.99, jump to ③;

[0072] This prioritizes the explicit solution format (ODE45 algorithm) to start the solution process; if that fails, then switch to another format.

[0073] ③ Construct a circular memory window to record the stiffness values ​​corresponding to the current solution step and the previous 4 solution steps (a total of 5 solution steps' stiffness values, such as...). Figure 2 (as shown), and calculate the stiffness average of these 5 solution steps. If the stiffness average of the 5 consecutive solution steps is less than 2000, then the ODE-IVP equation of that solution step is considered to be a non-rigid equation, and jump to ⑥; otherwise, jump to ④;

[0074] Here, the rigid judgment threshold M is set to 2000, which can be changed according to specific needs;

[0075] ④ Determine whether the Jacobian matrix J at the current moment has more than two eigenvalues ​​with negative real parts. If so, it is a rigid equation, and jump to ⑤; otherwise, it is considered a non-rigid equation, and jump to ⑥.

[0076] ⑤ Calculate the Jacobian matrix J at the current time step, and calculate the stiffness (corresponding to) the ODE-IVP equations. Figure 1 The "Choice 1" in the equation refers to the ratio of the maximum to the minimum value of the negative real parts of all Jacobian matrices J after taking the absolute value of their eigenvalues.

[0077] The method for calculating the stiffness of the ODE-IVP equations is as follows:

[0078]

[0079] Where Stiff is the stiffness value of the ODE-IVP equation; m is the number of eigenvalues ​​with negative real parts;

[0080] K is the index of the eigenvalue with the largest negative real part; j is the index of the eigenvalue with the smallest negative real part; λ represents the eigenvalue. k Represents the k-th eigenvalue. Reλ is the real part of the current eigenvalue λ, which is equivalent to Re(λ).

[0081] ⑥ Directly record the stiffness value as 1998 (corresponding to Figure 1 Option 2).

[0082] Step 1.3: Calculate and record the stiffness value corresponding to the current solution step according to the selected stiffness monitoring method.

[0083] Record the current rigidity value and refresh the cyclic memory window. Figure 2 (As shown) the corresponding values.

[0084] Step 2: Based on the current stiffness value, determine and execute the single-step solution according to the selection criteria of the adaptive solution algorithm.

[0085] Step 2.1: Based on the stiffness value and relevant solver state information, determine the selection criteria for the solution algorithm, perform a single-step simulation according to the algorithm selection result, and record the result data:

[0086] If the current stiffness value Stiff > 2000, then the current ODE-IVP equation is considered to be a typical stiff equation, and the DISPRK22 method is used to solve it.

[0087] This method is a level 2, order 2 symplectic RK method, also known as the level 2, order 2 diagonal implicit symplectic block Runge-Kutta method.

[0088] If the previous stiffness value Stiff ≤ 2000, then the current ODE-IVP equations are considered to be typical non-rigid equations, and the ODE45 method is used to solve them. The ODE45 method is a Runge-Kutta method that is implemented in both 4th and 5th order, and is generally referred to as the "Dormand-Prince" method.

[0089] Step 2.2: Step-by-step simulation until the entire simulation process is completed.

[0090] Specific implementation examples are given below.

[0091] DETEST is introduced as a typical test case for stiff equations. This test case is a classic example for examining the impact of changes in the stiffness of the equations on the computational process. Specifically, DETEST is the following "van der Bohr" equation:

[0092]

[0093] The parameter η in the model fluctuates sinusoidally at a frequency of 1Hz between 0 and 1000. The above test model was solved using ODE45, DISPRK22, and the method proposed in this invention, respectively. Each solution method was tested five times, and the time taken was recorded, as shown in Table 1.

[0094] Table 1. Summary of time consumption for different solution methods (unit: seconds)

[0095] Number of tests ODE45 calculation time DISPRK22 computation time The calculation time of this invention first 6.42 3.83 2.72 The second 6.85 3.39 2.65 The third 6.82 3.31 2.82 Fourth 6.98 3.1 2.87 Fifth 6.95 3.24 2.98 average 6.804 3.374 2.81

[0096] The statistical results show that the rigid / non-rigid adaptive switching method proposed in this invention has the shortest solution time. For the same test model, the solution time of the method proposed in this invention is 41.3% of that of the ODE45 method alone and 83.3% of that of the DISPRK22 method alone. The computational efficiency is significantly better than that of using a single fixed solution method alone.

[0097] The stiffness of the equation varies with the model as shown in the diagram. Figure 3As shown, the equation exhibits significant rigidity at 132.1s, 198.1s, 264s, 329.1s, and 395s, which correspond precisely to the points where the calculation results change drastically.

[0098] The adaptive switching process of the two algorithms with the stiffness of the equation is as follows: Figure 4 As shown, switching the algorithm to 0 (false) indicates using DISPRK22 (implicit / suitable for rigid equations), while switching it to 1 (true) indicates using ode45 (explicit / suitable for non-rigid equations). The curves in the figure show that when the rigidity is high, the method proposed in this invention will switch the solution method from 1 to 0, i.e., use the DISPRK22 method for solving the problem.

[0099] The embodiments given above are preferred examples for implementing the present invention, and the present invention is not limited to the above embodiments. Any non-essential additions or substitutions made by those skilled in the art based on the technical features of the present invention are within the protection scope of the present invention.

Claims

1. A method for adaptive switching between rigid and non-rigid algorithms for solving system simulation models, characterized in that, Includes the following operations: 1) Begin solving the simulation model in the system simulation and obtain the corresponding Jacobian matrix based on its ODE-IVP equation; 2) Generate random numbers based on probability density, and preliminarily judge the rigidity of the ODE-IVP equation based on the random number sampling rate: Generate a random number x in [0,1] with the average distribution as the probability density function; This random number x is uncorrelated with the ODE-IVP equation and the Jacobian matrix; If x ≤ 0.99, then no one was selected. Otherwise, you will be selected. If it is not selected, it is considered to be a non-rigid equation. If selected, the rigidity is dynamically monitored based on the cyclic memory window: a cyclic memory window is constructed to record the rigidity values ​​corresponding to the current solution step and the previous few solution steps, and the average rigidity of the solution steps within the cyclic memory window is calculated; If the average rigidity of the recurrent memory window corresponding to multiple consecutive solution steps is less than the rigidity judgment threshold M If the ODE-IVP equations for the current solution step are not rigid, then the Jacobian matrix of the current solution step is considered to be non-rigid; otherwise, the Jacobian matrix of the current solution step is determined. J Does it have two or more eigenvalues ​​with negative real parts? If there is no equation, it is considered a non-rigid equation; if there is an equation, it is a rigid equation, and the Jacobian matrix at the current time step is calculated. J， And calculate the stiffness value of the ODE-IVP equation; 3) Record the stiffness value corresponding to the current solution step, where non-rigid equations are denoted as... M -2, the rigid equation then calculates the rigidity value and refreshes the corresponding value in the loop memory window; 4) Based on the current stiffness value, adaptively select the solution method and perform a single-step solution: If the current rigidity value Stiff > M If so, the implicit scheme method is used to solve it; If the front stiffness value Stiff ≤ M If so, the explicit scheme method should be used to solve it; 5) Step-by-step simulation until the entire simulation process is completed.

2. The adaptive switching method for rigid / non-rigid algorithms for solving system simulation models as described in claim 1, characterized in that, The circular memory window is constructed as follows: The circular memory window records the stiffness values ​​corresponding to the current solution step and the previous four solution steps; If the average stiffness value corresponding to 5 consecutive solution steps is less than M Then the stiffness value is denoted as M -2; otherwise, check the rigidity of the current solution step: Determine the Jacobian matrix of the current solution step J If there are two or more eigenvalues ​​with negative real parts, calculate the stiffness value for the current solution step; otherwise, denote the stiffness value as... M -2.

3. The adaptive switching method for rigid / non-rigid algorithms for solving system simulation models as described in claim 1, characterized in that, If the current rigidity value Stiff > M Then the DISPRK22 method is used to solve it; If the front stiffness value Stiff ≤ M If so, the ODE45 method is used to solve it.

Images