Rigid / non-rigid algorithm adaptive switching method for solving system simulation model

Abstract

The invention discloses a rigid / non-rigid algorithm adaptive switching method for solving a system simulation model, and in the solving process of the system simulation model, the rigidity of an ODE-IVP equation is monitored in real time, and a proper solving method is automatically selected according to the current rigidity value. According to the method, in order to improve the solving efficiency as much as possible, the rigidity of the current equation is not calculated once in each solving step, and dynamic sampling inspection is carried out on the model according to a certain probability, so that the rigid change information of the model cannot be lost, and the calculation efficiency can be considered.

Description

Technical Field

[0001] The present invention belongs to the technical field of system simulation, and relates to a method for adaptively switching between a stiff / non-stiff algorithm for solving a system simulation model. Background Art

[0002] A system simulation model can be expressed as an initial value problem of an N-dimensional ordinary differential equation in the real number field, that is, ODE-IVP (Ordinary Differential Equation - Initial Value Problem):

[0003]

[0004] For ODE-IVP, the stiffness of its equation is an important index. If there are simultaneously solution components that change extremely fast (with a small time scale) and extremely slow (with a large time scale) in the solution components of ODE-IVP, this characteristic is mathematically called stiffness (Stiff), and the ODE-IVP equation describing this change process is called a stiff equation.

[0005] The stiffness of the equation does not depend on the numerical method for solving this ODE-IVP equation. However, precisely because of this property, the numerical solution of some large stiff equations is very difficult, and some implicit format methods have to be adopted. But the implicit format solution method requires numerical iteration, resulting in low computational efficiency and being difficult to meet the requirements for occasions with high real-time requirements for model solution.

[0006] In the field of system modeling and simulation, after engineers complete the model construction according to mathematical and physical principles, they cannot intuitively obtain the stiffness of the built model, so it is very difficult to select a suitable solution algorithm. Currently, the solution of system simulation models usually selects a certain solution method before the calculation starts, and then this method is used throughout the entire calculation process. If the model itself has little stiffness, but the engineer selects an implicit format algorithm that can solve stiff equations, this will waste unnecessary efficiency; if the model itself has large stiffness, but the engineer selects an explicit format algorithm that cannot solve stiff equations, the equation is very difficult to converge and often cannot be solved smoothly.

[0007] In addition, due to the complexity of the physical process, 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 applicable when a different set of parameter inputs is used. Summary of the Invention

[0008] The problem solved by the present invention is to propose a method for adaptive switching between rigid / non-rigid models in the field of system simulation. During the solution process of the system simulation model, the stiffness of the ODE-IVP equation can be monitored in real time, and a suitable solution method can be automatically selected according to the current stiffness value, ultimately improving the solution efficiency of the model.

[0009] The present invention is realized through the following technical solutions:

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

[0011] 1) Start solving the simulation model in the system simulation, and obtain the corresponding Jacobian matrix according to its ODE-IVP equation;

[0012] 2) Generate random numbers based on the probability density, and preliminarily judge the stiffness of the ODE-IVP equation based on the random number sampling rate. If not sampled, it is considered a non-rigid equation; if sampled, dynamically monitor the stiffness based on the cyclic memory window:

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

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

[0015] If not, it is considered a non-rigid equation; if so, it is a rigid equation, then calculate the Jacobian matrix J at the current moment, and calculate the stiffness value of the ODE-IVP equation;

[0016] 3) Record the stiffness value corresponding to the current solution step, where non-rigid equations are all recorded as M - 2, and rigid equations are the calculated stiffness values, and refresh the corresponding values in the cyclic memory window;

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

[0018] If the current stiffness value Stiff > M, use the implicit format method for solution;

[0019] If the previous stiffness value Stiff ≤ M, use the explicit format method for solution;

[0020] 5) Step forward in the simulation until the entire simulation process ends.

[0021] The operation of preliminarily judging the stiffness of the ODE-IVP equation based on the random number sampling rate is as follows:

[0022] Generate a random number x in [0,1] with the probability density function of uniform distribution;

[0023] If x ≤ 0.99, it means not drawn; otherwise, it means drawn.

[0024] The construction of the cyclic memory window is as follows:

[0025] The cyclic memory window records the stiffness values corresponding to the current solution step and the previous 4 solution steps;

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

[0027] Judge 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; if not, record the stiffness value as M - 2.

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

[0029] If the previous stiffness value Stiff ≤ M, the ODE45 method is used for solution.

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

[0031] The present invention provides a stiffness / non-stiffness model adaptive switching method applicable to the field of system simulation. During the solution process of the system simulation model, it can monitor the stiffness of the ODE-IVP equation in real time and automatically select an appropriate solution method according to the current stiffness value, significantly improving the solution efficiency of the system simulation model.

[0032] In order to improve the solution efficiency as much as possible, the present invention does not calculate the stiffness of the current equation at each solution step, but dynamically "spot-checks" the model according to a certain probability, which can not only ensure that the model stiffness change information is not lost, but also take into account the calculation efficiency.

[0033] The present invention is applicable to the solution of models with multiple disciplines cross-linked and coupled in multi-disciplinary simulation, especially in the simulation scenario where a large time-scale model (such as hydraulic, heat transfer, fuel, pneumatic, etc.) and a small time-scale model (such as motor, electronic control, etc.) are included in the same simulation project. BRIEF DESCRIPTION OF THE DRAWINGS

[0034] Figure 1 It is a schematic flow chart of the present invention;

[0035] Figure 2Schematic diagram for applying a cyclic memory window to record the stiffness values of the ODE-IVP equation for five consecutive solution steps;

[0036] Figure 3 Schematic diagram for the change of equation stiffness with the model;

[0037] Figure 4 Schematic diagram for the adaptive switching of the solution method with the equation stiffness. Detailed implementation manners

[0038] Aiming at the deficiencies of the solution methods for traditional system simulation models, the present invention focuses on the dynamic real-time monitoring of the stiffness of the ODE-IVP equation and the automatic selection of the solution method, and proposes an adaptive switching method for stiff / non-stiff models in the field of system simulation.

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

[0040] 1) Start solving the simulation model in the system simulation, and obtain the corresponding Jacobian matrix according to its ODE-IVP equation;

[0041] 2) Generate a random number based on the probability density, and preliminarily judge the stiffness of the ODE-IVP equation based on the random number sampling rate. If not sampled, it is considered a non-stiff equation; if sampled, dynamically monitor the stiffness based on a cyclic memory window:

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

[0043] If the average stiffness of the cyclic memory window corresponding to consecutive solution steps is less than the stiffness judgment threshold M, it is considered that the ODE-IVP equation of this solution step is a non-stiff equation; otherwise, judge whether the Jacobian matrix J of the current solution step has more than two eigenvalues with negative real parts;

[0044] If not, it is considered a non-stiff equation; if so, it is a stiff equation. Then calculate the Jacobian matrix J at the current moment and calculate the stiffness value of the ODE-IVP equation;

[0045] 3) Record the stiffness value corresponding to the current solution step, where non-stiff equations are all recorded as M - 2, and stiff equations are the calculated stiffness values, and refresh the corresponding values in the cyclic memory window;

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

[0047] If the current stiffness value Stiff > M, use the implicit format method for solution;

[0048] If the previous stiffness value Stiff≤M, the explicit format method is used to solve;

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

[0050] The improvements made by the present invention include:

[0051] First, the present invention proposes a "dotting 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 rigidity of the ODE-IVP equation can be obtained by calculating the ratio of the negative real parts of the maximum and minimum eigenvalues ​​of the state matrix (obtained by linearizing the ODE-IVP equation at a specific time). If the relevant calculation of the state matrix and its eigenvalues ​​is performed once in each solution step, the purpose of dynamically monitoring the rigidity of the ODE-IVP equation can be achieved, but this will bring a large amount of calculation and affect the overall solution efficiency. In order to take into account the overall solution efficiency and avoid missing any significant changes in the rigidity of the ODE-IVP equation as much as possible, the present invention proposes a rigidity dynamic monitoring method based on probability dotting.

[0053] 2. Based on rigid dynamic monitoring, the present invention proposes the selection of adaptive solution method

[0054] Select ODE45 (Dormand-Prince) and DISPRK22 (Diagonally Implicit Symplectic P-Runge-Kutta with 2 nd stage and 2 nd The 2nd-order diagonally implicit symplectic block Runge-Kutta method is used as the basic algorithm, and the automatic switching of explicit / implicit format algorithms is realized by identifying the rigidity of the ODE-IVP equation and the relevant states of the solution process.

[0055] Specific embodiments are given below to illustrate each part.

[0056] like Figure 1 As shown, the rigid / non-rigid algorithm adaptive switching method for solving the system simulation model proposed by the present invention includes the following operations:

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

[0058] Step 1.1: Start solving the simulation model in the system simulation and obtain the corresponding Jacobian matrix according to 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 [J] can be obtained according to its ODE-IVP equation dim×dim , that is:

[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 calculation time, and dim is the dimension of y;

[0063] For example: y = [y 1 , y 2 ,... y dim , y' = [y 1 ', y' 2 ,... y' dim ;

[0064] Step 1.2: Perform stiffness monitoring according to the "punching method"

[0065] If the stiffness of the ODE-IVP equation is monitored and recorded at each integration step, the following two problems will arise:

[0066] The first problem is that the computational amount increases; because obtaining the specific stiffness also consumes time, frequent calculations will inevitably slow down the overall solution efficiency;

[0067] The second problem: numerical fluctuations will cause violent fluctuations in stiffness, making it difficult to select a solution algorithm. Due to factors such as numerical errors, the monitored stiffness may be discontinuous and may fluctuate violently. In order to obtain a relatively stable stiffness value and avoid introducing untrue stiffness values to misjudge the algorithm switching, it is necessary to perform a certain "filtering" process on the stiffness obtained by real-time monitoring to eliminate interference "noise".

[0068] The present invention proposes the "punching method" to solve the above two problems, and the specific operation is as follows:

[0069] ① Generate a random number x in [0,1] with the probability density function of uniform distribution;

[0070] This random number x is not related to the ODE-IVP equation and the Jacobian matrix; the purpose of setting the random number x is that at the beginning of the calculation, regardless of the stiffness of the equation, it preferentially enters the explicit format solution (ODE45 algorithm). The explicit format calculation efficiency is generally relatively high. Even if it is not suitable, this algorithm is generally relatively stable. The implicit format solution algorithm is prone to non-convergence and reporting errors (occasionally);

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

[0072] Prioritize the explicit format (ODE45 algorithm) in this way to start the solution process first; if it doesn't work, then switch.

[0073] ③Construct a cyclic 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, as Figure 2 shown), and calculate the average stiffness value of these 5 solution steps. If the average stiffness value corresponding to 5 consecutive solution steps is less than 2000, it is considered that the ODE-IVP equation of this solution step is a non-stiff equation, and jump to ⑥; otherwise, jump to ④;

[0074] Here, the stiffness judgment threshold M is 2000, which can be changed according to specific requirements;

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

[0076] ⑤Calculate the Jacobian matrix J at the current moment, and calculate the stiffness of the ODE-IVP equation (corresponding to Figure 1 "Option 1" in it), that is, the ratio of the maximum value to the minimum value of the absolute values of the negative real parts of all eigenvalues of the Jacobian matrix J.

[0077] The calculation method of the stiffness of the ODE-IVP equation is:

[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 serial number with the maximum absolute value of the negative real part; j is the serial number with the minimum absolute value of the negative real part; λ represents the eigenvalue, λ k represents the kth eigenvalue, and Reλ is to take the real part of the current eigenvalue λ, equivalent to Re(λ);

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

[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 stiffness value and refresh the corresponding values in the cyclic memory window ( Figure 2 shown).

[0084] Step 2: Combine the current stiffness value, and judge and execute a single-step solution based on the adaptive solution algorithm selection criterion

[0085] Step 2.1: According to the stiffness value and relevant solver status information, clarify the criterion for selecting 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, it is considered that the current ODE-IVP equation is a typical stiff equation, and the DISPRK22 method is used for solution;

[0087] This method is a 2-stage 2nd-order symplectic-preserving RK method, called the 2-stage 2nd-order diagonal implicit symplectic-preserving block Runge-Kutta method;

[0088] If the previous stiffness value Stiff ≤ 2000, it is considered that the current ODE-IVP equation is a typical non-stiff equation, and the ODE45 method is used for solution. The ODE45 method is a Runge-Kutta method jointly implemented by the 4th order and the 5th order, generally called the "Dormand-Prince" method.

[0089] Step 2.2: Perform a step-by-step simulation until the entire simulation process ends.

[0090] Specific embodiments are given below.

[0091] DETEST is introduced as a test case for typical stiff equations. This test case is a classic one for examining the influence of equation stiffness variation on the calculation process. Specifically, DETEST is the following "Van der Pol" equation:

[0092]

[0093] The parameter η in this model fluctuates sinusoidally at a frequency of 1 Hz between 0 and 1000. The above test model is solved by applying ODE45, DISPRK22, and the method proposed by the present invention respectively. Each solution method is tested five times, and their respective time consumptions are statistically counted, as shown in Table 1.

[0094] Table 1 Summary of time consumption comparison of different solution methods (unit: second)

[0095] Number of tests Computation time of ODE45 Computation time of DISPRK22 Computation time of the present invention First time 6.42 3.83 2.72 Second time 6.85 3.39 2.65 Third time 6.82 3.31 2.82 Fourth time 6.98 3.1 2.87 Fifth time 6.95 3.24 2.98 Average 6.804 3.374 2.81

[0096] It can be seen from the statistical results that the stiffness / non-stiffness adaptive switching method proposed by the present invention has the shortest time consumption. For the same test model, the time consumption of the method proposed by the present invention for solution is 41.3% of that of using the ODE45 method alone and 83.3% of that of using the DISPRK22 method alone. The calculation efficiency is significantly better than that of using a single fixed solution method alone.

[0097] The schematic diagram of the equation stiffness changing with the model is as Figure 3As shown, it can be seen that the stiffness of the equation is relatively large at four moments: 132.1s, 198.1s, 264s, 329.1s, and 395s, and these moments exactly correspond to the positions where the calculation results change violently.

[0098] The adaptive switching process of the two algorithms with the stiffness of the equation is as Figure 4 shown. Among them, when the algorithm switches to 0 (i.e., false), it means using DISPRK22 (implicit / suitable for stiff equations), and when the algorithm switches to 1 (i.e., true), it means using ode45 (explicit / suitable for non-stiff equations). It can be seen from the curves in the figure that at the moments with relatively large stiffness, the method proposed in the present invention will switch the solution method from 1 to 0, that is, use the DISPRK22 method for solution.

[0099] The above-given embodiments are preferred examples for implementing the present invention, and the present invention is not limited to the above embodiments. Any non-essential addition or replacement made by those skilled in the art according to the technical features of the technical solution of the present invention shall fall within the protection scope of the present invention.

Claims

1. A rigid / non-rigid algorithm adaptive switching method for solving a system simulation model, characterized in that: The following operations are included: 1) Start solving the simulation model in the system simulation and obtain the corresponding Jacobian matrix according to 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. If no random number is selected, it is considered to be a non-rigid equation; if a random number is selected, the rigidity is dynamically monitored based on the cyclic memory window: Construct a loop memory window, record the stiffness values ​​corresponding to the current solution step and the previous solution steps, and calculate the stiffness average value of the solution steps in the loop memory window; If the average rigidity value of the loop memory window corresponding to multiple consecutive solution steps is less than the rigidity judgment threshold M, 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; If there is no such equation, it is considered a non-stiff equation; if there is such an equation, it is a stiff equation, then the Jacobian matrix J at the current moment is calculated, and the stiffness value of the ODE-IVP equation is calculated; 3) Record the rigidity value corresponding to the current solution step, where the non-rigid equations are all recorded as M-2, and the rigid equations are the calculated rigidity values, and refresh the corresponding values ​​in the loop memory window; 4) Combined with the current stiffness value, adaptively select the solution method and perform a single-step solution: If the current stiffness value Stiff>M, the implicit method is used for solution; If the previous stiffness value Stiff≤M, the explicit format method is used to solve; 5) Step the simulation until the entire simulation process is completed.

2. The method for adaptively switching between rigid and non-rigid algorithms for solving a system simulation model according to claim 1, characterized in that: The rigidity of the ODE-IVP equation is preliminarily judged based on the random number sampling rate. The operation is as follows: Generate a random number x in [0,1] using the average distribution as the probability density function; If x≤0.99, then the result is not a winner; otherwise, the result is a winner.

3. The method for adaptively switching between rigid and non-rigid algorithms for solving system simulation models according to claim 1, characterized in that: The construction of the cyclic memory window is: The loop memory window records the stiffness values ​​corresponding to the current solution step and the previous four solution steps; If the average rigidity value corresponding to five consecutive solution steps is less than M, the rigidity value is recorded as M-2; otherwise, the rigidity of the current solution step is determined as follows: Determine whether the Jacobian matrix J of the current solution step has more than two eigenvalues ​​with negative real parts. If so, calculate the rigidity value of the current solution step; if not, record the rigidity value as M-2.

4. The method for adaptively switching between rigid and non-rigid algorithms for solving a system simulation model according to claim 1, characterized in that: If the current stiffness value Stiff>M, the DISPRK22 method is used for solution; If the previous stiffness value Stiff≤M, the ODE45 method is used to solve.

Images