GPU-based robot simulation engine constraint solution calculation method and system

Abstract

The invention discloses a GPU (Graphics Processing Unit)-based robot simulation engine constraint solution calculation method and system. According to the method, firstly, a dynamic tree is converted into a sparse matrix form, and an inertia force term and a constraint acceleration term are calculated through parallel sparse matrix-vector multiplication and used for a subsequent iterative solution process. In each iteration, a solver firstly traverses a constraint set in parallel, judges a constraint type and calculates correlation quantities such as a constraint force, then calculates a target function gradient, solves a linear system based on a block inertia matrix divided by a dynamic tree to obtain a gradient in an inertia space, and then executes parallel line search to obtain an optimal forward step length. And updating the acceleration and the correlation vector for the next iteration. According to the method, the parallel computing characteristic of the GPU is fully utilized, the constraint solution acceleration method suitable for the GPU architecture is designed, and the computing efficiency of a robot simulation engine in a large-scale multi-instance scene is remarkably improved.

Description

TECHNICAL FIELD

[0001] The application belongs to the field of machine physical simulation, and particularly relates to a GPU-based robot simulation engine constraint solving calculation method and system. BACKGROUND

[0002] With the continuous improvement of robot intelligence and autonomous control level, robot simulation has become an important tool for algorithm verification, reinforcement learning training, and system optimization. By running a large number of virtual instances of robot and environment interaction in the simulation environment, the sample efficiency of policy training can be significantly improved, thereby shortening the debugging cycle and cost of real robots. Especially in the training scenario based on reinforcement learning or imitation learning, simultaneously running a large number of simulation instances (i.e., large-scale multi-instance parallel simulation) can generate more training data in the same time and speed up the policy convergence. Such computing tasks are naturally suitable for acceleration using the high parallelism of a graphics processing unit (GPU).

[0003] A traditional simulation engine is usually composed of multiple core modules, including: a collision detection module: responsible for detecting the contact situation between objects, calculating potential contact points, normal directions, and penetration depths; a constraint solving module: establishing and solving constraint equations according to contact information and joint constraint conditions, and correcting the calculated acceleration to meet physical constraints; a time integration module: updating the velocity and position state of the system according to the relationship between force and acceleration, and completing time step advancement.

[0004] Among these modules, the constraint solving module is the computing core of physical simulation, which is responsible for handling a large number of rigid body joint constraints, contact constraints, and friction constraints in the robot system, ensuring the mechanical consistency and system stability during simulation. Typical constraint solving methods include projection Gauss-Seidel (PGS) iteration method, Newton method, conjugate gradient method, etc., the essence of which is multiple iteration solving and projection operation on sparse constraint matrices. This process requires a large number of floating point operations and memory accesses, accounting for the majority of the total simulation computation.

[0005] However, existing mainstream simulation engines (such as MuJoCo, Bullet, ODE, etc.) are mainly implemented based on CPU architecture, which has limited parallelism, insufficient memory bandwidth, and poor computational scalability. When the number of simulation instances increases, the time-consuming of constraint solving increases nearly linearly, making it difficult to meet the demand of large-scale real-time simulation. Although a few simulation engines based on GPU implementation (such as IsaacGym) have higher performance, they are not open source and are difficult to promote and research. Therefore, designing an efficient constraint solving parallel computing method for GPU architecture is of great significance to improve the efficiency of large-scale robot training and complex physical system simulation. SUMMARY

[0006] The purpose of this invention is to overcome the performance bottleneck of existing robot simulation engines in large-scale, multi-instance scenarios, and to propose a GPU-based constraint-solving computation method and system for robot simulation engines. This method fully utilizes the parallel computing architecture of GPUs, modularly reconstructing and parallelizing the constraint-solving process, significantly improving computational efficiency while ensuring simulation accuracy and stability. This method can be applied to scenarios such as robot reinforcement learning training, swarm robot simulation, and parallel testing of intelligent agents.

[0007] To achieve the above-mentioned objectives, the present invention specifically adopts the following technical solution: In a first aspect, the present invention provides a constraint solving calculation method for a GPU-based robot simulation engine, comprising the following steps: S1. In the GPU parallel computing environment, obtain the dynamic parameters of the robot in the current simulation scenario, namely the inertia matrix determined by the mass and moment of inertia of each rigid body of the robot, the constraint Jacobian matrix determined by the kinematic constraints and contact constraints of the robot, and the acceleration vector representing the motion state of the robot in the generalized coordinates. Calculate the product of the inertia matrix and the acceleration vector in parallel to obtain the inertial force term, and calculate the product of the constraint Jacobian matrix and the acceleration vector to obtain the constraint acceleration term. S2. Using the inertial force term and constraint acceleration term as the initial values ​​for iteration, the constraint is solved through multiple iterations. When the preset iteration stopping condition is reached, the final acceleration vector is output, and the constraint solution is completed. In each iteration, all active constraints in the current simulation time step are grouped according to their constraint type, and a thread is assigned to each constraint. Each thread calculates the constraint force required to satisfy the constraint in the constraint space in parallel according to its corresponding constraint type, and obtains the constraint force vector. Then, the product of the constraint force vector and the transposed constraint Jacobian matrix is ​​calculated in parallel to obtain the force in the generalized coordinate space, which is used to update the robot's motion state. S3. In the generalized coordinate space, calculate the objective function gradient of the constraint force vector based on the inertial force term, and use the inertial matrix to transform the objective function gradient into the gradient in the inertial space; S4. Obtain the gradient in the inertial space corresponding to the current simulation time step and use it as the current gradient. At the same time, obtain the search direction of the previous simulation time step and update it according to the conjugate gradient method to obtain the search direction vector of the current simulation time step. S5. Generate a series of candidate step sizes in the current search direction using line search, and then calculate the objective function value corresponding to each candidate step size in parallel. Calculate the step size corresponding to the minimum objective function value through parallel reduction operation as the optimal forward step size. Then update the acceleration, inertial force term and constraint acceleration term according to this step size for the next iteration.

[0008] Based on the above scheme, each step can be implemented in the following preferred manner.

[0009] As a preferred embodiment of the first aspect mentioned above, in S1, when calculating the product of the inertia matrix and the acceleration vector in parallel, the dynamic tree representation is pre-converted into a sparse matrix CSR representation, and a thread is allocated for each row. Each thread sequentially accesses the non-zero values ​​in each row to calculate the product.

[0010] As a preferred embodiment of the first aspect above, in S2, the iteration is terminated when the upper limit of the number of iterations is reached, or the change in the objective function is less than the preset iteration termination error threshold, or the absolute value of the gradient of the objective function is less than the preset iteration termination error threshold.

[0011] As a preferred embodiment of the first aspect mentioned above, in S3, the gradient of the objective function is obtained by subtracting the unconstrained force input to the system and the force acting in the generalized coordinate space from the inertial force term in sequence.

[0012] As a preferred embodiment of the first aspect mentioned above, in S3, when calculating the gradient in inertial space, a block-based solution mode is adopted. First, a system of linear equations is constructed with the inertial matrix as the coefficient matrix and the gradient of the objective function as the constant term vector. Then, the inertial matrix is ​​decomposed by Cholesky to obtain the decomposed matrix. A thread block is assigned to each diagonal block matrix in the decomposed matrix. Each thread block solves a subsystem. The number of threads in each thread block is the maximum value of the degrees of freedom of the corresponding subsystem. After all thread blocks are solved in parallel, the gradient components of each subsystem are obtained. The gradient components of all subsystems constitute the gradient in inertial space. Each subsystem is a system of linear equations corresponding to a diagonal block matrix.

[0013] As a preferred embodiment of the first aspect mentioned above, in S4, the search direction vector of the current simulation time step... The update method is as follows: in, This represents the gradient in inertial space corresponding to the current simulation time step; This represents the direction update coefficient for the current simulation time step in the conjugate gradient method; Indicates the first The search direction vector for each simulation time step; Indicates the first The search direction vector for each simulation time step.

[0014] As a preferred embodiment of the first aspect mentioned above, the specific process of S5 is as follows: S51. Calculate the product of the inertia matrix and the search direction vector of the current simulation time step and use it as the first correlation term; calculate the product of the constraint Jacobian matrix and the search direction vector of the current simulation time step and use it as the second correlation term. S52. Based on the preset step size sampling error threshold generate There are candidate step sizes with values ​​between 0 and 1, and two adjacent candidate step sizes satisfy... and , It is the number of all candidate step sizes; They represent the first One candidate step size; S53. Assign a thread to each candidate step size. Each thread calculates the objective function value corresponding to the candidate step size. Calculate in parallel to obtain the minimum value among all the objective function values ​​corresponding to the candidate step sizes. Take the candidate step size corresponding to the minimum value as the optimal step size. S54. Combine the optimal forward step length with the first... Multiply the search direction vector of the simulation time step by the first vector, and then multiply it by the second vector of the simulation time step. The accelerations at each simulation time step are summed to obtain the first... The acceleration at each simulation time step; multiply the optimal step length by the first correlation term, and then by the... The inertial force terms at each simulation time step are added together to obtain the first... The inertial force term for each simulation time step; multiply the optimal forward step length by the second correlation term, and then multiply by the first... The constraint acceleration terms at each simulation time step are summed to obtain the first... Constrained acceleration terms for each simulation time step.

[0015] As a preferred embodiment of the first aspect mentioned above, in S53, a quadratic approximation method is used to estimate the objective function value corresponding to the candidate step size. The specific calculation method is as follows: in, Indicates candidate step size The objective function estimate; Indicates the first The acceleration vector for each simulation time step; Indicates the first The search direction vector for each simulation time step; Indicates the first The objective function value for each simulation time step; Indicates the gradient operator; superscript Indicates matrix transpose; This represents the operator for extracting the Hessian matrix.

[0016] Secondly, the present invention provides a GPU-based robot simulation engine constraint solving calculation system, comprising: The data acquisition module is used to acquire the dynamic parameters of the target robot; The result acquisition module is used to perform constraint solving calculation based on the dynamic parameters of the target robot according to the robot simulation engine constraint solving calculation method described in any of the first aspects above, and output the acceleration vector of the target robot to update the motion state of the target robot.

[0017] Thirdly, the present invention provides a computer program product, including a computer program / instruction, which, when executed by a processor, can implement the constraint solving calculation method for a GPU-based robot simulation engine as described in any of the solutions in the first aspect above.

[0018] Compared with the prior art, the present invention has the following advantages: 1. GPU Parallel Solving Framework Design: Based on the characteristics of GPU architecture, a parallel solving process using the conjugate gradient method suitable for GPU was designed to achieve thread-level parallel solving of constraint equations.

[0019] 2. Inertia Matrix Calculation Optimization: A parallel multiplication algorithm based on sparse matrix structure is proposed, and the inertia matrix access mode is optimized to effectively reduce memory access latency.

[0020] 3. Block-based linear system solution strategy: A block-based inertia matrix partitioning method based on dynamic tree is adopted to achieve block-based parallel solution on GPU, thereby improving data locality and computational efficiency.

[0021] 4. Parallel line search mechanism: Design a parallel line search method suitable for GPU architecture, and accelerate the convergence of conjugate gradient iteration by calculating the optimal step size through discrete sampling and parallel reduction.

[0022] 5. Significant performance improvement: This method significantly improves the parallel efficiency of the constraint solving stage and reduces simulation frame latency in large-scale multi-instance simulation scenarios, providing reliable computational support for efficient robot training. Attached Figure Description

[0023] Figure 1 This is a flowchart of the method of the present invention; Figure 2 This is a system block diagram of the method of the present invention; Figure 3 This is a schematic diagram of 100 humanoid test simulation scenarios in this embodiment. Detailed Implementation

[0024] To make the above-mentioned objects, features, and advantages of the present invention more apparent and understandable, the specific embodiments of the present invention will be described in detail below with reference to the accompanying drawings. Many specific details are set forth in the following description to provide a thorough understanding of the present invention. However, the present invention can be practiced in many other ways different from those described herein, and those skilled in the art can make similar modifications without departing from the spirit of the present invention. Therefore, the present invention is not limited to the specific embodiments disclosed below. Technical features in the various embodiments of the present invention can be combined accordingly without mutual conflict.

[0025] In the description of this invention, it should be understood that the terms "first" and "second" are used only for descriptive purposes and should not be construed as indicating or implying relative importance or implicitly specifying the number of indicated technical features. Therefore, a feature defined with "first" and "second" may explicitly or implicitly include at least one of those features.

[0026] like Figure 1 As shown, in a preferred embodiment of the present invention, the conjugate gradient method is used as the constraint solving method. The constraint solving calculation method of the GPU-based robot simulation engine includes the following steps S1 to S5. The specific implementation process of each step is described in detail below.

[0027] S1. In a GPU parallel computing environment, obtain the robot's dynamic parameters in the current simulation scenario, which are the inertia matrices determined by the mass and moment of inertia of each rigid body of the robot. The constraint Jacobian matrix determined by the robot's kinematic and contact constraints. and the acceleration vector representing the robot's motion state in generalized coordinates. The product of the inertia matrix and the acceleration vector is calculated in parallel to obtain the inertial force term. Simultaneously, the product of the constraint Jacobian matrix and the acceleration vector is calculated to obtain the constraint acceleration term. .

[0028] It should be noted that in S1 of this invention, when calculating the product of the inertia matrix and the acceleration vector in parallel, the dynamic tree representation is converted into the CSR representation of the sparse matrix in advance, and a thread is allocated for each row. Each thread accesses the non-zero values ​​in each row in turn to calculate the product.

[0029] It should be noted that in robot dynamics, the inertia matrix contains the robot's mass distribution and inertial information, used to calculate the forces required for motion; it is a... Symmetric sparse matrix, This represents the total number of degrees of freedom in the simulation scene. The constraint Jacobian matrix, used in a constraint system, indicates how constraints affect the system's velocity or acceleration; its size is... This reflects the partial derivative of the constraint equation with respect to the generalized acceleration. This represents the number of system constraints. The acceleration vector can include linear acceleration and angular acceleration.

[0030] Traditional algorithms accumulate the contribution of each degree of freedom by traversing the dynamics tree from bottom to top, calculating the product of the inertia matrix and the acceleration vector. However, this bottom-up traversal approach is unsuitable for the parallel computing strategies of GPUs. To implement inertia matrix multiplication on a GPU, the dynamics tree needs to be pre-converted into a sparse matrix representation. A dynamics tree can be represented using a sparse matrix of length... An array representation of the dynamics tree, where each element records the index of the parent node for each degree of freedom, and the root node stores the element -1. The corresponding matrix representation of this dynamics tree is... The matrix is ​​a sparse matrix, where the column index of the non-zero element in each row corresponds to the preceding node of the corresponding degree of freedom. The preceding node of a degree of freedom is a node in the path from that node to the root node in the dynamics tree (excluding the node itself). In most cases, the inertia matrix is ​​a sparse matrix. In this embodiment of the invention, the sparse matrix specifically adopts the CSR (Compressed Sparse Row) format. In the parallel inertia matrix multiplication, a thread is allocated for each row, and each thread sequentially accesses the non-zero values ​​in each row to calculate the product. The calculation of the product of the constraint Jacobian matrix and the acceleration vector can be performed in parallel according to the construction form of the Jacobian matrix, which is existing technology and will not be elaborated further. The finally calculated inertial force term and constraint acceleration term together constitute the initial data basis for subsequent iterative solutions.

[0031] S2. Using the inertial force term and constraint acceleration term as the initial values ​​for iteration, the constraint is solved through multiple iterations. When the preset iteration stopping condition is reached, the final acceleration vector is output, and the constraint solution is completed. In each iteration, all active constraints within the current simulation time step are grouped according to their constraint type, and a thread is assigned to each constraint. Each thread calculates in parallel the constraint forces required to satisfy the constraint in the constraint space according to its corresponding constraint type, thus obtaining the constraint force vector. Then, the product of the constraint force vector and the transposed constraint Jacobian matrix is ​​calculated in parallel to obtain the forces in the generalized coordinate space, which are used to update the robot's motion state.

[0032] It should be noted that in S2, constraint solving is performed through iterative loops. The iteration terminates when the upper limit of the number of iterations is reached, or the change in the objective function is less than the preset iteration termination error threshold, or the absolute value of the gradient of the objective function is less than the preset iteration termination error threshold.

[0033] It's important to note that in S2, preprocessing is required before GPU computation begins to prepare data for parallel computing. This involves grouping all active constraints (such as contact constraints, hinge joint constraints, ball joint constraints, etc.) within the current simulation time step according to their constraint type. All data related to the same type of constraint (such as geometric parameters, friction coefficients, current violation values, etc.) is reorganized from their scattered state in memory and copied to a contiguous memory region in the GPU's global memory to construct a data structure suitable for GPU parallel access. This way, each constraint type corresponds to a contiguous data block, minimizing branching between threads. After preprocessing, each thread reads its corresponding portion of data from the contiguous global memory based on its constraint type, and then independently calls the dedicated solver for that constraint type (e.g., calculating a normal impulse or a tangential friction force). Based on the corresponding physical formula, it outputs a constraint force scalar value (or a small force vector). Thus, the constraint force scalar values ​​calculated by all threads are automatically written into the constraint force vector in global memory. The corresponding position. After all constraint force calculations are completed, another highly optimized parallel computing kernel (such as CUDA-based sparse matrix-vector multiplication) is launched to further parallelize the calculation of the product of the constraint force vector and the transposed constraint Jacobian matrix. , where superscript This represents the transpose of the matrix. Since the constraint Jacobian matrix itself is sparse, and its structure is determined by the constraint type and system topology, techniques such as parallel reduction can be used to efficiently calculate the equivalent forces acting on the entire robot multibody system in the generalized coordinate space.

[0034] S3. In generalized coordinate space, according to the inertial force term Calculate the gradient of the objective function of the constraint force vector, and use the inertia matrix to transform the gradient of the objective function into the gradient in inertial space.

[0035] It should be noted that in this invention, S3 first calculates the gradient of the objective function. The actual calculation formula depends on the chosen objective function. In this embodiment, a specific method for calculating the gradient of the objective function is provided for reference: It is obtained by successively subtracting the unconstrained force input to the system and the force acting in the generalized coordinate space from the inertial force term.

[0036] in, It is an unconstrained force input by the system; It is a force acting in a generalized coordinate space.

[0037] The parallel computation of the above formula belongs to the parallel computation of basic vector operations, which is existing technology, so it will not be elaborated further.

[0038] After obtaining the gradient of the objective function, it's important to consider that it cannot be directly used for efficient updating of the constraint forces. To find the most effective update direction, the gradient needs to be transformed into "inertial space," which requires calculating the gradient in inertial space based on the gradient of the objective function. The gradient direction represents the constraint adjustment direction that can most effectively reduce constraint violations.

[0039] It should be noted that in S3 of this invention, the gradient in the inertial space is calculated. At that time, a block-based solution approach was adopted. First, a calculation based on the inertia matrix was constructed. As the coefficient matrix, it is derived from the gradient of the objective function. Linear equations as constant term vectors Then, the inertia matrix is ​​decomposed by Cholesky to obtain the decomposed matrix. A thread block is assigned to each diagonal block matrix in the decomposed matrix. Each thread block solves a subsystem. The number of threads in each thread block is the maximum number of degrees of freedom of the corresponding subsystem. After all thread blocks are solved in parallel, the gradient components of each subsystem are obtained. The gradient components of all subsystems constitute the gradient in inertial space. Each subsystem is a system of linear equations corresponding to a diagonal block matrix.

[0040] In this embodiment, the specific process of obtaining the gradient in inertial space is as follows: Calculate the gradient in inertial space This is equivalent to solving the following linear system, where the gradient is the solution vector of the linear system: In large-scale multi-instance simulation scenarios, the dynamics tree of each instance is independent, resulting in a diagonally partitioned structure for the inertia matrix, as shown in the following equation. Each diagonal block matrix corresponds to the inertia matrix representation of a single instance: in, This represents the number of independent dynamical trees, i.e., the number of instances of the system. , , They represent the first The, the The, the A diagonal block matrix; When used as a subscript, it represents the subsystem index divided according to the dynamic tree, and its value is consistent with the diagonal block matrix index, that is, it can also represent the diagonal block matrix index.

[0041] To solve the above linear system, the inertia matrix needs to be decomposed using Cholesky decomposition beforehand. The decomposed matrix still consists of diagonal block matrices, the first... A diagonal block matrix can be represented in the following form: in, The first term obtained from Cholesky decomposition is... A lower triangular matrix; The first term obtained from Cholesky decomposition is... A diagonal matrix.

[0042] In the parallel implementation of the solution, a thread block is allocated for each diagonal block matrix. Each thread block solves a subsystem, and the number of threads in each thread block is the maximum number of degrees of freedom of all subsystems. Here, a subsystem is a system of linear equations corresponding to a diagonal block matrix. Taking the diagonal block matrix as an example, its corresponding system of linear equations (i.e., the first...) The subsystems are shown in the following equation: in, For the first Gradient components of each subsystem; The first in the gradient of the objective function Each component.

[0043] It should also be noted that there are already relatively mature methods for solving the subsystems in existing technologies. Solving the subsystems involves solving an upper triangular matrix and a lower triangular matrix. Specifically, the solution is performed row by row: for the upper triangular matrix, it is solved in reverse order from the last row; for the lower triangular matrix, it is solved in ascending order from the first row. Finally, the solutions for all subsystems can be combined to form the complete gradient in inertial space.

[0044] In this embodiment, the solution process of the subsystem is briefly described below. Solving the subsystem involves solving an upper triangular matrix, a lower triangular matrix, and a diagonal matrix. The specific solution sequence is as follows: 1) Solve the upper triangular system of equations ; 2) Solve the diagonal system of equations ; 3) Solve the lower trigonometric equation system ; in, These are all intermediate variables used in the solution process.

[0045] For upper triangular matrices, a back-substitution method is used to solve them in reverse order from the last row; for lower triangular matrices, a forward-substitution method is used to solve them in forward order from the first row. When solving each row, the thread corresponding to each non-zero element participates in the computation in parallel. When solving the upper and lower triangular equation systems, thread synchronization within the thread block is required after solving each row, ensuring that threads participating in solving the same subsystem are within the same thread block. For diagonal matrices, element-wise scalar division of the right-hand vector with respect to the diagonal element vectors is performed in parallel.

[0046] S4. Obtain the gradient in the inertial space corresponding to the current simulation time step and use it as the current gradient. At the same time, obtain the search direction of the previous simulation time step and update it according to the conjugate gradient method to obtain the search direction vector of the current simulation time step.

[0047] It should be noted that in S4, the search direction vector of the current simulation time step The update method is as follows: in, This represents the gradient in inertial space corresponding to the current simulation time step; This represents the direction update coefficient for the current simulation time step in the conjugate gradient method; Indicates the first The search direction vector for each simulation time step; Indicates the first The search direction vector for each simulation time step.

[0048] It should also be noted that the direction update coefficients in the conjugate gradient method There are various calculation formulas; in embodiment S4 of this invention, a commonly used Polak-Ribiere calculation formula is selected. The details are as follows: in, Indicates the first Gradient transpose of the objective function at each simulation time step; Indicates the first The gradient in inertial space corresponding to each simulation time step; Indicates the first Gradient transpose of the objective function at each simulation time step Indicates the first The simulation time step and the first The gradient difference in inertial space at each simulation time step z.

[0049] It should also be noted that in the first simulation time step of this invention, the search direction is initialized as a negative gradient in inertial space.

[0050] The above gives the update formula for the conjugate gradient method. The corresponding parallel computation involves the parallel computation of vector addition, subtraction and dot multiplication operations, which is a prior art.

[0051] S5. Generate a series of candidate step sizes in the current search direction using line search, and then calculate the objective function value corresponding to each candidate step size in parallel. Calculate the step size corresponding to the minimum objective function value through parallel reduction operation as the optimal forward step size. Then update the acceleration, inertial force term and constraint acceleration term according to this step size for the next iteration.

[0052] It should be noted that the specific process of S5 in this invention is as follows: S51. Calculate the inertia matrix Search direction vector at the current simulation time step The product of and as the first relevant term Calculate the constrained Jacobian matrix Search direction vector at the current simulation time step The product of and as the second relevant term .

[0053] S52. Based on the preset step size sampling error threshold generate Candidate step sizes with values ​​between 0 and 1 Two adjacent candidate step sizes satisfy and , It is the number of all candidate step sizes; They represent the first One candidate step size.

[0054] In this embodiment S52, for any Using a simple and practical uniform sampling method, candidate step sizes are generated according to the following formula, and the step sizes satisfy the error condition and order condition in the above formula: S53. Assign a thread to each candidate step size. Each thread calculates the objective function value corresponding to the candidate step size. Calculate in parallel the minimum objective function value among all candidate step sizes and take the candidate step size corresponding to the minimum value as the optimal forward step size.

[0055] In S53 of this invention, a quadratic approximation method is used to estimate the objective function value corresponding to the candidate step size. The specific calculation method is as follows: in, Indicates candidate step size The objective function estimate; Indicates the first The acceleration vector for each simulation time step; Indicates the first The search direction vector for each simulation time step; Indicates the first The objective function value for each simulation time step; Indicates the gradient operator; superscript Indicates matrix transpose; This represents the operator for extracting the Hessian matrix.

[0056] In this embodiment, the objective function value corresponding to the candidate step size can be expressed as follows: in, , as well as These are three values ​​that can be pre-calculated. The specific calculation formula depends on the selected objective function. An embodiment of this invention provides a reference calculation formula as follows: in, It is the mass matrix of constraints, which is a diagonal matrix, where each diagonal element represents the mass of the corresponding constraint. It is an unconstrained input from the system.

[0057] It should be noted that in S53 of this invention, the minimum value of the objective function corresponding to all candidate step sizes is obtained through parallel computation. This step can be computed using a parallel reduction operation, simply by defining the binary operation in the reduction operation as taking the minimum of two values. This is existing technology and will not be elaborated further. Therefore, this embodiment obtains the optimal step size, denoted as... : S54. Optimal forward step length With the Multiply the search direction vector of the simulation time step by the first vector, and then multiply it by the second vector of the simulation time step. Acceleration at each simulation time step Add them together to get the first one. Acceleration at each simulation time step ; Optimal forward length Related to the first item After multiplying, then with the first... Inertial force term at each simulation time step Add them together to get the first one. Inertial force term at each simulation time step ; Optimal forward length Related to the second item After multiplying, then with the first... Constraint acceleration terms for each simulation time step Add them together to get the first one. Constraint acceleration terms for each simulation time step .

[0058] In this embodiment S54, the acceleration, inertial force term, and constraint acceleration term are updated according to the optimal forward step length formula to obtain the first... Acceleration at each simulation time step Inertial force term and constrained acceleration terms : To better demonstrate the specific implementation and technical effects of the present invention, the GPU-based robot simulation engine constraint solving calculation method shown in steps S1 to S5 of the above preferred implementation is applied to a specific example.

[0059] Example The specific implementation process of the constraint solving calculation method based on the GPU robot simulation engine used in this embodiment is as described above and will not be repeated here.

[0060] To demonstrate the technical effects of this invention, this embodiment tested the constraint solution time for different CPU and GPU versions under different scenarios using the open-source robot simulation engine mujoco. Specific data is shown in Table 1. The test simulation scenarios used are as follows: Figure 3 As shown, the difference lies in the number of humanoids in the scene. Figure 3 The simulation scenario examples with 100 humanoids are given. The number of humanoids in the scenarios used for testing ranges from 1 to 200. The specific number of humanoids tested is shown in Table 1.

[0061] Table 1. Comparison of constraint solution time between CPU version and GPU version Therefore, it can be seen that the present invention can significantly improve the computational efficiency of constraint solving steps in complex scenarios with a large number of humanoid instances.

[0062] It should also be noted that the GPU-based robot simulation engine constraint solving calculation method in the above embodiments can essentially be executed by a computer program or module. Therefore, similarly, based on the same inventive concept, another preferred embodiment of the present invention also provides a GPU-based robot simulation engine constraint solving calculation system corresponding to the GPU-based robot simulation engine constraint solving calculation method provided in the above embodiments, such as... Figure 2 As shown, it includes: The data acquisition module is used to acquire the dynamic parameters of the target robot; The result acquisition module is used to perform constraint solving based on the dynamic parameters of the target robot according to the constraint solving calculation method of the robot simulation engine described in the above embodiment, and output the acceleration vector of the target robot to update the motion state of the target robot.

[0063] It should be noted that after obtaining the acceleration vector of the target robot, the joint velocities or overall velocity of the robot can be updated using numerical integration methods (such as the Euler method) based on the acceleration vector and the simulation time step. The updated velocities can then be used to further update the robot's position (including joint angles and end effector positions).

[0064] It is understood that the GPU-based robot simulation engine constraint solving calculation method described in S1-S5 above can essentially be implemented by a computer program. Therefore, based on the same inventive concept, another preferred embodiment of the present invention also provides a computer program product corresponding to the GPU-based robot simulation engine constraint solving calculation method provided in the above embodiments, which includes a computer program / instructions. When the computer program / instructions are executed by a processor, they can implement the GPU-based robot simulation engine constraint solving calculation method as described in the above embodiments.

[0065] It is understood that the processors mentioned above can be general-purpose processors, including central processing units (CPUs), network processors (NPs), etc.; they can also be digital signal processors (DSPs), application-specific integrated circuits (ASICs), field-programmable gate arrays (FPGAs), or other programmable logic devices, discrete gate or transistor logic devices, or discrete hardware components.

[0066] It should also be noted that those skilled in the art will understand that, for the sake of convenience and brevity, the specific working process of the system described above can be referred to the corresponding process in the foregoing method embodiments, and will not be repeated here. In the embodiments provided in this application, the division of steps or modules in the system and method is merely a logical functional division, and there may be other division methods in actual implementation. For example, multiple modules or steps may be combined or integrated together, and a module or step may also be split.

[0067] The embodiments described above are merely preferred embodiments of the present invention and are not intended to limit the invention. Those skilled in the art can make various changes and modifications without departing from the spirit and scope of the invention. Therefore, all technical solutions obtained through equivalent substitution or transformation fall within the protection scope of the present invention.

Claims

1. A constraint solving calculation method based on a GPU-based robot simulation engine, characterized in that, Includes the following steps: S1. In the GPU parallel computing environment, obtain the dynamic parameters of the robot in the current simulation scenario, namely the inertia matrix determined by the mass and moment of inertia of each rigid body of the robot, the constraint Jacobian matrix determined by the kinematic constraints and contact constraints of the robot, and the acceleration vector representing the motion state of the robot in the generalized coordinates. Calculate the product of the inertia matrix and the acceleration vector in parallel to obtain the inertial force term, and calculate the product of the constraint Jacobian matrix and the acceleration vector to obtain the constraint acceleration term. S2. Using the inertial force term and constraint acceleration term as the initial values ​​for iteration, the constraint is solved through multiple iterations. When the preset iteration stopping condition is reached, the final acceleration vector is output, and the constraint solution is completed. In each iteration, all active constraints in the current simulation time step are grouped according to their constraint type, and a thread is assigned to each constraint. Each thread calculates the constraint force required to satisfy the constraint in the constraint space in parallel according to its corresponding constraint type, and obtains the constraint force vector. Then, the product of the constraint force vector and the transposed constraint Jacobian matrix is ​​calculated in parallel to obtain the force in the generalized coordinate space, which is used to update the robot's motion state. S3. In the generalized coordinate space, calculate the objective function gradient of the constraint force vector based on the inertial force term, and use the inertial matrix to transform the objective function gradient into the gradient in the inertial space; S4. Obtain the gradient in the inertial space corresponding to the current simulation time step and use it as the current gradient. At the same time, obtain the search direction of the previous simulation time step and update it according to the conjugate gradient method to obtain the search direction vector of the current simulation time step. S5. Generate a series of candidate step sizes in the current search direction using line search, and then calculate the objective function value corresponding to each candidate step size in parallel. Calculate the step size corresponding to the minimum objective function value through parallel reduction operation as the optimal forward step size. Then update the acceleration, inertial force term and constraint acceleration term according to this step size for the next iteration.

2. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 1, characterized in that, In S1, when calculating the product of the inertia matrix and the acceleration vector in parallel, the dynamic tree representation is first converted into a sparse matrix CSR representation, and a thread is allocated for each row. Each thread accesses the non-zero values ​​in each row in turn to calculate the product.

3. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 1, characterized in that, In S2, the iteration terminates when the upper limit of the number of iterations is reached, or the change in the objective function is less than the preset iteration termination error threshold, or the absolute value of the gradient of the objective function is less than the preset iteration termination error threshold.

4. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 1, characterized in that, In S3, the gradient of the objective function is obtained by subtracting the unconstrained force input to the system and the force acting in the generalized coordinate space from the inertial force term.

5. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 4, characterized in that, In S3, a block-based solution approach is used to calculate the gradient in inertial space. First, a system of linear equations is constructed, with the inertia matrix as the coefficient matrix and the gradient of the objective function as the constant term vector. Then, the inertia matrix is ​​decomposed using Cholesky decomposition to obtain the decomposed matrix. A thread block is assigned to each diagonal block of the decomposed matrix, and each thread block solves a subsystem. The number of threads in each thread block is the maximum number of degrees of freedom of the corresponding subsystem. After all thread blocks are solved in parallel, the gradient components of each subsystem are obtained. The gradient components of all subsystems constitute the gradient in inertial space. Each subsystem is a system of linear equations corresponding to a diagonal block matrix.

6. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 1, characterized in that, In S4, the search direction vector of the current simulation time step The update method is as follows: ； in, This represents the gradient in inertial space corresponding to the current simulation time step; This represents the direction update coefficient for the current simulation time step in the conjugate gradient method; Indicates the first The search direction vector for each simulation time step; Indicates the first The search direction vector for each simulation time step.

7. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 1, characterized in that, The specific process of S5 is as follows: S51. Calculate the product of the inertia matrix and the search direction vector of the current simulation time step and use it as the first correlation term; calculate the product of the constraint Jacobian matrix and the search direction vector of the current simulation time step and use it as the second correlation term. S52. Based on the preset step size sampling error threshold generate There are candidate step sizes with values ​​between 0 and 1, and two adjacent candidate step sizes satisfy... and , It is the number of all candidate step sizes; They represent the first One candidate step size; S53. Assign a thread to each candidate step size. Each thread calculates the objective function value corresponding to the candidate step size. Calculate in parallel to obtain the minimum value among all the objective function values ​​corresponding to the candidate step sizes. Take the candidate step size corresponding to the minimum value as the optimal step size. S54. Combine the optimal forward step length with the first... Multiply the search direction vector of the simulation time step by the first vector, and then multiply it by the second vector of the simulation time step. The accelerations at each simulation time step are summed to obtain the first... The acceleration at each simulation time step; multiply the optimal step length by the first correlation term, and then by the... The inertial force terms at each simulation time step are added together to obtain the first... The inertial force term for each simulation time step; multiply the optimal forward step length by the second correlation term, and then multiply by the first... The constraint acceleration terms at each simulation time step are summed to obtain the first... Constrained acceleration terms for each simulation time step.

8. The constraint solving calculation method for a GPU-based robot simulation engine as described in claim 7, characterized in that, In S53, a quadratic approximation method is used to estimate the objective function value corresponding to the candidate step size. The specific calculation method is as follows: ； in, Indicates candidate step size The objective function estimate; Indicates the first The acceleration vector for each simulation time step; Indicates the first The search direction vector for each simulation time step; Indicates the first The objective function value for each simulation time step; Indicates the gradient operator; superscript Indicates matrix transpose; This represents the operator for extracting the Hessian matrix.

9. A GPU-based robot simulation engine constraint solving calculation system, characterized in that, include: The data acquisition module is used to acquire the dynamic parameters of the target robot; The result acquisition module is used to perform constraint solving calculation based on the dynamic parameters of the target robot according to the robot simulation engine constraint solving calculation method according to any one of claims 1 to 8, and output the acceleration vector of the target robot for updating the motion state of the target robot.

10. A computer program product comprising a computer program / instructions, characterized in that, When the computer program / instruction is executed by the processor, it can implement the constraint solving calculation method of the GPU-based robot simulation engine as described in any one of claims 1 to 8.

Images