Quadruped robot control method and device

Abstract

The invention discloses a quadruped robot control method and device.The method comprises the steps that quadruped robot joint angle information is obtained, and the quadruped robot joint angle information comprises hip joint adduction and abduction angles, hip joint flexion and extension angles and knee joint bending angles; calculating a rotation error vector according to the joint angle information of the quadruped robot; according to the rotation error vector, a model prediction control objective function is constructed, and the optimization objective of the model prediction control objective function is that the error between a prediction model and an expected state is minimum; according to the model prediction control objective function, quadratic programming solving processing is carried out, and a foot end force prediction value is obtained; generating a joint torque command by using a virtual force model control method; and robot control is conducted according to the foot end force predicted value and the joint torque command. The quadruped robot control method realizes quadruped robot control, improves environmental fitness and anti-rollover capability, and reduces calculation complexity and response delay. The method can be widely applied to the technical field of robot control.

Description

Technical Field

[0001] The present invention relates to the technical field of robot control, and particularly relates to a control method and device for a quadruped robot. Background Art

[0002] A quadruped robot is a four-legged bionic robot based on bionics and composed of four walking legs. The traditional zero moment point control method can make the horizontal projection of the robot's center of gravity always fall inside the support polygon formed by all feet. However, in a side slope environment, the stability margin decreases significantly, and the environmental adaptability is low. At the same time, a large number of floating-point operations need to be performed in each control cycle, with high computational complexity and large response delay. The traditional method based on a central pattern generator forms a stable oscillation pattern through the mutual inhibition and excitation between neurons, and can autonomously generate coordinated periodic rhythm signals to control movement without external rhythm input. However, the response time of the lateral force is long, and the response delay is large. At the same time, when the applied lateral force is large, the anti-interference torque threshold decreases significantly, and the anti-rollover ability is weak.

[0003] In summary, the technical problems existing in the related art need to be improved. Summary of the Invention

[0004] The embodiments of the present invention provide a control method and device for a quadruped robot, effectively improving the environmental adaptability and anti-rollover ability, and reducing the computational complexity and response delay.

[0005] On the one hand, the embodiments of the present invention provide a control method for a quadruped robot, including the following steps:

[0006] Obtain the joint angle information of the quadruped robot, where the joint angle information of the quadruped robot includes the hip adduction / abduction angle, the hip flexion / extension angle, and the knee bending angle;

[0007] Calculate the rotation error vector according to the joint angle information of the quadruped robot;

[0008] Construct a model predictive control objective function according to the rotation error vector, where the optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state;

[0009] Perform quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the foot end force;

[0010] Generate joint torque commands using the virtual force model control method;

[0011] Perform robot control according to the predicted value of the foot end force and the joint torque commands.

[0012] In some embodiments, calculating a rotation error vector based on the joint angle information of the quadruped robot includes:

[0013] Constructing a joint angle vector based on the joint angle information of the quadruped robot;

[0014] Calculating the end-effector spatial coordinates corresponding to each leg based on the joint angle vector;

[0015] Calculating an end-effector position combination matrix based on multiple end-effector spatial coordinates;

[0016] Calculating a correction matrix based on a correction coefficient and the end-effector position combination matrix;

[0017] Constructing a forward kinematics function corresponding to each leg based on the correction matrix;

[0018] Taking the partial derivative of the forward kinematics function with respect to the joint angle to obtain a Jacobian matrix;

[0019] Setting the yaw angle, pitch angle, and roll angle of the quadruped robot based on the Jacobian matrix;

[0020] Constructing a basic rotation matrix about the yaw angle based on the yaw angle;

[0021] Constructing a basic rotation matrix about the pitch angle based on the pitch angle;

[0022] Constructing a basic rotation matrix about the roll angle based on the roll angle;

[0023] Combining the basic rotation matrix about the yaw angle, the basic rotation matrix about the pitch angle, and the basic rotation matrix about the roll angle to obtain an overall rotation matrix;

[0024] Constructing an error matrix based on the overall rotation matrix and an expected rotation matrix;

[0025] Performing a logarithmic mapping process on the overall rotation matrix based on the error matrix;

[0026] Constructing the rotation error vector based on multiple overall rotation matrices after the logarithmic mapping process.

[0027] In some embodiments, constructing a model predictive control objective function based on the rotation error vector includes:

[0028] Constructing a translational motion model based on the total mass of the robot, the centroid acceleration, the foot end reaction force, and the gravitational acceleration;

[0029] Constructing a rotational motion model based on the inertia matrix of the robot about the centroid, the rotation error vector, the angular acceleration, the vector from the robot centroid to the leg contact point, and the foot end reaction force;

[0030] Construct a state vector based on the centroid position, centroid velocity, rigid body attitude parameters, and angular velocity;

[0031] Construct a translational part dynamic model based on the translational motion model and the state vector;

[0032] Construct a rotational part dynamic model based on the rotational motion model and the state vector;

[0033] Construct a state matrix based on the translational part dynamic model;

[0034] Construct a first input matrix based on the rotational part dynamic model;

[0035] Construct an overall continuous model based on the state matrix, the first input matrix, the state vector, and the control input vector;

[0036] Calculate the discretized state transition equation according to the discrete time step and the basic equation of the Euler integration method;

[0037] Construct a substitution model equation based on the overall continuous model and the discretized state transition equation;

[0038] Construct a discrete state transition matrix and a control influence matrix based on the substitution model equation;

[0039] Construct a discrete time model based on the discrete state transition matrix and the control influence matrix;

[0040] Perform an approximation process on the translational part according to the discrete time model to obtain a translational approximation equation set;

[0041] Perform an approximation process on the rotational part according to the discrete time model to obtain a rotational approximation equation set;

[0042] Construct the model predictive control objective function based on the translational approximation equation set, the rotational approximation equation set, the prediction model, and the desired state;

[0043] In some embodiments, the performing quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the end-effector force includes:

[0044] Construct a compression equation based on the end-effector reaction force;

[0045] Construct a standard quadratic programming objective function based on the compression equation and the model predictive control objective function;

[0046] Construct end-effector constraints, where the end-effector constraints include end-effector force limit constraints and friction constraints;

[0047] Combine the standard quadratic programming objective function and the foot-end constraints to obtain a standard quadratic programming problem;

[0048] Use a quadratic programming solver to solve the standard quadratic programming problem to obtain the predicted value of the foot-end force.

[0049] In some embodiments, constructing a compression equation according to the foot-end reaction force includes:

[0050] Set a control variable according to the foot-end reaction force;

[0051] Calculate a compensation acceleration according to the centroid position, centroid velocity, desired position, and desired velocity;

[0052] Calculate an angular acceleration according to the gain matrix of the attitude error, the gain matrix of the angular velocity error, the attitude error vector, the desired angular velocity, and the current angular velocity;

[0053] Use the compensation acceleration and the angular acceleration as state tracking targets;

[0054] Construct a friction cone constraint according to the foot-end force and the ground friction coefficient;

[0055] Perform a linear transformation on the friction cone constraint to obtain a system of linear inequalities, and use the system of linear inequalities as an energy consumption target;

[0056] According to the control variable, combine the state tracking target and the friction cone constraint to obtain a first programming problem, where the parameters of the first programming problem include the foot-end force, the weighted matrix of the state error, and the weighted matrix of the control input;

[0057] Construct a relationship equation between the robot state and the control input according to the first programming problem, the continuous-time system matrix, and the second input matrix;

[0058] Discretize the relationship equation between the robot state and the control input according to the second time step to obtain a discretized model;

[0059] Perform a compression transformation on the discretized model to obtain the compression equation.

[0060] In some embodiments, constructing a standard quadratic programming objective function according to the compression equation and the model predictive control objective function includes:

[0061] Rewrite the model predictive control objective function according to the quadratic penalty of the state error term and the control input term to obtain a first objective function;

[0062] Rewrite the first objective function according to the compression equation and the desired state to obtain a second objective function;

[0063] Expand the second objective function to obtain a standard quadratic form objective function;

[0064] Calculate the overall control quantity according to the inertia matrix of the robot base in the body coordinate system and the overall rotation matrix;

[0065] Construct a mapping matrix according to the sum of the end - effector forces of each leg, the rotational effects generated by the end - effector forces of each leg, and the skew - symmetric matrix function;

[0066] Rewrite the standard quadratic form objective function according to the overall control quantity and the mapping matrix to obtain a third objective function;

[0067] Expand the third objective function to obtain a fourth objective function;

[0068] Perform a form conversion on the fourth objective function to obtain the standard quadratic programming objective function.

[0069] In some embodiments, constructing the end - effector constraints includes:

[0070] Construct the end - effector force limit constraint according to the upper limit and lower limit of the end - effector force;

[0071] Construct a friction constraint matrix by using the friction pyramid approximation according to the fundamental matrix of friction;

[0072] Construct the friction constraint according to the friction constraint matrix, the upper bound of the vertical component, and the lower bound of the vertical component.

[0073] In some embodiments, using the virtual force model control method to generate joint torque commands includes:

[0074] In horizontal trajectory planning, perform interpolation according to the initial position, the target position, and the cubic Bezier curve interpolation equation to obtain a horizontal trajectory;

[0075] Differentiate the horizontal trajectory to obtain the horizontal velocity and horizontal acceleration;

[0076] In vertical trajectory planning, calculate the velocity and acceleration in the rising stage according to the initial coordinates, the target coordinates, the lifting height, the local time variable in the rising stage, and the cubic Bezier curve interpolation equation;

[0077] Calculate the velocity and acceleration in the descending stage according to the initial coordinates, the target coordinates, the lifting height, the local time variable in the descending stage, and the cubic Bezier curve interpolation equation;

[0078] Calculate the expected value of the swing trajectory based on the horizontal speed, horizontal acceleration, speed in the ascending stage, acceleration in the ascending stage, speed in the descending stage, and acceleration in the descending stage, where the expected value of the swing trajectory includes the expected foot end position and the expected foot end speed;

[0079] Calculate the corresponding foot end position of each leg in the robot coordinate system based on the center of mass position, the expected foot end position, and the fixed offset between the foot end and the geometric reference of the robot body;

[0080] Calculate the corresponding foot end speed of each leg in the robot coordinate system based on the center of mass speed and the expected foot end speed;

[0081] Calculate the current foot end speed based on the Jacobian matrix and the joint angular velocity;

[0082] Calculate the virtual control force based on the corresponding foot end position of each leg in the robot coordinate system, the corresponding foot end speed of each leg in the robot coordinate system, the current foot end speed, and the current foot end position;

[0083] Generate the joint torque command based on the virtual control force.

[0084] On the other hand, an embodiment of the present invention provides a quadruped robot control device, including:

[0085] A first module for obtaining the joint angle information of the quadruped robot, where the joint angle information of the quadruped robot includes the hip abduction / adduction angle, the hip flexion / extension angle, and the knee bending angle;

[0086] A second module for calculating the rotation error vector based on the joint angle information of the quadruped robot;

[0087] A third module for constructing a model predictive control objective function based on the rotation error vector, where the optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state;

[0088] A fourth module for performing quadratic programming solution processing based on the model predictive control objective function to obtain the predicted value of the foot end force;

[0089] A fifth module for generating the joint torque command using the virtual force model control method;

[0090] A sixth module for controlling the robot based on the predicted value of the foot end force and the joint torque command.

[0091] On the other hand, an embodiment of the present invention provides a computer device, including:

[0092] At least one processor;

[0093] At least one memory for storing at least one program;

[0094] When the at least one program is executed by the at least one processor, the at least one processor implements the method.

[0095] On the other hand, an embodiment of the present invention provides a computer-readable storage medium storing a computer program, and when the computer program is executed by a processor, the method is implemented.

[0096] The beneficial effects of the present invention are as follows:

[0097] In an embodiment of the present invention, first, the joint angle information of the quadruped robot is obtained, then, according to the joint angle information of the quadruped robot, the rotation error vector is calculated, according to the rotation error vector, the model predictive control objective function is constructed, and then, according to the model predictive control objective function, quadratic programming solution processing is performed to obtain the predicted value of the foot end force, the virtual force model control method is used to generate the joint torque command, and finally, according to the predicted value of the foot end force and the joint torque command, the robot is controlled, so that the control of the quadruped robot can be realized by combining model predictive control and virtual force model control, thereby improving the environmental adaptability and anti-rollover ability, and reducing the computational complexity and response delay.

[0098] Other features and advantages of the present invention will be described in the following specification, and, in part, will be obvious from the specification, or will be understood by implementing the present invention. The objectives and other advantages of the present invention can be realized and obtained by the structures specifically pointed out in the specification and the drawings. Description of the Drawings

[0099] In order to more clearly illustrate the technical solutions in the embodiments of the present application, the following will briefly introduce the drawings required for the description of the embodiments. Obviously, the following drawings are only some embodiments of the present application. For those of ordinary skill in the art, other drawings can be obtained based on these drawings without creative efforts.

[0100] Figure 1 It is a flowchart of a quadruped robot control method according to an embodiment of the present invention;

[0101] Figure 2 It is a schematic structural diagram of a quadruped robot control device according to an embodiment of the present invention;

[0102] Figure 3 It is a schematic hardware structure diagram of a computer device according to an embodiment of the present invention. Detailed Embodiments

[0103] In order to make the objectives, technical solutions, and advantages of this application more clearly understood, the following further elaborates on this application in conjunction with the accompanying drawings and embodiments. It should be understood that the specific embodiments described herein are only used to explain this application and are not used to limit this application. When the following description involves the accompanying drawings, unless otherwise indicated, the same numbers in different drawings represent the same or similar elements. The implementation manners described in the following exemplary embodiments do not represent all implementation manners consistent with the embodiments of this application. They are merely examples of devices and methods that are consistent with some aspects of the embodiments of this application as detailed in the appended claims.

[0104] It can be understood that the terms "first", "second", etc. used in this application can be used herein to describe various concepts, but unless otherwise specified, these concepts are not limited by these terms. These terms are only used to distinguish one concept from another. For example, without departing from the scope of the embodiments of this application, the first information can also be referred to as the second information, and similarly, the second information can also be referred to as the first information. Depending on the context, the words "if", "when" as used herein can be interpreted as "when...", "while...", or "in response to determining".

[0105] The terms "at least one", "multiple", "each", "any one", etc. used in this application, at least one includes one, two, or more than two, multiple includes two or more than two, each refers to each of the corresponding multiple, and any one refers to any one of the multiple.

[0106] Unless otherwise defined, all technical and scientific terms used herein have the same meaning as commonly understood by those skilled in the technical field to which this application belongs. The terms used herein are only for the purpose of describing the embodiments of this application and are not intended to limit this application.

[0107] Before elaborating on the embodiments of this application in detail, first, some nouns and terms involved in the embodiments of this application are explained. The nouns and terms involved in the embodiments of this application are applicable to the following explanations.

[0108] Zero moment point control (ZMP): It is a concept used for analyzing and controlling dynamic stability, especially in the gait planning of legged robots. The ZMP is the point of application of the ground reaction force such that the tipping moment at this point of the robot is zero. This means that in an ideal situation, the robot will not topple. By keeping the ZMP within the support polygon of the robot (usually the contact area of the feet), the robot can maintain stability when walking or running. ZMP control is widely applied in the gait control of bipedal robots such as ASIMO to ensure its balance during dynamic movements.

[0109] Central Pattern Generators (CPG): It is a biologically inspired control model that mimics the neural networks in organisms to generate rhythmic movements. CPG is typically used to generate continuous and periodic movement patterns, such as walking, swimming, and other forms of biological movement. The advantages of CPG lie in its inherent stability and robustness to external disturbances, enabling it to generate natural movements in complex environments. By adjusting the parameters of CPG, different movement patterns can be achieved, suitable for a variety of robotic applications.

[0110] Virtual model control (VMC): It is a control method that influences the physical movement of a robot by introducing virtual forces and torques. VMC uses a virtual model to simulate the desired dynamic behavior and converts it into actual control inputs. This method provides a flexible framework that can adjust the movement and posture of the robot according to task requirements. A key advantage of VMC is its intuitiveness, allowing designers to directly think about how to achieve complex motion control through the virtual model.

[0111] Model Predictive Control (MPC): It is an advanced control strategy that uses a mathematical model of the system to predict future behavior and selects appropriate control inputs through an optimization process. MPC solves an optimization problem in each control cycle to minimize a cost function while satisfying system constraints. The advantages of MPC are its ability to handle multi-input multi-output systems and effectively manage constraint conditions, making it widely used in fields such as industrial process control, autonomous vehicles, and robot control. By continuously updating the prediction model, MPC can adapt to the dynamic changes of the system.

[0112] Quadratic Programming (QP) optimization: It is an optimization method for solving a quadratic objective function while satisfying linear constraint conditions. In this system, QP optimization is used to distribute the foot-end forces of each support leg so that the overall resultant force and torque generated by these foot-end forces can approximate the desired control objective as much as possible while satisfying physical constraints (such as friction conditions and force limits).

[0113] Proportional-Derivative Control (PD control): It is a classical feedback control method. Its basic idea is to decompose the output of the controller into two parts, including a proportional (P) control term and a derivative (D) control term. In the proportional (P) control term, this part is proportional to the current error of the system, that is, u P = K P e, where e represents the error between the current state and the desired state, and KP is the proportional gain. Proportional control can quickly respond to the error and make the system state approach the target state. In the derivative (D) control term, this part is proportional to the change rate of the error, that is where is the derivative of the error, and K1 is the derivative gain. The derivative control plays a damping role and can suppress the oscillations and overshoots that may occur during the response process of the system, thereby improving the stability of the system.

[0114] In related technologies, a quadruped robot is a four-legged bionic robot based on bionics and composed of four walking legs. The traditional zero moment point control method can make the horizontal projection of the robot's center of gravity always fall within the support polygon formed by all feet. When the robot is in a balanced state, the reaction forces and torques exerted by each limb of the robot on the ground will make a certain point (i.e., the zero moment point) in a "zero-overturn state", that is, this point satisfies that the sum of the torques in the horizontal plane is zero, so as to maintain stable walking. In terms of sensors, ZMP control relies on six-axis force sensors on the soles of the feet (such as the force / torque data installed at each support point), and calculates the real-time ZMP position through weighted averaging. However, in a side-slope environment, the stability margin decreases significantly, and the environmental adaptability is low. When the terrain inclination angle exceeds 18° (the lateral critical angle), the area of the support polygon shrinks rapidly, and the stability margin decreases by more than 40%. Experiments show that in a 15° side-slope environment, the ZMP stability margin decreases by more than 38%. At the same time, a large number of floating-point operations need to be performed in each control cycle, the computational complexity is high, and the response delay is large. More than about 300 floating-point operations need to be performed in each control cycle, resulting in the real-time control frequency being limited to below 50Hz. The response delay is usually in the range of 80–150ms. Especially when an external force is suddenly applied (such as a lateral impact of 5N), the delay may reach 150ms. In addition, under complex bumpy terrains and external disturbances, the rollover critical inclination angle and the maximum anti-lateral force data under ZMP control are 18° and about 5N (flat ground state) respectively, while the lateral instantaneous force may be only 1.8–2.5N under dynamic conditions, and the anti-interference ability is weak.

[0115] Traditional central pattern generator (CPG)-based methods form stable oscillation patterns through the mutual inhibition and excitation between neurons, enabling the autonomous generation of coordinated periodic rhythm signals without external rhythm input and the adaptive adjustment of time-consuming models to control movements such as walking and swimming. Its core mechanisms include the pacemaker / follower mechanism and the recursive inhibition mechanism, forming stable oscillation patterns through the mutual inhibition and excitation between neurons. Commonly used CPG models include the Matsuoka oscillator, the Hopf oscillator, etc. These oscillators can be coupled with each other to form a CPG network, and different gait switches and adaptations can be achieved by adjusting parameters (such as the time constant τ, oscillation amplitude, coupling weight, etc.). Research on parameter optimization shows that on average, 17.3 iterative experiments are required to achieve effective parameter adjustment. Recent research combines deep reinforcement learning methods to optimize CPG oscillator parameters in real time, enabling the robot to achieve reliable adaptation in different terrains according to environmental feedback. However, the response time to lateral forces is long, and the response delay is large. The average response time to a sudden application of a 5N lateral force is about 200±25 ms, which is about 33% higher than the response delay of the zero moment point (ZMP). At the same time, when the applied lateral force is large, the anti-interference torque threshold drops significantly, and the anti-rollover ability is weak. When the applied lateral force exceeds 3.2N, the anti-interference torque threshold of the CPG drops significantly; under the action of a 5N lateral force, the threshold drops to 45% of the nominal value, and an irreversible instability phenomenon occurs at 6.5N. Generally, CPG control performs outstandingly in dynamic adaptation, but due to the time-consuming internal signal reconstruction and parameter adaptation, its response speed and stability still lag behind the requirements under extreme disturbance conditions.

[0116] In the intuitive control (VMC) method, intuitive control (VMC) is based on establishing the difference feedback between the target state and the actual state of the robot based on a virtual model, and generating control commands through real-time adjustment. Although intuitive control has the characteristic of fast response, it also has obvious defects: in gravel terrain, the foot-end trajectory error rate can reach 28%–35%, and when the support polygon area is less than 0.12 m 2 , the error rate is higher, exceeding 40%, and the decision error rate is relatively high; when encountering a 3.2N lateral force, the average lateral error of the foot-end trajectory of intuitive control is about 12.5 cm±2.3 cm, and the longitudinal error is about 8.7 cm±1.8 cm, and the influence of external disturbances is significant; when the obstacle density is higher than 3 per m 2 or the slope reaches 18°, intuitive control is prone to failure, the error rate is significantly higher than expected, and the recovery time is also about 1.2 seconds longer than that of CPG, reaching the limit in complex environments.

[0117] In view of this, in the embodiments of the present invention, under complex terrains and with external force interference for a quadruped robot, the robot is controlled by combining the support phase model predictive control (MPC) and the swing phase intuitive control (VMC), thereby improving the environmental adaptability and anti-rollover ability, and reducing the computational complexity and response delay.

[0118] A quadruped robot control method provided by an embodiment of the present application relates to the technical field of robot control. The quadruped robot control method provided by the embodiment of the present application can be applied to a terminal, can also be applied to a server, or can be software running on a terminal or a server. In some embodiments, the terminal can be a smart phone, a tablet computer, a laptop computer, a desktop computer, a smart speaker, a smart watch, a vehicle-mounted terminal, etc., but is not limited thereto; the server side can be configured as an independent physical server, can also be configured as a server cluster or a distributed system composed of multiple physical servers, and can also be configured as a cloud server providing basic cloud computing services such as cloud services, cloud databases, cloud computing, cloud functions, cloud storage, network services, cloud communications, middleware services, domain name services, security services, CDN, and big data and artificial intelligence platforms. The server can also be a node server in a blockchain network; the software can be an application that implements a quadruped robot control method, etc., but is not limited to the above forms.

[0119] The present application can be used in many general or special computer system environments or configurations. For example: personal computers, server computers, handheld or portable devices, tablet devices, multiprocessor systems, microprocessor-based systems, set-top boxes, programmable consumer electronic devices, network PCs, minicomputers, mainframe computers, distributed computing environments including any of the above systems or devices, and so on. The present application can be described in the general context of computer-executable instructions executed by a computer, such as program modules. Generally, program modules include routines, programs, objects, components, data structures, etc. that perform specific tasks or implement specific abstract data types. The present application can also be practiced in a distributed computing environment where tasks are performed by remote processing devices connected through a communication network. In a distributed computing environment, program modules can be located in local and remote computer storage media including storage devices.

[0120] The following specifically explains the embodiments of the present application in conjunction with the drawings:

[0121] Figure 1 is an optional flowchart of a quadruped robot control method provided by an embodiment of the present application, Figure 1 The method in can include but is not limited to steps S101 to S106.

[0122] Step S101, obtain the joint angle information of the quadruped robot. The joint angle information of the quadruped robot includes the hip adduction / abduction angle, the hip flexion / extension angle, and the knee bend angle;

[0123] Step S102, calculate the rotation error vector according to the joint angle information of the quadruped robot;

[0124] Step S103: Construct a model predictive control objective function based on the rotational error vector. The optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state.

[0125] Step S104: Perform quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the foot end force.

[0126] Step S105: Generate joint torque commands using the virtual force model control method.

[0127] Step S106: Perform robot control based on the predicted value of the foot end force and the joint torque commands.

[0128] Steps S101 to S106 illustrated in the embodiments of the present application achieve the control of the quadruped robot, improve the environmental adaptability and anti-rollover ability, and reduce the computational complexity and response delay.

[0129] In step S101 of some embodiments, the joint angle information of the quadruped robot can be obtained through the robot information library. It can also be obtained by other means, which is not limited thereto. Among them, the joint angle information of the quadruped robot may include the hip adduction / abduction angle, hip flexion / extension angle, and knee bending angle.

[0130] In some embodiments, in step S102, calculating the rotational error vector based on the joint angle information of the quadruped robot may include, but is not limited to, the following steps:

[0131] Construct a joint angle vector based on the joint angle information of the quadruped robot;

[0132] Calculate the end-effector spatial coordinates corresponding to each leg according to the joint angle vector;

[0133] Calculate the end-effector position combination matrix according to multiple end-effector spatial coordinates;

[0134] Calculate the correction matrix according to the correction coefficient and the end-effector position combination matrix;

[0135] Construct the forward kinematic function corresponding to each leg according to the correction matrix;

[0136] Take the partial derivative of the forward kinematic function with respect to the joint angle to obtain the Jacobian matrix;

[0137] Set the yaw angle, pitch angle, and roll angle of the quadruped robot according to the Jacobian matrix;

[0138] Construct the basic rotation matrix about the yaw angle according to the yaw angle;

[0139] Construct the basic rotation matrix about the pitch angle according to the pitch angle;

[0140] Construct a basic rotation matrix about the roll angle according to the roll angle;

[0141] Combine the basic rotation matrix about the yaw angle, the basic rotation matrix about the pitch angle, and the basic rotation matrix about the roll angle to obtain the overall rotation matrix;

[0142] Construct an error matrix according to the overall rotation matrix and the desired rotation matrix;

[0143] Perform a logarithmic mapping process on the overall rotation matrix according to the error matrix;

[0144] Construct a rotation error vector according to the multiple overall rotation matrices after the logarithmic mapping process.

[0145] In some embodiments, the forward kinematics of the quadruped robot can be utilized to calculate the positions of the end effectors of the four legs of the quadruped robot by using the joint angles through forward kinematics. First, a joint angle vector can be constructed according to the joint angle information of the quadruped robot, where the expression of the joint angle vector is: q = [q1, q2, q3, q4, q5, q6, q7, q8, q9, q 10 , q 11 , q 12 T ​, where \(q1\) represents the hip adduction / abduction angle (adjusting the lateral position of the leg in the horizontal plane), \(q2\) represents the hip flexion / extension angle (controlling the forward and backward swing in the sagittal plane), \(q3\) represents the knee flexion angle (controlling the degree of leg extension), \(q(1:3)\) corresponds to the first leg, \(q(4:6)\) corresponds to the second leg, \(q(7:9)\) corresponds to the third leg, and \(q(10:12)\) corresponds to the fourth leg. Exemplarily, the robotic mechanical structure parameters have been determined in the design, and the specific parameters include 0.082, 0.175826, 0.192981, 0.00268261, 0.065, 0.085, with the unit of m. Among them, 0.085 is the length of the short link closest to the fuselage of the leg; 0.175826 can be regarded as the effective length of the "thigh" part of the quadruped robot, used to describe the main link from the hip joint to the knee joint; 0.192981 represents the length of the "calf" part, used to describe the link from the knee joint to the foot end; 0.065 is used for overall translation compensation when calculating the end position, making the calculation result more consistent with the actual installation position; 0.082 is the short link from the hip joint to the thigh connection. According to the joint angle vector, calculate the end spatial coordinates corresponding to each leg. Exemplarily, taking the first leg as an example, its coordinates in space include: in the x coordinate, \(x1 = 0.082\sin(q2)-0.175826\cos(q2)+0.192981\cos(q2 + q3)-0.00268261\sin(q2 + q3)+0.065\); in the y coordinate, \(y1 = 0.082\cos(q2)\sin(q1)-0.085\cos(q1)+0.175826\sin(q1)\sin(q2)-0.00268261\cos(q2 + q3)\sin(q1)-0.192981\sin(q1)\sin(q2 + q3)\); in the z coordinate, \(z1 = 0.00268261\cos(q1)\cos(q2 + q3)-0.082\cos(q1)\cos(q2)-0.175826\cos(q1)\sin(q2)-0.085\sin(q1)+0.192981\cos(q1)\sin(q2 + q3)\). For other legs, the formula is similar, but the angle index and some constant symbols can be adjusted accordingly.

[0146] Then, according to multiple end spatial coordinates, calculate the end position combination matrix. Among them, the end position combination matrix is: In the formula, \((x i , y i , z i ) are the coordinates of the \(i\)-th leg, where \(i = 1, 2, 3, 4\). According to the correction coefficient and the end position combination matrix, calculate the correction matrix. Among them, the correction matrix is: Based on the rsf, the rbf corrects the translation of each leg in the x and y coordinates to correct the deviation between the actual installation position of the foot end and the ideal model. The selection of the correction parameters of the correction matrix is based on the error relative to the completely ideal model in actual assembly or simulation model assembly. For example, after actual assembly, through laser measurement or calibration experiments, it is found that there is a fixed deviation between the theoretical foot end position and the actual position, or the error generated during simulation modeling. Adding correction parameters can make the calculation more in line with the actual assembly.

[0147] Then, according to the correction matrix, the forward kinematic function corresponding to each leg is constructed, and the partial derivative of the forward kinematic function with respect to the joint angle is obtained to get the Jacobian matrix. Exemplarily, the Jacobian matrix is used to describe the relationship between the end effector velocity and the joint velocity. After obtaining the forward kinematic expression, the partial derivatives of the relevant joint angles of each position equation are calculated, and all the partial derivatives are organized into a matrix, where each column corresponds to the partial derivative result of a joint angle, and each row corresponds to the components of the end effector in the x, y, and z directions. It can be to establish the forward kinematic expression of the end of each leg of the quadruped robot in the world coordinate system, denoted as The end position combination matrix rsf is the discrete calculation result of r(q), which directly reflects the end position calculated by the forward kinematic model at the current q. The Jacobian matrix J(q) is defined as the partial derivative of the forward kinematic function with respect to the joint angle: where each column corresponds to the partial derivative of a joint angle, and each row corresponds to the components of the end in the x, y, and z directions.

[0148] Then, according to the Jacobian matrix, the yaw angle, pitch angle, and roll angle of the quadruped robot are set. According to the yaw angle, the basic rotation matrix around the yaw angle is constructed. According to the pitch angle, the basic rotation matrix around the pitch angle is constructed. According to the roll angle, the basic rotation matrix around the roll angle is constructed. Exemplarily, the desired rotation matrix can be set as R d , the actual overall rotation matrix is R, the yaw angle of the quadruped robot is ψ, the pitch angle is θ, and the roll angle is φ. The basic rotation matrices around each coordinate axis are constructed. The basic rotation matrix around the z-axis (yaw angle) is: where ψ = RPY(3). The basic rotation matrix around the y-axis (pitch angle) is: where θ = RPY(2). The basic rotation matrix around the x-axis (roll angle) is: The basic rotation matrix around the yaw angle, the basic rotation matrix around the pitch angle, and the basic rotation matrix around the roll angle are combined to obtain the overall rotation matrix. Among them, combined in the ZYX order, the overall rotation matrix is: R = R z R y R x . According to the overall rotation matrix and the desired rotation matrix, an error matrix is constructed. Among them, the error matrix is: where R d is the desired rotation matrix.

[0149] Finally, based on the error matrix, a logarithmic mapping process is performed on the overall rotation matrix, and a rotation error vector is constructed based on multiple overall rotation matrices after the logarithmic mapping process. Exemplarily, to calculate the rotation error vector, a function for calculating the rotation error vector through the rotation matrix can be constructed. In attitude control, it is often necessary to compare the difference between the desired attitude and the actual attitude. The constructed function is used to obtain the rotation error vector ω. This vector describes the rotation required from the current attitude to the desired attitude, with its direction being the rotation axis and its length being the rotation angle. The obtained rotation error vector can be used as the input of the controller to generate control signals. For example, in a VMC controller or model predictive control, the rotation error vector is often used to calculate the angular velocity correction to help the robot reach the desired attitude faster. By performing a logarithmic mapping on the rotation matrix, the non-linear rotation problem is transformed into a linear error expression, thereby simplifying the design and solution of the control algorithm. If R is a rotation matrix, its logarithmic mapping is defined as: where is the skew-symmetric matrix of the rotation axis. The constructed rotation error vector is: where is the scaling factor, which adjusts the unscaled rotation vector to the true rotation vector, that is, the rotation axis u multiplied by the rotation angle θ. When θ is very small, to avoid numerical instability caused by dividing by sin(θ) close to zero, ω / 2 is directly used as an approximation. At this time, the rotation angle is very small, and it can be considered that sin(θ)≈θ, so the scaling factor is approximately 1 / 2. The function finally outputs the rotation error vector ω, whose direction is the rotation axis and the modulus is the rotation angle θ.

[0150] In some embodiments, in step S103, constructing the model predictive control objective function according to the rotation error vector may include but is not limited to the following steps:

[0151] Construct a translational motion model based on the total mass of the robot, the centroid acceleration, the foot end reaction force, and the gravitational acceleration;

[0152] Construct a rotational motion model based on the inertia matrix of the robot about the centroid, the rotation error vector, the angular acceleration, the vector from the robot centroid to the leg contact point, and the foot end reaction force;

[0153] Construct a state vector based on the centroid position, the centroid velocity, the rigid body attitude parameters, and the angular velocity;

[0154] Construct a translational part dynamic model based on the translational motion model and the state vector;

[0155] Construct a dynamic model for the rotational part based on the rotational motion model and the state vector;

[0156] Construct a state matrix according to the translational part dynamic model;

[0157] Construct a first input matrix according to the rotational part dynamic model;

[0158] Construct an overall continuous model based on the state matrix, the first input matrix, the state vector, and the control input vector;

[0159] Calculate the discretized state transition equation according to the discrete time step and the basic equation of Euler integration method;

[0160] Construct a substitution model equation based on the overall continuous model and the discretized state transition equation;

[0161] Construct a discrete state transition matrix and a control influence matrix according to the substitution model equation;

[0162] Construct a discrete-time model based on the discrete state transition matrix and the control influence matrix;

[0163] Perform an approximation on the translational part according to the discrete-time model to obtain an approximate translational equation set;

[0164] Perform an approximation on the rotational part according to the discrete-time model to obtain an approximate rotational equation set;

[0165] Construct a model predictive control objective function based on the approximate translational equation set, the approximate rotational equation set, the prediction model, and the desired state.

[0166] In some embodiments, a quadruped robot MPC (Model Predictive Control) controller can be designed. To achieve real-time gait control, in this embodiment, the quadruped robot is simplified to a single rigid body model. This model regards the robot as a mass point and a rigid body, ignoring the dynamic effects of each joint and the swinging legs. In this embodiment, the Newton-Euler equations are used to describe the translational and rotational motions of the robot respectively, and the end positions of each leg are calculated through forward kinematics and the Jacobian matrix to provide necessary geometric information for dynamic modeling. To reduce the model complexity and computational burden, it can be assumed that the whole robot is simplified to a rigid body, whose center of mass position p and attitude are represented by the rotation matrix R, and the independent dynamic effects of joints and each swinging leg are ignored. First, a translational motion model can be constructed according to the total mass of the robot, the acceleration of the center of mass, the ground reaction force at the foot end, and the gravitational acceleration. The translational motion of the center of mass can be described by Newton's law, and the translational motion model is: In the formula, m is the total mass of the robot, is the acceleration of the center of mass, f iThe end - effector reaction force applied to the $i$-th leg, and $\mathbf{g}$ is the gravitational acceleration vector. It can be understood that the end - point positions $\mathbf{p}$ of each leg can be calculated based on the joint angles $\mathbf{q}$ i After that, based on this, the direction and magnitude of the force applied by each leg are determined through the contact state of the supporting legs, and the Jacobian matrix is further used to map the joint velocity to the end - effector velocity. According to the inertia matrix of the robot about the center of mass, the rotation error vector, the angular acceleration, the vector from the center of mass of the robot to the leg contact point, and the end - effector reaction force, a rotational motion model is constructed. Exemplarily, the rotational motion model can be described by the Euler's equation as: In the formula, $I$ xx $=\iiint$ V $(y$ 2 $+z$ 2 )$\rho dv$, $I$ yy $=\iiint$ V $(x$ 2 $+z$ 2 )$\rho dv$, $I$ zz $=\iiint$ V $(y$ 2 $+x$ 2 )$\rho dv$, $I$ xy $=\iiint$ V $(xy)\rho dv$, $I$ xz $=\iiint$ V $(xz)\rho dv$, $I$ yz $=\iiint$ V $(yz)\rho dv$, $I$ is the inertia matrix of the robot about the center of mass, $\omega$ is the angular velocity of the rigid body of the robot, which is the actual angular motion quantity of the robot and appears in the Euler's equation of rigid - body rotation, and is used to describe and predict the rotational dynamics of the rigid body under the action of external torque. The error angular velocity represents the deviation between the current attitude and the desired attitude, and is used to generate a compensation signal in PD control to make the robot adjust its attitude to track the desired target. $\dot{\omega}$ is the derivative of the angular velocity of the robot's rigid body with respect to time, that is, the angular acceleration, $\mathbf{r}$ i is the vector from the center of mass of the robot to the contact point of the $i$-th leg, $\mathbf{f}$ i is the end - effector reaction force applied by this leg. Moreover, the outer contour of the body of the quadruped robot is approximately a symmetric and regular cuboid, so the non - diagonal elements in the inertia tensor here can be neglected. The simplified inertia matrix is: The end - effector positions of each leg are obtained through forward kinematics, and then each $\mathbf{r}$ i is obtained by using the rigid - body model, and then the rotational equation is constructed. For the attitude error, the logarithm of the rotation matrix can be mapped to a rotation vector, so as to obtain the difference between the current attitude and the desired attitude.

[0167] Then, a state vector is constructed based on the centroid position, centroid velocity, rigid body attitude parameters, and angular velocity. The state vector is as follows: In the formula, represents the centroid position, represents the centroid velocity, represents the rigid body attitude parameters (such as the rotation vector obtained through logarithmic mapping), represents the angular velocity. Based on the translational motion model and the state vector, a translational part dynamic model is constructed. The translational part dynamic model can be described according to Newton's second law as: Where, f i (t) represents the forces applied at the end of each leg, m is the total mass of the robot, and g is the gravitational acceleration vector. Based on the rotational motion model and the state vector, a rotational part dynamic model is constructed. The rotational part dynamic model can be described based on the rigid body Euler equation (after linearization) as: Where, represents the angular velocity, r i is the vector from the centroid to the contact point of the i-th leg, and I is the inertia matrix of the robot base in the body coordinate system.

[0168] Then, according to the translational part dynamic model, a state matrix is constructed. According to the rotational part dynamic model, a first input matrix is constructed. And according to the state matrix, the first input matrix, the state vector, and the control input vector, an overall continuous model is constructed. The overall continuous model described in state space form is: In the formula, A c is the state matrix, B c is the first input matrix, and u(t) is the control input vector. According to the discrete time step and the basic equation of the Euler integration method, the discretized state transition equation is calculated. The basic equation of the Euler integration method is: The discretized state transition equation is: In the formula, x k = x(kΔt), u k = u(kΔt), and Δt is the discrete time step. According to the overall continuous model and the discretized state transition equation, a substitution model equation is constructed. The discretized state transition equation can be substituted into the overall continuous model to obtain the substitution model equation: x k+1 = x k + Δt[A c x k + B c u k . According to the substitution model equation, a discrete state transition matrix and a control influence matrix are constructed. The discrete state transition matrix is: A = I + Ac For Δt, the control influence matrix is: B = B c Δt. Based on the discrete state transition matrix and the control influence matrix, a discrete-time model is constructed, where the discrete-time model is: x k+1 = Ax k + Bu k , where A is the discrete state transition matrix and B is the control influence matrix.

[0169] Finally, based on the discrete-time model, an approximation process is performed on the translational part to obtain a translational approximation equation set, where the translational approximation equation set is: Based on the discrete-time model, an approximation process is performed on the rotational part to obtain a rotational approximation equation set, where the rotational approximation equation set is: Based on the translational approximation equation set, the rotational approximation equation set, the prediction model, and the desired state, a model predictive control objective function is constructed, where the optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state. Exemplarily, in the objective function optimized by MPC, the error between the prediction model X = Sx0 + TU and the desired state X ref is minimized, and the model predictive control objective function is: It can be understood that the translational and rotational approximation equations ensure that X contains information on both the translational motion of the center of mass and the rotational motion of the attitude, so that the MPC optimization not only tracks the position and velocity of the center of mass, but also considers the change in the attitude (rotation) state. In addition, the inequality constraints are also based on these state predictions to ensure that the predicted trajectory satisfies physical constraints (such as friction, force limits, etc.).

[0170] In some embodiments, in step S104, according to the model predictive control objective function, a quadratic programming solution process is performed to obtain the predicted value of the foot force, which may include but is not limited to steps S201 to S205:

[0171] Step S201: Construct a compression equation based on the foot reaction force;

[0172] Step S202: Construct a standard quadratic programming objective function based on the compression equation and the model predictive control objective function;

[0173] Step S203: Construct foot constraints, where the foot constraints include foot force limit constraints and friction constraints;

[0174] Step S204: Combine the standard quadratic programming objective function and the foot constraints to obtain a standard quadratic programming problem;

[0175] Step S205: Use a quadratic programming solver to solve the standard quadratic programming problem to obtain the predicted value of the foot force.

[0176] In some embodiments, in the design process of the quadruped robot mpc controller, a compression equation can be constructed based on the foot-end reaction force first, and then a standard quadratic programming objective function can be constructed according to the compression equation and the model predictive control objective function. Next, foot-end constraints are constructed, including foot-end force limit constraints and friction constraints. Finally, the standard quadratic programming objective function and the foot-end constraints are combined to obtain a standard quadratic programming problem, and a quadratic programming solver is used to solve the standard quadratic programming problem to obtain the predicted value of the foot-end force. Exemplarily, in the QP optimization process, the standard quadratic programming objective function and the foot-end constraints can be combined to obtain the standard quadratic programming problem as follows: In the construction process of each matrix and vector, a prediction table can be constructed first. According to the discrete dynamics model, matrices Q and T are constructed such that X = Sx0 + TU. Then, the objective function is constructed, and H and f are constructed using the state error X - X ref and the control input penalty: H = 2(T T S d T + W d ), f = 2T T S d (Sx0 - X ref ). Then, the equality constraints (if explicitly included) are constructed, and A control and b control are generated according to the state transition equation. The inequality constraints are constructed, and combined with the friction pyramid constraints, force limits, etc., and they are arranged into the form of A ineq U ≤ b ineq . Finally, the QP solver is called, and H, f, A control , b control , A ineq , b ineq are all passed to the qpOASES solver to obtain the optimal predicted value of the foot-end force U * .

[0177] In some embodiments, in step S201, constructing the compression equation according to the foot-end reaction force may include but is not limited to the following steps:

[0178] Set the control quantity according to the foot-end reaction force;

[0179] Calculate the compensation acceleration according to the centroid position, centroid velocity, desired position, and desired velocity;

[0180] Calculate the angular acceleration according to the gain matrix of the attitude error, the gain matrix of the angular velocity error, the attitude error vector, the desired angular velocity, and the current angular velocity;

[0181] Take the compensation acceleration and the angular acceleration as the state tracking targets;

[0182] Construct the friction cone constraint according to the foot-end force and the ground friction coefficient;

[0183] Perform a linear transformation on the frictional cone constraint to obtain a system of linear inequalities, which is used as the energy consumption objective.

[0184] According to the control variables, combine the state tracking objective and the frictional cone constraint to obtain a first programming problem. The parameters of the first programming problem include the foot-end force, the weighted matrix of the state error, and the weighted matrix of the control input.

[0185] Based on the first programming problem, the continuous-time system matrix, and the second input matrix, construct the relationship equation between the robot state and the control input.

[0186] According to the second time step, discretize the relationship equation between the robot state and the control input to obtain the discretized model.

[0187] Perform a compression transformation on the discretized model to obtain the compression equation.

[0188] In some embodiments, the control variables can be set first according to the foot-end reaction force. Exemplarily, the control variable u k is selected as the foot-end reaction forces (or equivalent components) exerted by the four legs at the discrete time k, that is, where, represents the force vector exerted by the i-th leg at time k. A single u k usually has a dimension of 12, corresponding to the (f x , f y , f z ) components of the four legs. The main control objective of this embodiment is to make the actual state (position, velocity, attitude) of the robot track the desired state as much as possible, and at the same time minimize the control energy consumption on the premise of ensuring physical feasibility. For this purpose, two types of objectives can be defined, including the state tracking objective and the energy consumption objective. The compensation acceleration can be calculated according to the centroid position, centroid velocity, desired position, and desired velocity. Exemplarily, assume that the desired state at the discrete time k is x ref,k , and the actual state is obtained from the discrete dynamic equation as x k . Then it is desired to minimize the state error ||x k - x ref,k ||. To enable the robot's centroid to smoothly track the expected trajectory, a state tracking control law based on PD feedback can be designed. Let r represent the actual centroid position of the robot, v represent the actual centroid velocity of the robot, and r ref represent the desired position, v ref represent the desired velocity, then the compensation acceleration a is defined as: a = K pcom (r ref - r) + K dcom (v ref - v), where the proportional gain Kpoom For reducing the position error, differential gain K dcom For attenuating the velocity difference, overall forming a virtual spring-damping system. This compensated acceleration can not only correct the centroid motion error in real time, but also generate a reference force by combining with the gravity term (i.e., F = m(a + g)), and then be input as the target dynamic quantity in the MPC optimization problem to achieve the coordinated control of the overall dynamic behavior of the robot.

[0189] Then, according to the gain matrix of the attitude error, the gain matrix of the angular velocity error, the attitude error vector, the desired angular velocity, and the current angular velocity, calculate the angular acceleration, and use the compensated acceleration and the angular acceleration as the state tracking targets. Exemplarily, to achieve precise control of the robot's attitude (roll, pitch, yaw), a similar PD control strategy can be adopted in the attitude space. Through the logarithmic mapping of the rotation matrix, extract the attitude error vector qw, and combine it with the angular velocity error, and define the angular acceleration aw as: aw = K pbase (qw) + K d_base (w ref - w), where qw represents the rotation vector error ω between the current attitude and the desired attitude. Here, qw is described to distinguish the actual angular velocity, w ref is the desired angular velocity (set to 0, indicating that it is desired to maintain a stable attitude), w is the actual angular velocity, K p＿base is the attitude error, K d＿base is the gain matrix of the angular velocity error. When aw is mapped through the inertia matrix I, i.e., τ = Iaw, the attitude control torque τ can be obtained, which is used to compensate and correct the deviation of the robot base in attitude. Together with the force F obtained from the centroid translation part, these forces and torques are finally synthesized into the overall control quantity b control = [F; τ], and play a role in the subsequent optimization and control process.

[0190] In the calculation of the energy consumption target, a friction cone constraint can be constructed according to the foot-end force and the ground friction coefficient, and the friction cone constraint is linearly transformed to obtain a system of linear inequalities, where the system of linear inequalities is used as the energy consumption target. Exemplarily, to avoid excessive or ineffective consumption of the foot-end force, it is necessary to impose a penalty on the control input u k , that is, try to make the force distribution meet the motion requirements while ensuring that the magnitude of the force is appropriate. Usually, a quadratic penalty term of ‖u k ‖ is added to the cost function. Combining these two types of targets forms a quadratic cost function, which will be specifically introduced in the subsequent optimization problem section. To prevent the foot-end from sliding during the support process, a friction constraint of the foot-end force needs to be imposed. Let the foot-end force f i,k = [f x,i , f y,i , f z,i ​T (with the z - axis vertically upward), the ground friction coefficient is μ. In an ideal situation, the friction cone constraint can be expressed as: Converting the friction cone constraint into a number of linear inequalities gives: where α = tan(θ) is the tangent value of the cone surface opening angle.

[0191] Then, according to the control quantity, the state tracking target and the friction cone constraint are combined to obtain the first planning problem. Among them, the parameters of the first planning problem include the foot - end force, the weighted matrix of the state error, and the weighted matrix of the control input. Exemplarily, the state tracking target and the friction cone constraint can be combined to obtain a quadratic programming problem with equality and inequality constraints, that is, the first planning problem. Let U denote the vector formed by splicing the control inputs at all time steps within the prediction window. For a quadruped robot, the control input at each moment is usually the foot - end force u k of each leg. Then, for the prediction window of N steps, Let S be the weighted matrix of the state error, which is used to penalize the deviation between the predicted state and the desired state. Its dimension is the same as that of the state vector x. Let W be the weighted matrix of the control input, which is used to penalize the magnitude of the control input. Its dimension is the same as that of a single control input u k is the same. Let x ref (also denoted as x ref,k ) represent the desired state at each time step within the prediction window. The state usually includes information such as the position of the robot's center of mass, velocity, and attitude. It is obtained through gait planning or the upper - layer controller and serves as the state tracking target. Since at each discrete moment, based on PD control, a compensation acceleration a and a compensation angular acceleration aw are generated and added to the reference state advancement at the next moment, so that x ref,k+1 is adaptively adjusted according to the current error. Therefore, when constructing ||x k+1 - x ref，k+1 ||, the effect of PD compensation is already implicit.

[0192] To predict the system state online, it is necessary to first discretize the continuous - dynamics model. According to the first planning problem, the continuous - time system matrix, and the second input matrix, a relationship equation between the robot state and the control input can be constructed. Among them, the relationship equation between the robot state and the control input is: where, is the system state (including the position of the center of mass, velocity, attitude, etc.), is the control input (the foot-end force of the four legs or its equivalent components). To predict the future state at each discrete moment, it is necessary to discretize the above continuous system to a time step Δt. The relationship equation between the robot state and the control input can be discretized according to the second time step to obtain a discretized model; let k represent the discrete time index, then the discretized model is: x k+1 = Ax k + Bu k , where A = exp(A c Δt), Perform a compression transformation on the discretized model to obtain a compression equation. After recursion, all future states can be written in a compressed form, and the obtained compression equation is: X = Qx0 + TU, where X is the concatenated vector of states within the prediction window, Q is the matrix representing the influence of the initial state on the future state, and T is the matrix representing the influence of the control input on the future state.

[0193] In some embodiments, in step S202, according to the compression equation and the model predictive control objective function, constructing a standard quadratic programming objective function may include, but is not limited to, the following steps:

[0194] Rewrite the model predictive control objective function according to the quadratic penalty of the state error term and the control input term to obtain a first objective function;

[0195] Rewrite the first objective function according to the compression equation and the desired state to obtain a second objective function;

[0196] Expand the second objective function to obtain a standard quadratic form objective function;

[0197] Calculate the overall control quantity according to the inertia matrix of the robot base in the body coordinate system and the overall rotation matrix;

[0198] Construct a mapping matrix according to the sum of the foot-end forces of each leg, the rotational effect generated by the foot-end forces of each leg, and the skew-symmetric matrix function;

[0199] Rewrite the standard quadratic form objective function according to the overall control quantity and the mapping matrix to obtain a third objective function;

[0200] Expand the third objective function to obtain a fourth objective function;

[0201] Perform a form transformation on the fourth objective function to obtain a standard quadratic programming objective function.

[0202] In some embodiments, first, rewrite the model predictive control objective function according to the quadratic penalty of the state error term and the control input term to obtain a first objective function, where the first objective function is: According to the compression equation and the desired state, rewrite the first objective function to obtain the second objective function. The compressed expression X = Qx0 + TU can be concatenated with the desired state X ref , and the second objective function is obtained as: J(U) = (Qx0 + TU - X ref ) T W d (Qx0 + TU - X ref ) + U + S d U, where S d and W d are both block diagonal matrices, and S and W act on the state and control input at each moment respectively. Expand the second objective function to obtain the standard quadratic form objective function as: where H = 2(T T S d T + W d ), f = 2T T S d (Qx0 - X ref ), and C is a constant term. According to the inertia matrix of the robot base in the body coordinate system and the overall rotation matrix, calculate the overall control quantity. Exemplarily, when actually constructing the objective function, in addition to the state tracking term, the end-effector force u and the overall control quantity b control are related through the mapping matrix A control . The overall control quantity is: where τ is the torque required by the base, τ = I world aw, I world = RIR T , I is the inertia matrix of the robot base in the body coordinate system, and R is the overall rotation matrix.

[0203] Then, according to the sum of the end-effector forces of each leg, the rotational effect generated by the end-effector forces of each leg, and the skew-symmetric matrix function, construct the mapping matrix. Exemplarily, when defining the mapping matrix A control , the end-effector forces U of each supporting leg can be mapped to the overall translational force and rotational torque. This matrix is divided into two parts: The first part directly takes the sum of the end-effector forces of each leg and multiplies it by the flag factor flag, where flag(i) is 1 or 0, and is used to indicate whether the i-th leg is in the support phase or is activated in the current QP optimization process. When flag(i) = 1, it means that the leg is in the support state and its end-effector force should participate in the optimization; while when flag(i) = 0, it means that the leg is in the swing phase and its end-effector force does not participate in the QP optimization of the current support phase. The second part uses the Skew matrix to represent the rotational effect generated by the end-effector forces of each leg. Formally, the matrix A control is defined as: where rbf(:,i) is the end position of the i-th leg after forward kinematics calculation and correction, and Skew(·) represents taking the skew-symmetric matrix of this matrix. Using A control , theoretically the overall output is: A control U, and the foot end force U and b control satisfy the relationship: A control U≈b control , where A control contains both the translational part (direct action) and the rotational part (describing the rotational effect through the Skew matrix of the foot end position).

[0204] Then, according to the overall control quantity and the mapping matrix, the standard quadratic objective function is rewritten to obtain the third objective function. Among them, the third objective function is: Among them, represents the weighted two-norm, and the weight matrix S is selected as a diagonal matrix, which is diag([1,1,10,50,30,10]), used to balance the relative importance of the translational and rotational parts. For a quadruped robot, the stability of some components (such as roll, pitch, yaw) is crucial because they directly determine whether it will tip over or become unbalanced. In contrast, the x and y horizontal translational errors can be tolerated moderately, so a lower weight can be given in the weighted matrix. Therefore, the weight order of the diagonal matrix is x, y, roll, pitch, yaw. Since in the current QP construction, the optimization problem is established for the matching between the foot end force and the desired overall control quantity, and the influence of the initial state usually appears as a constant term and does not affect the solution of the optimization variables, so Qx0 can be ignored. represents the regularization term for U, and the weight matrix W is 0.001I 12x12 , which is a 12×12 identity matrix, indicating that the same regularization is applied to the foot end force components of the four legs (3 components for each leg, a total of 12 components). The advantage of this design is that the control inputs in all directions are treated equally, and the penalty size is only adjusted by the scalar 0.001, making the optimization problem numerically well-conditioned. α is the regularization factor (which can take the value 0.01).

[0205] Finally, the third objective function is expanded to obtain the fourth objective function, and the fourth objective function is transformed in form to obtain the standard quadratic programming objective function. Exemplarily, the expanded fourth objective function is: J(U)=(A control U - b control ) T S(A control U - b control ) + αU T WU, which can also be described as: In the formula, The standard quadratic programming objective function obtained after the formal transformation is as follows: In this embodiment, a compensation acceleration generated by PD control is added to the centroid translational part to correct the desired trajectory or the desired force. At the same time, to achieve precise control of the robot's attitude, a similar PD control strategy is adopted in the attitude space to generate a compensation angular acceleration. When constructing the MPC optimization problem, let the reference state x ref,k+1 or the reference force F ref,k include a compensation term, so that the objective function of the QP can adaptively respond to the current error. The online optimization implemented thereby can not only track static or pre-planned trajectories, but also dynamically compensate for the centroid position and velocity deviation and the attitude deviation caused by model uncertainty or external disturbances.

[0206] In some embodiments, in step S203, constructing the foot end constraints may include, but are not limited to, the following steps:

[0207] Construct a foot end force limit constraint according to the upper limit and the lower limit of the foot end force;

[0208] Construct a friction constraint matrix by using the friction pyramid approximation according to the basic matrix of friction;

[0209] Construct a friction constraint according to the friction constraint matrix, the upper bound of the vertical component, and the lower bound of the vertical component.

[0210] In some embodiments, foot end constraints may be constructed, where the foot end constraints include a foot end force limit constraint and a friction constraint. The inequality constraints mainly include friction constraints, friction constraints, and other kinematic or dynamic limitations. In the friction constraint, for each leg at each moment, using the friction pyramid approximation, the constraint can be written as: -αf z,i,k ≤f x,i,k ≤αf z,i,k , -αf z,i,k ≤f y,i,k ≤αf z,i,k , f z,i,k ≥0. In the friction constraint, there are maximum and minimum limits for the foot end force, written as: u min ≤u k ≤u max . In other kinematic or dynamic limitations, such as the centroid position or attitude range, leg length constraints, etc., they can also be linearized and written in the form of inequalities. After organizing the overall inequality constraints, they can be expressed as: A ineq U≤b ineq , where each row of A ineq corresponds to a linear constraint, and b ineq is the corresponding upper bound. In this embodiment, the foot end force limit constraint can be constructed first according to the upper limit and the lower limit of the foot end force. Exemplarily, for the foot end force component f of each legi = [f i,x , f i,y , f i,z T , the upper limit and lower limit of the foot-end force can be set, and the foot-end force boundary constraint is: -flag(i)F max ≤ f i,j ≤ flag(i)F max , j = 1, 2, 3, where F max is the upper limit of the foot-end force, which is a sufficiently large value (e.g., 100000), and flag(i) ∈ {0, 1} is an indicator variable indicating whether the leg is in the supporting state. In this way, on the supporting leg, the foot-end force can have a sufficiently large value range, while for the leg not participating in the support (flag(i) = 0), its contribution is masked.

[0211] Then, according to the fundamental matrix of friction, a friction constraint matrix is constructed using the friction pyramid approximation, and based on the friction constraint matrix, the upper bound of the vertical component, and the lower bound of the vertical component, a friction constraint is constructed. Exemplarily, to prevent the foot-end from slipping, it is necessary to satisfy the friction cone constraint (i.e., the friction constraint): However, since the friction cone constraint is non-linear, the friction pyramid approximation method can be used to transform it into a set of linear inequalities. The fundamental matrix is defined as: where uf = 0.5 is the approximation parameter. For each leg, after multiplying by the indicator variable flag(i), the overall friction constraint matrix is constructed as: And the upper bound vector lbA and the lower bound vector ubA are set respectively. For example, for the first row constraint of each leg, lbA = -flag(i) × 100000 and ubA = 0 are set; for the vertical component, the upper bound is flag(i) × f max (such as 160) and the lower bound is flag(i) × 10, and these values are determined according to the actual hardware and friction characteristics. After integrating the friction cone constraint and the foot-end force boundary constraint, it can be uniformly written in the standard form: A ineq f ≤ b ineq , where This expression ensures that when solving the QP problem, the optimization variable f must simultaneously satisfy the upper and lower bounds of the foot-end force and the friction constraint, thus ensuring the stability and anti-slip property of the robot during the support phase.

[0212] In some embodiments, in step S105, using the virtual force model control method to generate the joint torque command may include, but is not limited to, the following steps:

[0213] In horizontal trajectory planning, interpolation is performed according to the initial position, the target position, and the cubic Bezier curve interpolation equation to obtain the horizontal trajectory;

[0214] ​Derive the horizontal trajectory to obtain the horizontal velocity and horizontal acceleration;

[0215] In vertical trajectory planning, calculate the velocity and acceleration in the rising phase according to the initial coordinates, target coordinates, lifting height, local time variable in the rising phase, and cubic Bezier curve interpolation equation;

[0216] According to the initial coordinates, target coordinates, lifting height, local time variable in the descending phase, and cubic Bezier curve interpolation equation, calculate the velocity and acceleration in the descending phase;

[0217] Calculate the expected value of the swing trajectory according to the horizontal velocity, horizontal acceleration, velocity in the rising phase, acceleration in the rising phase, velocity in the descending phase, and acceleration in the descending phase. The expected value of the swing trajectory includes the expected foot end position and the expected foot end velocity;

[0218] Calculate the corresponding foot end position of each leg in the robot coordinate system according to the position of the center of mass, the expected foot end position, and the fixed offset between the foot end and the geometric reference of the robot body;

[0219] Calculate the corresponding foot end velocity of each leg in the robot coordinate system according to the center of mass velocity and the expected foot end velocity;

[0220] Calculate the current foot end velocity according to the Jacobian matrix and the joint angular velocity;

[0221] Calculate the virtual control force according to the corresponding foot end position of each leg in the robot coordinate system, the corresponding foot end velocity of each leg in the robot coordinate system, the current foot end velocity, and the current foot end position;

[0222] Generate the joint torque command according to the virtual control force.

[0223] In some embodiments, in the main control process, the system determines which legs are in the swing phase. The VMC controller calculates the swing progress of each leg according to information such as the global phase, offset, duration, etc. If a leg is in the swing phase, a swing phase controller is required for leg motion planning and virtual model control. Virtual model control (VMC) realizes compliant control by introducing a "virtual spring-damper" system between the foot end and the expected position. Its essence is a PD control: F vmc =K p (r des -r)+K d (v des -v), where K p corresponds to the spring stiffness (proportional gain), K d corresponds to the damping coefficient (differential gain), r des is the expected position, v desis the desired speed, r is the actual foot-end position, and v is the actual foot-end speed. F vmc represents the virtual foot-end force calculated in the virtual model control (VMC). The basic idea is to construct a virtual spring-damper system to generate a desired feedback force for the swinging leg or the supporting leg during trajectory tracking. This virtual force is converted into joint torques through mapping (such as the joint Jacobian matrix) to drive the motors to achieve smooth swinging of the foot-end. In the swinging trajectory planning, to achieve smooth and continuous movement of the foot-end during the swinging phase of the quadruped robot, a trajectory planning method based on Bezier curves is adopted. This method not only ensures smooth changes in the horizontal plane (x, y) trajectory but also designs a segmented planning scheme for the vertical direction (z) leg-lifting requirement, including horizontal trajectory planning and vertical trajectory planning.

[0224] In the horizontal trajectory planning, according to the initial position, target position, and the cubic Bezier curve interpolation equation, interpolation is performed to obtain the horizontal trajectory, and the horizontal trajectory is differentiated to obtain the horizontal speed and horizontal acceleration. Exemplarily, for the movement of the foot-end in the horizontal direction, we assume its initial position is the target position is During the entire swinging period, let the normalized time parameter t ∈ [0, 1] (i.e., phase), and cubic Bezier curve interpolation is used. The interpolation formula is: B(t) = t 3 + 3t 2 (1 - t), and the actually generated horizontal trajectory is: p xy (t) = p xy,init +(t 3 + 3t 2 (1 - t))(p xy,final - p xy,init ), since t 3 + 3t 2 (1 - t) = 3t 2 - 2t 3 , differentiating the horizontal trajectory curve, the speed is obtained as: the acceleration is obtained as: To map the speed and acceleration to the real-time scale, the speed can also be divided by the swinging time T, and the acceleration is divided by T 2 .

[0225] In the vertical trajectory planning, since the foot-end needs to lift and then descend during the swinging phase, a segmented strategy is adopted for the vertical trajectory planning. Let the initial z coordinate be z init , and the target z coordinate be z final, and the height of lifting the leg is h. Then the planning is divided into an ascending stage and a descending stage. In the ascending stage (0 ≤ t < 0.5), the velocity and acceleration in the ascending stage can be calculated according to the initial coordinates, target coordinates, leg-lifting height, local time variable in the ascending stage, and the cubic Bezier curve interpolation equation. Exemplarily, let the local time variable be τ = 2t, τ ∈ [0, 1]. In the ascending stage, the foot end rises from z init to z init + h. Using the same interpolation formula as the horizontal planning: p z (t) = z init +(3τ 2 -2τ 3 )h. The corresponding velocity in the ascending stage is: The acceleration in the ascending stage is: where T is the swing time. In the descending stage (0.5 ≤ t ≤ 1), the velocity and acceleration in the descending stage can be calculated according to the initial coordinates, target coordinates, leg-lifting height, local time variable in the descending stage, and the cubic Bezier curve interpolation equation. Exemplarily, let the local time variable be τ = 2t - 1, τ ∈ [0, 1]. In the descending stage, the foot end descends from z init + h to z final . The planning formula is: p z (t)=(z init + h)+(3τ 2 -2τ 3 )(z final -(z init + h)). The corresponding velocity in the descending stage is: The acceleration in the descending stage is:

[0226] Then, according to the horizontal velocity, horizontal acceleration, velocity in the ascending stage, acceleration in the ascending stage, velocity in the descending stage, and acceleration in the descending stage, the expected value of the swing trajectory is calculated. The expected value of the swing trajectory includes the expected foot end position and the expected foot end velocity. Exemplarily, by judging the value of t, two formulas in the ascending stage and the descending stage are respectively used for the z component, and the calculation results are assigned to the third component of the output vector. The horizontal and vertical components are spliced to obtain the complete swing leg trajectory, whose position, velocity, and acceleration are all continuous, meeting the smoothness requirements of the robot gait planning. Finally, the output expected foot end position is: The expected foot end velocity is: The expected foot end acceleration is:

[0227] To describe the trajectory of the swinging leg in the robot body coordinate system, the desired position and velocity in the world coordinate system can be transformed into the robot coordinate system. The corresponding foot position of each leg in the robot coordinate system can be calculated based on the centroid position, the desired foot position, and the fixed offset between the foot and the geometric reference of the robot body. And the corresponding foot velocity of each leg in the robot coordinate system can be calculated based on the centroid velocity and the desired foot velocity. Exemplarily, denote the current centroid position of the robot as and the rotation matrix R associates the robot body coordinate system with the world coordinate system. Then the foot position of the i-th leg in the robot coordinate system is: In the formula, is the foot position of the i-th leg in the robot coordinate system, represents the preset fixed offset between the foot and the geometric reference of the robot body. Similarly, the corresponding foot velocity of each leg in the robot coordinate system is: where is the desired centroid velocity, and v body is the centroid velocity of the robot.

[0228] Finally, based on the Jacobian matrix and the joint angular velocity, the current foot velocity is calculated. Exemplarily, the current foot position is obtained by the forward kinematics module, while the current foot velocity is calculated by the product of the Jacobian matrix and the joint angular velocity. The calculation formula for the current foot velocity is: In the formula, J (i) (q) is the Jacobian matrix of the i-th leg, is the joint angular velocity, which can also be the relevant velocity variable. Based on the corresponding foot position of each leg in the robot coordinate system, the corresponding foot velocity of each leg in the robot coordinate system, the current foot velocity, and the current foot position, the virtual control force is calculated. Exemplarily, in the swing phase, the virtual spring-damping effect is achieved through PD control. For the i-th swinging leg, its virtual control force has the following calculation formula: In the formula, K p,cart is the gain of the foot position (spring stiffness), which acts as a spring and generates a restoring force proportional to the position deviation. K d,cart is the gain of the velocity feedback (damping coefficient), which acts as a damper and suppresses the oscillation caused by the velocity error. Based on the virtual control force, the joint torque command is generated. Exemplarily, the virtual control force can be mapped through the joint Jacobian matrix to generate the joint space control command. The specific mapping formula is: In the formula, τ (i) is the joint torque command of the i-th root in the swing phase.

[0229] In some embodiments, in step S106, robot control can be performed based on the predicted foot-end force value and the joint torque command. A control signal can be generated according to the predicted foot-end force value, and the corresponding joint motors can be directly driven in combination with the joint torque command to achieve compliant tracking of the foot-end, thereby realizing the control of the quadruped robot, improving the environmental adaptability and anti-rollover ability, and reducing the computational complexity and response delay.

[0230] In some embodiments, anti-interference can be analyzed, including external force disturbances and complex terrain adaptation. In external force disturbances, when the robot is subjected to external impacts or wind disturbances, the compensation acceleration generated by PD feedback will quickly reflect the deviation between the center of mass and the desired state. While considering this deviation, MPC control redistributes the foot-end forces of the supporting legs through online QP solving to ensure that the robot can stably recover to the equilibrium state. At the same time, the VMC module of the swinging leg can also quickly adjust the foot-end trajectory through virtual forces to reduce the vibration caused by external disturbances. In complex terrain adaptation, when encountering rough or irregular terrain, the smooth trajectory generated by the foot-end trajectory planning module provides a compliant motion reference for the robot; while the PD control in the VMC module realizes compliant tracking locally, effectively buffering the instantaneous impact caused by terrain changes. At the same time, the friction and force limit constraints in MPC optimization ensure that the foot-end forces of the supporting legs are always within a reasonable range, preventing slipping and instability, thereby achieving robust adaptation to complex terrain.

[0231] In some embodiments, the anti-rollover ability in complex terrains is greatly improved in this embodiment. When predicting future states, the MPC in the stance phase of this embodiment effectively shortens the control delay through real-time optimization, ensuring a stable response under sudden disturbances. The VMC in the swing phase has the ability to quickly correct the foot-end. Even in complex environments, the error is controlled within a reasonable range, reducing abnormal fluctuations caused by sensor lag or intuitive control errors, and taking into account both the response speed and robustness. The hybrid control strategy adopted in this embodiment achieves a balance between low-energy consumption operation and high motion stability. The energy consumption ratio under complex terrains is significantly lower than that of traditional methods (such as within an energy consumption ratio of 1.3:1), optimizing the energy consumption; at the same time, the response delay is controlled below 200 ms, breaking through the traditional control cycle limit. The MPC in the stance phase of this embodiment uses an accurate dynamic model for state prediction and ensures stability under sudden external disturbances through anti-interference constraints; when the environmental complexity is high, the VMC in the swing phase ensures gait continuity through a hybrid adjustment of virtual forces and positions, reducing the decision error rate caused by model errors, and enhancing the system robustness and adaptive ability.

[0232] The beneficial effects of implementing the embodiments of the present invention include: The embodiments of the present invention first obtain the joint angle information of the quadruped robot, then calculate the rotation error vector according to the joint angle information of the quadruped robot, construct a model predictive control objective function according to the rotation error vector, perform quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the foot end force, generate a joint torque command using the virtual force model control method, and finally perform robot control according to the predicted value of the foot end force and the joint torque command, so as to be able to combine model predictive control and virtual force model control to achieve quadruped robot control, thereby improving the environmental adaptability and anti-rollover ability, and reducing the computational complexity and response delay.

[0233] As Figure 2 shown, the embodiments of the present invention also provide a quadruped robot control device, including:

[0234] The first module 801 is used to obtain the joint angle information of the quadruped robot, and the joint angle information of the quadruped robot includes the hip adduction / abduction angle, the hip flexion / extension angle, and the knee flexion angle;

[0235] The second module 802 is used to calculate the rotation error vector according to the joint angle information of the quadruped robot;

[0236] The third module 803 is used to construct a model predictive control objective function according to the rotation error vector, and the optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state;

[0237] The fourth module 804 is used to perform quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the foot end force;

[0238] The fifth module 805 is used to generate a joint torque command using the virtual force model control method;

[0239] The sixth module 806 is used to perform robot control according to the predicted value of the foot end force and the joint torque command.

[0240] The content in the above method embodiments is applicable to the device embodiments of the present invention. The functions specifically implemented by the device embodiments of the present invention are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those of the above method embodiments.

[0241] As Figure 3 shown, the embodiments of the present invention also provide a computer device, including:

[0242] At least one processor 901;

[0243] At least one memory 902, used to store at least one program;

[0244] When at least one program is executed by at least one processor such that the at least one processor implements Figure 1 the method shown.

[0245] The content in the above method embodiments is applicable to the device embodiments of the present invention. The functions specifically implemented by the device embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those of the above method embodiments.

[0246] The embodiments of the present invention further provide a computer-readable storage medium. The computer-readable storage medium stores a computer program, and when the computer program is executed by a processor, it implements Figure 1 the method shown.

[0247] The content in the above method embodiments is applicable to the storage medium embodiments of the present invention. The functions specifically implemented by the storage medium embodiments are the same as those of the above method embodiments, and the beneficial effects achieved are also the same as those of the above method embodiments.

[0248] The preferred embodiments of the embodiments of the present application have been described above with reference to the accompanying drawings. This does not limit the scope of the rights of the embodiments of the present application. Any modifications, equivalent replacements, and improvements made by those skilled in the art without departing from the scope and essence of the embodiments of the present application shall fall within the scope of the rights of the embodiments of the present application.

Claims

1. A quadruped robot control method, characterized in that It includes the following steps: Obtain the joint angle information of the quadruped robot, where the joint angle information of the quadruped robot includes the hip adduction / abduction angle, the hip flexion / extension angle, and the knee bending angle; Calculate the rotation error vector according to the joint angle information of the quadruped robot; Construct a model predictive control objective function according to the rotation error vector, and the optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state; Perform quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the foot end force; Generate joint torque commands using the virtual force model control method; Perform robot control according to the predicted value of the foot end force and the joint torque commands.

2. The method according to claim 1, wherein The calculating the rotation error vector according to the joint angle information of the quadruped robot includes: Construct a joint angle vector according to the joint angle information of the quadruped robot; Calculate the corresponding end-effector spatial coordinates of each leg according to the joint angle vector; Calculate the end-effector position combination matrix according to the multiple end-effector spatial coordinates; Calculate the correction matrix according to the correction coefficient and the end-effector position combination matrix; Construct the forward kinematic function corresponding to each leg according to the correction matrix; Take the partial derivative of the forward kinematic function with respect to the joint angle to obtain the Jacobian matrix; Set the yaw angle, pitch angle, and roll angle of the quadruped robot according to the Jacobian matrix; Construct the basic rotation matrix about the yaw angle according to the yaw angle; Construct the basic rotation matrix about the pitch angle according to the pitch angle; Construct the basic rotation matrix about the roll angle according to the roll angle; Combine the basic rotation matrix about the yaw angle, the basic rotation matrix about the pitch angle, and the basic rotation matrix about the roll angle to obtain the overall rotation matrix; Construct an error matrix according to the overall rotation matrix and the desired rotation matrix; Perform logarithmic mapping processing on the overall rotation matrix according to the error matrix; Construct the rotation error vector according to the multiple overall rotation matrices after logarithmic mapping processing.

3. The method according to claim 1, wherein The constructing the model predictive control objective function according to the rotation error vector includes: Construct a translational motion model according to the total mass of the robot, the centroid acceleration, the foot end reaction force, and the gravitational acceleration; Construct a rotational motion model according to the inertia matrix of the robot about the centroid, the rotation error vector, the angular acceleration, the vector from the robot centroid to the leg contact point, and the foot end reaction force; Construct a state vector according to the centroid position, the centroid velocity, the rigid body attitude parameters, and the angular velocity; Construct a translational part dynamic model according to the translational motion model and the state vector; Construct a rotational part dynamic model according to the rotational motion model and the state vector; Construct a state matrix according to the translational part dynamic model; Construct a first input matrix according to the rotational part dynamic model; Construct an overall continuous model according to the state matrix, the first input matrix, the state vector, and the control input vector; Calculate the discretized state transition equation according to the discrete time step and the basic equation of the Euler integration method; Construct a substitution model equation according to the overall continuous model and the discretized state transition equation; Construct a discrete state transition matrix and a control influence matrix according to the substitution model equation; Construct a discrete-time model according to the discrete state transition matrix and the control influence matrix; Perform an approximation process on the translational part according to the discrete-time model to obtain a translational approximation equation set; Perform an approximation process on the rotational part according to the discrete-time model to obtain a rotational approximation equation set; Construct the model predictive control objective function according to the translational approximation equation set, the rotational approximation equation set, the prediction model, and the desired state; 4. The method according to claim 1, wherein Perform a quadratic programming solution process according to the model predictive control objective function to obtain the predicted value of the foot-end force, including: Construct a compression equation according to the foot-end reaction force; Construct a standard quadratic programming objective function according to the compression equation and the model predictive control objective function; Construct foot-end constraints, where the foot-end constraints include foot-end force limit constraints and friction constraints; Combine the standard quadratic programming objective function and the foot-end constraints to obtain a standard quadratic programming problem; Use a quadratic programming solver to solve the standard quadratic programming problem to obtain the predicted value of the foot-end force.

5. The method according to claim 4, characterized in that, The constructing a compression equation according to the foot-end reaction force includes: Set the control quantity according to the foot-end reaction force; Calculate the compensation acceleration according to the centroid position, centroid velocity, desired position, and desired velocity; Calculate the angular acceleration according to the gain matrix of the attitude error, the gain matrix of the angular velocity error, the attitude error vector, the desired angular velocity, and the current angular velocity; Take the compensation acceleration and the angular acceleration as the state tracking target; Construct a friction cone constraint according to the foot-end force and the ground friction coefficient; Perform a linear transformation on the friction cone constraint to obtain a system of linear inequalities, and the system of linear inequalities is used as the energy consumption target; Combine the state tracking target and the friction cone constraint according to the control quantity to obtain a first programming problem, and the parameters of the first programming problem include the foot-end force, the weighted matrix of the state error, and the weighted matrix of the control input; Construct an equation for the relationship between the robot state and the control input according to the first programming problem, the continuous-time system matrix, and the second input matrix; Discretize the equation for the relationship between the robot state and the control input according to the second time step to obtain a discretized model; Perform a compression transformation on the discretized model to obtain the compression equation.

6. The method according to claim 4, characterized in that, The constructing a standard quadratic programming objective function according to the compression equation and the model predictive control objective function includes: Rewrite the model predictive control objective function according to the quadratic penalty of the state error term and the control input term to obtain a first objective function; Rewrite the first objective function according to the compression equation and the desired state to obtain a second objective function; Expand the second objective function to obtain a standard quadratic form objective function; Calculate the overall control quantity according to the inertia matrix of the robot base in the body coordinate system and the overall rotation matrix; Construct a mapping matrix according to the sum of the foot-end forces of each leg, the rotational effect generated by the foot-end forces of each leg, and the skew-symmetric matrix function; Rewrite the standard quadratic form objective function according to the overall control quantity and the mapping matrix to obtain a third objective function; Expand the third objective function to obtain a fourth objective function; Perform a form transformation on the fourth objective function to obtain the standard quadratic programming objective function.

7. The method according to claim 4, wherein The construction of the foot end constraints includes: Construct the foot end force boundary constraint according to the upper limit and lower limit of the foot end force; Construct a friction constraint matrix by using the friction pyramid approximation according to the basic matrix of friction; Construct the friction constraint according to the friction constraint matrix, the upper bound of the vertical component, and the lower bound of the vertical component.

8. The method according to claim 1, wherein The generation of the joint torque command by using the virtual force model control method includes: In horizontal trajectory planning, perform interpolation according to the initial position, the target position, and the cubic Bezier curve interpolation equation to obtain a horizontal trajectory; Differentiate the horizontal trajectory to obtain the horizontal velocity and the horizontal acceleration; In vertical trajectory planning, calculate the velocity and acceleration in the rising stage according to the initial coordinates, the target coordinates, the lifting height, the local time variable in the rising stage, and the cubic Bezier curve interpolation equation; Calculate the velocity and acceleration in the descending stage according to the initial coordinates, the target coordinates, the lifting height, the local time variable in the descending stage, and the cubic Bezier curve interpolation equation; Calculate the expected value of the swing trajectory according to the horizontal velocity, the horizontal acceleration, the velocity in the rising stage, the acceleration in the rising stage, the velocity in the descending stage, and the acceleration in the descending stage, where the expected value of the swing trajectory includes the expected foot end position and the expected foot end velocity; Calculate the corresponding foot end position of each leg in the robot coordinate system according to the centroid position, the expected foot end position, and the fixed offset between the foot end and the geometric reference of the robot body; Calculate the corresponding foot end velocity of each leg in the robot coordinate system according to the centroid velocity and the expected foot end velocity; Calculate the current foot end velocity according to the Jacobian matrix and the joint angular velocity; Calculate the virtual control force according to the corresponding foot end position of each leg in the robot coordinate system, the corresponding foot end velocity of each leg in the robot coordinate system, the current foot end velocity, and the current foot end position; Generate the joint torque command according to the virtual control force.

9. A four-legged robot control device, characterized in that, Includes: A first module for obtaining the joint angle information of the quadruped robot, where the joint angle information of the quadruped robot includes the hip adduction / abduction angle, the hip flexion / extension angle, and the knee bending angle; A second module for calculating the rotation error vector according to the joint angle information of the quadruped robot; A third module for constructing a model predictive control objective function according to the rotation error vector, where the optimization objective of the model predictive control objective function is to minimize the error between the prediction model and the desired state; A fourth module for performing quadratic programming solution processing according to the model predictive control objective function to obtain the predicted value of the foot end force; A fifth module for generating the joint torque command by using the virtual force model control method; A sixth module for controlling the robot according to the predicted value of the foot end force and the joint torque command.

10. A computer device, characterized in that, Includes: At least one processor; At least one memory for storing at least one program; When the at least one program is executed by the at least one processor such that the at least one processor implements the method according to any one of claims 1-8.

Images