A method for dynamic estimation of gait muscle strength based on muscle synergy

Abstract

This invention relates to a gait muscle strength dynamic estimation method based on muscle synergy, belonging to the fields of biomechanical analysis and rehabilitation exoskeleton control technology. The steps include: acquiring gait data of the subject and constructing a personalized musculoskeletal digital twin model; constructing a smooth muscle-tendon dynamic model; constructing an optimal control problem incorporating hierarchical anatomical synergistic constraints; discretizing the problem into a nonlinear programming problem using the direct collocation method; generating an initial solution by implementing a warm-start strategy based on static optimization and a hierarchical degradation mechanism; and solving and outputting the gait cycle muscle strength curve using a solver. This invention, by introducing neural control priors and employing a hierarchical warm-start mechanism, eliminates non-physiological high-frequency computational artifacts while ensuring high dynamic tracking accuracy, and solves the problem of getting trapped in local minima in high-dimensional solutions. Therefore, it provides a high-fidelity and robust methodological basis for clinical rehabilitation and exoskeleton controller planning for specific patients.

Description

Technical Field

[0001] This invention relates to the field of biomechanical analysis and rehabilitation wearable robot control technology, and in particular to a method for dynamic estimation of gait muscle strength based on muscle coordination. Background Technology

[0002] Gait analysis, by quantifying the kinematic and dynamic characteristics of human walking, reveals the control strategies of the neuromuscular system and is a core research method in biomechanics, clinical rehabilitation engineering, and wearable robot (such as rehabilitation exoskeletons) control. However, due to ethical constraints and technological bottlenecks, the contractile force of deep human muscles is difficult to measure directly using non-invasive methods. Furthermore, the human lower limb motor system exhibits significant "muscle redundancy," meaning the number of muscle actuators driving the joints far exceeds the number of degrees of freedom in the system. This biological characteristic leads to inherent ill-posedness in inverse dynamics analysis; the force distribution scheme of a single muscle cannot be uniquely determined solely based on observed trajectories and external reaction forces. Therefore, while ensuring physiological rationality, analyzing the dynamic distribution of muscle forces during the gait cycle using computational models and optimization theory remains a crucial problem that urgently needs to be addressed in this field.

[0003] Currently, optimization methods are mainly used to search for optimal muscle recruitment strategies, including static optimization and traditional dynamic optimization (such as computational muscle control, CMC). Firstly, static optimization methods are widely used due to their simplicity, but they assume that muscle force depends only on instantaneous torque requirements, ignoring the time dependence of muscle contraction dynamics, resulting in low prediction accuracy in periodic movements such as gait. Secondly, while traditional dynamic optimization methods (such as CMC) consider the time-varying characteristics of muscles, they often overfit the inverse dynamic residuals, leading to non-physiological high-frequency oscillations (computational artifacts) in the prediction results. This phenomenon stems from the fact that simple mechanical goal optimization ignores the inherent "muscle synergy" mechanism of the central nervous system—the physiological fact that the nervous system controls muscle groups through modular signals—causing the algorithm to tend to produce "on-off" activation patterns that do not conform to human physiology. Finally, in high-fidelity musculoskeletal modeling applications specific to individual subjects, differences in individual anatomical parameters lead to high numerical rigidity in the optimization problem. This makes the constructed large-scale nonlinear programming (NLP) problem extremely sensitive to initial guesses, extremely difficult to solve, and prone to getting trapped in local minima or causing the algorithm to fail to converge. In summary, existing gait muscle strength calculation frameworks still have significant limitations in integrating physiological prior knowledge and solving the computational robustness of high-dimensional models, making it difficult to directly provide high-fidelity and high-robust underlying algorithmic support for clinical gait abnormality assessment and personalized torque planning of exoskeleton controllers. Summary of the Invention

[0004] This invention aims to address the technical problems in existing technologies, such as the generation of non-physiological high-frequency computational artifacts in traditional dynamic optimization algorithms when analyzing muscle redundancy, and the high numerical rigidity, easy trapping in local minima, or convergence failure in personalized high-fidelity musculoskeletal models during nonlinear solutions. The invention provides a gait muscle strength dynamic estimation method based on muscle synergy.

[0005] To solve the above-mentioned technical problems, the technical solution of the present invention is as follows:

[0006] A gait muscle strength dynamic estimation method based on muscle coordination, whose data acquisition device includes reflective markers on the body surface 1, a three-dimensional force platform 2, an infrared optical motion capture camera 3, and a host processor 4, performs the following steps based on the acquired data:

[0007] Step S1: Obtain the subject's gait data and construct a personalized musculoskeletal digital twin model;

[0008] Step S2: Construct a smooth muscle-tendon dynamics model and establish the core physical path constraints of the state variables;

[0009] Step S3: Construct the optimal control problem with fusion hierarchical anatomical co-constraints;

[0010] Step S4: Numerical discretization based on the direct collocation method;

[0011] Step S5: Implement a hierarchical hot start strategy based on static optimization priors;

[0012] Step S6: Solve using a nonlinear solver and output digital twin muscle strength index.

[0013] In the above technical solution, step S1, acquiring the subject's gait data and constructing a personalized musculoskeletal digital twin model, specifically includes:

[0014] The gait motion information of the subject under a specific walking pattern is collected synchronously by a data acquisition device. The motion information includes three-dimensional kinematic trajectory and ground reaction force data.

[0015] Based on a general human musculoskeletal model, the collected data is used to perform geometric scaling to adjust the optimal muscle fiber length, tendon relaxation length, and maximum isometric muscle strength in order to construct a personalized musculoskeletal model.

[0016] By combining inverse kinematics and inverse dynamics calculations, the joint angles, net torque requirements, and physical boundary conditions of the muscle torque arms of each joint in the lower limb during the gait cycle are obtained.

[0017] In the above technical solution, step S2, constructing a smooth muscle-tendon dynamics model, specifically includes:

[0018] The mechanical properties of each muscle are described using a Hill-type three-element muscle-tendon unit model;

[0019] By introducing a nonlinear hyperbolic tangent function to smooth the activation and deactivation time constants, a continuously differentiable activation dynamics differential equation is constructed.

[0020] Based on the continuous exponential and Gaussian functions to express the active force-length, passive force-length, and force-velocity characteristics of muscles, a differential equation for muscle contraction dynamics is constructed to form path constraints for the dynamic optimization problem.

[0021] In the above technical solution, step S3, which involves constructing an optimal control problem with fused hierarchical anatomical collaborative constraints, specifically includes:

[0022] The muscle force distribution during the gait cycle is set as a multi-objective optimal control problem in the time domain;

[0023] The objective function of the optimal control problem is set to consist of weighted sub-objectives, specifically including: a metabolic energy consumption minimization term, a reserve torque minimization term, a motion smoothness regularization term, and an anatomical hierarchical cooperative constraint term;

[0024] The anatomical hierarchical synergistic constraint term includes a three-level penalty mechanism: anatomical bundle penalty, same-function agonist penalty, and functional synergistic pair penalty;

[0025] By explicitly introducing the anatomical hierarchical co-containment constraint term into the objective function, functionally related muscle groups are forced to exhibit a coordinated co-activation state in the form of soft constraints, thus simulating the modular control strategy of the central nervous system.

[0026] In the above technical solution, step S4 specifically includes:

[0027] The optimal control problem is discretized in a fully implicit manner using the direct collocation method based on Legendre-Gauss-Radau nodes;

[0028] By using Lagrange polynomials to approximate the state variables, the continuous differential equations are transformed into algebraic constraints on discrete nodes, thus converting the optimal control problem into a large-scale nonlinear programming problem.

[0029] In the above technical solution, step S5 specifically includes:

[0030] A quadratic programming algorithm is used to solve the static optimization problem with simplified cooperative constraints frame by frame;

[0031] During the frame-by-frame solution process, an automatic hierarchical degradation mechanism is set up and implemented, specifically as follows:

[0032] Prioritize performing fully constrained calculations that are subject to both torque balance and muscle physiological extremes;

[0033] If the current frame encounters a situation where there is no solution due to constraint conflicts, a downgrade is triggered, and the unidirectional constraints are relaxed to retry the calculation.

[0034] If no solution is found, a further degradation is triggered, and calculations are performed while only retaining the torque balance constraints.

[0035] The robustly obtained static optimization discrete solution is then smoothed by splines and interpolated onto a dynamically optimized collocation grid to generate physically feasible, high-quality initial guesses that are then passed to the nonlinear solver.

[0036] In the above technical solution, step S6 specifically includes:

[0037] Load the initial guess value, and call the nonlinear solver based on the interior point method to iteratively optimize the large-scale nonlinear programming problem until the preset convergence tolerance is reached;

[0038] Output and save the temporally continuous lower limb muscle force distribution curves and muscle activation levels that satisfy physiological boundaries within the gait cycle.

[0039] The present invention has the following beneficial effects:

[0040] This invention innovatively transforms the prior knowledge of "muscle synergy" in the central nervous system into a mathematical regularization term in optimal control, supplemented by direct collocation technology across the entire time domain. This method effectively filters out non-physiological high-frequency computational artifacts (oscillations) caused by overfitting to inverse dynamic residuals in traditional muscle control computation methods. The generated muscle recruitment pattern is more coordinated, and the peak muscle force distribution is more consistent with the real physiological mechanism of distributed load in the human body. Addressing the numerical rigidity problem commonly found in high-fidelity personalized digital twin models, this invention implements a hot-start strategy based on static optimization and automatic hierarchical degradation mechanisms. This strategy provides an absolutely physically feasible search starting point for high-dimensional nonlinear solvers, completely breaking through the computational bottleneck of traditional algorithms that are prone to getting trapped in local minima or failing to converge. While ensuring extremely high dynamic tracking accuracy (minimal joint torque reconstruction error), this invention significantly improves the physiological fidelity of predicted muscle force, thus providing high-fidelity and highly robust underlying algorithmic support for clinical rehabilitation assessment of specific patients, abnormal gait analysis, and personalized torque planning for wearable exoskeletons. Attached Figure Description

[0041] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0042] Figure 1 This is a schematic diagram of the overall system architecture of the gait muscle strength dynamic estimation method based on muscle coordination of the present invention.

[0043] Figure 2This is a schematic diagram illustrating the steps involved in constructing a personalized musculoskeletal model in the gait muscle strength dynamic estimation method based on muscle synergy of the present invention.

[0044] Figure 3 This is a schematic diagram illustrating the steps of the layered hot-start strategy in the gait muscle strength dynamic estimation method based on muscle synergy of the present invention.

[0045] Figure 4 This is a schematic diagram illustrating the steps of multi-level collaborative optimization in the gait muscle strength dynamic estimation method based on muscle coordination of the present invention.

[0046] Figure 5 This is a schematic diagram of the gait data acquisition device involved in the gait muscle strength dynamic estimation method based on muscle coordination of the present invention. Detailed Implementation

[0047] The core inventive idea of ​​this invention is:

[0048] The present invention provides a gait muscle strength dynamic estimation method based on muscle synergy. This method explicitly introduces the modular control characteristics of "muscle synergy" of the central nervous system as a regularization penalty term into the objective function of the optimal control problem (OCP). It then uses a direct collocation method based on Legendre-Gauss-Radau (LGR) nodes for full-time domain discretization. Finally, it designs a static optimization hot-start strategy based on a hierarchical degradation mechanism to overcome the numerical rigidity problem in solving large-scale nonlinear programming (NLP) problems, providing high-quality initial guesses for high-dimensional solvers. This achieves personalized, high physiological fidelity, and highly robust lower limb gait muscle strength dynamic estimation and prediction.

[0049] The present invention will now be described in detail with reference to the accompanying drawings.

[0050] like Figures 1 to 5 As shown, the gait muscle strength dynamic estimation method based on muscle coordination of the present invention has a data acquisition device including reflective markers on the body surface 1, a three-dimensional force platform 2, an infrared optical motion capture lens 3, and a host processor 4. Based on the acquired data, the following steps are performed:

[0051] Step S1: Obtain the subject's gait data and construct a personalized musculoskeletal digital twin model.

[0052] Personalized kinematic and kinetic data of the subjects were collected using a data acquisition device.

[0053] In this embodiment, combined with Figure 5The diagram shows a gait data acquisition device. The subject wears a special test suit with multiple reflective markers (1) on the body surface and walks at a constant speed on a test track with a three-dimensional force platform (2). Infrared optical motion capture lenses (3) installed around the space track the three-dimensional spatial coordinates of the reflective markers in real time; at the same time, the three-dimensional force platform (2) collects ground reaction force (GRF) and torque information.

[0054] Specifically, after the above data is transmitted to the workstation, it is smoothed and denoised using a fourth-order zero-phase-shift Butterworth low-pass filter with a cutoff frequency of 6Hz. The filtering process is calculated according to the following formula:

[0055]

[0056] Where x(n) represents the original acquired data, y(n) represents the filtered data, and b0, b1, b2, a1, a2 represent the filter coefficients determined according to the set cutoff frequency.

[0057] First, combined Figure 2 As shown, after filtering, a general OpenSim musculoskeletal model (such as the Gait2392 general musculoskeletal model) is imported, and geometric scaling is performed. The formula for proportional scaling of the optimal muscle fiber length is:

[0058]

[0059] in, This represents the optimal muscle fiber length for the personalized model. D represents the optimal muscle fiber length of the general model. subject and D generic These represent the distances between the subject and the corresponding skeletal segment markers in the general model.

[0060] Then, allometric growth scaling calibration was performed based on the subject's body weight to determine maximum isometric muscle strength:

[0061]

[0062] in, M represents the calibrated maximum isometric muscle strength. subject M represents the subject's weight. generic This indicates the standard weight set by the general model.

[0063] Next, inverse kinematics (IK) solutions are performed based on the constructed personalized model, followed by inverse dynamics (ID) solutions. The Newton-Euler equations for the inverse dynamics multibody system can be expressed as:

[0064]

[0065] in, Let M(q) represent the net torque vector of each joint, and M(q) represent the system mass matrix. Let G(q) represent the Coriolis force and centrifugal force terms, and let R(q) represent the gravity term. GRF The torque q represents the force generated by the ground reaction force. , These represent the angle, angular velocity, and angular acceleration of the joint, respectively.

[0066] Step S2: Construct a smooth muscle-tendon dynamics model and establish the core physical path constraints of the state variables.

[0067] Combination Figure 3 As shown, this embodiment uses the classic Hill-type three-element muscle-tendon unit model, and the mechanical equation for muscle generation is:

[0068]

[0069] Among them, F M 'a' represents muscle strength, 'a' represents muscle activation, and 'f' represents muscle strength. act Main power-length coefficient, f v The force-velocity coefficient, f pas The passive-length coefficient, and These are normalized muscle length and velocity, respectively. It is a feathery angle.

[0070] Specifically, to overcome the numerical discontinuity caused by the abrupt change in activation and deactivation time constants in traditional models, a nonlinear hyperbolic tangent function (tanh) is introduced to construct a smooth and differentiable activation dynamics equation:

[0071]

[0072]

[0073] Where e represents the neural activation signal, and a represents the muscle activation level. and These represent the time constants for activation and deactivation, with values ​​of 0.015s and 0.060s, respectively.

[0074] First, the tendon force-length property is described in exponential form as its nonlinear elasticity, calculated according to the following formula:

[0075]

[0076] in, k represents the normalized tendon length. TThe stiffness coefficient is 35, and c1, c2, and c3 represent constant coefficients of 0.200, 0.995, and 0.250, respectively.

[0077] Then, the active force-length characteristic is described using a continuous Gaussian function model and calculated according to the following formula:

[0078]

[0079] in, b represents the normalized fiber length. 1i b 2i b 3i b 4i This represents the Gaussian curve smoothing coefficient.

[0080] Then, the passive-length characteristic is described by an exponential function and calculated according to the following formula:

[0081]

[0082] Where, k pe e0 represents the exponential shape parameter with a value of 4.0, and e0 represents the passive strain limit with a value of 0.6.

[0083] Finally, the force-velocity characteristics are described using a hyperbolic logarithmic model and calculated according to the following formula:

[0084]

[0085] in, d1 represents the normalized fiber velocity, and d2, d3, and d4 represent the smoothing coefficients of the force-velocity curve, which are -0.318, -8.149, -0.374, and 0.886, respectively.

[0086] Step S3: Construct the optimal control problem with fusion hierarchical anatomical co-constraints;

[0087] Combination Figure 4 As shown, the muscle force distribution during the gait cycle is set as a multi-objective optimal control problem (OCP), aiming to simulate the coordinated control strategy of the central nervous system.

[0088] Specifically, the objective function J in the continuous time domain is calculated according to the following formula:

[0089]

[0090] Where t0 and t f The start and end times of the gait are represented by w1 to w4, which represent the adaptive weighting coefficients of each sub-objective.

[0091] In this embodiment, J1 is the metabolic energy minimization term, defined as the sum of squares of all muscle activations:

[0092]

[0093] Where, N mus a represents the total number of muscles in the model. m (t) represents the activation level of the m-th muscle, w m This represents the metabolic weighting coefficient for the corresponding muscle.

[0094] J2 is the reserve torque minimization term, used to compensate for model measurement errors and inverse dynamic residuals:

[0095]

[0096] Where, N dof Indicates the number of degrees of freedom of the model. T represents the reserve torque of the j-th degree of freedom. opt This represents the set ideal reserve peak torque.

[0097] J3 is the control signal smoothness regularization term, which effectively suppresses high-frequency computational oscillations by penalizing the time derivatives of variables such as activation degree.

[0098]

[0099] in, and w represents the time derivative of muscle activation and normalized tendon strength, respectively. a and w F This represents the smoothness regularization weight.

[0100] J4 is the anatomical layered synergistic constraint term, which is the core innovation of this invention. It includes anatomical bundles, agonist muscles of the same function, and functional synergistic pair constraints.

[0101]

[0102] Among them, S parts S agonist and S func These represent the three levels of anatomical sets, w parts w agonist wfunc represents the corresponding hierarchical penalty weight, a k，1 and a k，2 This indicates the activation level of muscles belonging to the same functional synergistic group.

[0103] Step S4: Numerical discretization based on the direct collocation method;

[0104] Combination Figure 4As shown, this embodiment uses Legendre-Gauss-Radau (LGR) nodes to implicitly discretize the time domain, transforming the infinite-dimensional optimal control problem into a high-dimensional nonlinear programming (NLP) problem.

[0105] First, the time domain Mapping to collocation interval :

[0106]

[0107] Where t represents physical time and τ represents the standard time of the LGR coordinate interval.

[0108] Then, the state variables are approximated using Lagrange polynomials, and the state differential equation is transformed into an algebraic equation at the collocation point:

[0109]

[0110] Where N represents the number of coordinate intervals, and D k，j Describes the differential matrix elements of the LGR collocation method. and f represents the values ​​of the state variable and the control variable at the discrete nodes, respectively. i This represents the system's dynamic function.

[0111] Step S5: Implement a hierarchical hot start strategy based on static optimization priors;

[0112] To address the high numerical rigidity of large-scale nonlinear programming, combined with Figure 3 As shown, this embodiment performs frame-by-frame layered static optimization (SO) to generate initial guess values.

[0113] First, the first level is fully constrained optimization: execution is subject to joint torque balance and muscle physiological extremes. Optimization with dual constraints. If the calculation converges, the result is saved; if there is no solution due to kinematic incompatibility, a first-level degradation is triggered.

[0114] Then, the second level is relaxation constraint optimization: relaxation variables are introduced to relax the physiological extreme value constraint to... A large penalty term is added to the objective function to force optimization. If the level still fails to converge, a second-level fallback solution is triggered.

[0115] Then, the third level is a fallback solution of torque balance only: directly removing the physiological extreme value constraints of muscles and retaining only the torque balance equation. In this state, there must be a mathematically feasible solution with a global minimum.

[0116] In the above embodiment, after all gait frames are calculated, the discrete solution is smoothed by Cubic Spline interpolation to generate high-quality hot-start initial guess values ​​aligned to the LGR collocation grid.

[0117] Step S6: Solve the problem using a nonlinear solver and output the digital twin muscle strength index;

[0118] The transformed large-scale nonlinear programming (NLP) model and the initial guesses are loaded into the interior-point solver (IPOPT) based on the CasaADi framework for iteration.

[0119] Specifically, IPOPT uses the primal-dual interior-point method to search for solutions that satisfy the Karush-Kuhn-Tucker (KKT) optimality conditions, and its error verification criterion is calculated according to the following formula:

[0120]

[0121] in, Represent the Lagrange function, Indicates the degree of constraint violation. This indicates the preset convergence tolerance. After iteration stops, the output is a sequence of dynamic muscle strength and activation levels for each muscle in the lower limbs, after eliminating non-physiological high-frequency oscillations.

[0122] This invention presents a gait muscle strength dynamic estimation method based on muscle synergy. It employs an infrared optical motion capture system and a three-dimensional force platform to collect kinematic and dynamic data from subjects, constructing a personalized musculoskeletal digital twin model. The modular control characteristics of the central nervous system's "muscle synergy" are explicitly introduced as a regularization penalty term into a multi-objective optimal control problem. Then, a direct collocation method based on Legendre-Gauss-Radau (LGR) nodes is used to discretize the continuous problem into a large-scale nonlinear programming model. Finally, when solving using a nonlinear solver, a static optimization warm-start strategy based on a hierarchical degradation mechanism is designed and implemented to generate high-quality initial guesses and complete iterative calculations, obtaining the subject's dynamic muscle strength and muscle activation sequence. This method differs from traditional dynamic optimization methods such as computational muscle control (CMC). While ensuring extremely high joint torque tracking accuracy, it effectively eliminates the high-frequency computational artifacts (oscillations) caused by overfitting to inverse dynamic residuals, which is common in traditional methods. By breaking through the numerical rigidity of high-dimensional personalized models at the underlying mathematical level, this invention not only significantly improves the convergence robustness of the full-time domain solution, but also makes the predicted muscle recruitment sequence smoother and highly consistent with the real neurophysiological laws. This provides high-fidelity and high-reliability underlying algorithm support for clinical abnormal gait analysis, rehabilitation effect evaluation, and personalized torque planning of exoskeleton robots.

[0123] Obviously, the above embodiments are merely illustrative examples for clear explanation and are not intended to limit the implementation. Those skilled in the art will recognize that other variations or modifications can be made based on the above description. It is neither necessary nor possible to exhaustively list all possible implementations here. However, obvious variations or modifications derived therefrom are still within the scope of protection of this invention.

Claims

1. A method for dynamic estimation of gait muscle strength based on muscle synergy, characterized in that, Based on the data acquisition device, which includes reflective markers on the body surface, a three-dimensional force platform, an infrared optical motion capture camera, and a host processor, the following steps are performed: Includes the following steps: Step S1: Obtain the subject's gait data and construct a personalized musculoskeletal digital twin model; Step S2: Construct a smooth muscle-tendon dynamics model and establish the core physical path constraints of the state variables; Step S3: Construct the optimal control problem with fusion hierarchical anatomical co-constraints; Step S4: Numerical discretization based on the direct collocation method; Step S5: Implement a hierarchical hot start strategy based on static optimization priors; Step S6: Solve using a nonlinear solver and output digital twin muscle strength index.

2. The gait muscle strength dynamic estimation method based on muscle synergy according to claim 1, characterized in that, In step S1, gait data of the subjects is acquired and a personalized musculoskeletal digital twin model is constructed, specifically including: Personalized kinematic and dynamic data of the subjects were acquired simultaneously using a motion capture system and a force platform, and the original marker trajectory and ground reaction force data were smoothed and denoised using a fourth-order zero-phase-shift Butterworth low-pass filter. Import a general musculoskeletal model, use the kinematic data to geometrically scale the general musculoskeletal model, adaptively adjust the optimal muscle fiber length and tendon relaxation length, and calibrate the maximum isometric muscle strength according to the subject's weight to construct a personalized musculoskeletal digital twin model. Based on the personalized musculoskeletal digital twin model, inverse kinematics and inverse dynamics calculations are performed to extract the torque arm, tendon length, and contraction velocity of each muscle during the gait cycle as state boundary conditions.

3. The gait muscle strength dynamic estimation method based on muscle synergy according to claim 1, characterized in that, In step S2, a smooth muscle-tendon dynamics model is constructed, specifically including: The mechanical properties of lower limb muscles were described using a Hill-type three-element muscle-tendon unit model. A nonlinear hyperbolic tangent function is introduced to smoothly transition the muscle activation time constant and deactivation time constant, and a continuously differentiable activation kinetic differential equation is constructed. The active-length characteristic, passive-length characteristic, and force-velocity characteristic of muscle contractile elements are expressed using continuous exponential functions and Gaussian functions, respectively.

4. The gait muscle strength dynamic estimation method based on muscle synergy according to claim 1, characterized in that, In step S3, the optimal control problem with fusion hierarchical anatomical co-constraints is constructed, specifically including: The muscle force distribution during the gait cycle is set as a multi-objective optimal control problem in the time domain. Its objective function consists of four weighted sub-terms: a metabolic energy consumption minimization term, a reserve torque minimization term, a control signal smoothness regularization term, and an anatomical hierarchical cooperative constraint term. The metabolic energy consumption minimization term is simulated by minimizing the sum of squares or cubes of all muscle activations; the control signal smoothness regularization term suppresses high-frequency jumps by penalizing the time derivatives of muscle activation and normalized tendon force. The anatomical hierarchical synergistic constraint terms include anatomical bundle penalties, agonist muscle penalties, and antagonistic or synergistic function pair penalties, which force functionally related muscle groups to exhibit a coordinated co-activation state through soft constraints.

5. The gait muscle strength dynamic estimation method based on muscle synergy according to claim 1, characterized in that, In step S4, the numerical discretization process based on the direct collocation method is as follows: Implicit discretization of the time grid is performed using the direct collocation method based on Legendre-Gauss-Radau nodes, mapping the time domain to the standard collocation interval; By using Lagrange polynomials to approximate the continuous trajectories of state and control variables, the continuous state differential equation is transformed into algebraic path constraints at discrete collocation points, thereby converting the optimal control problem into a large-scale nonlinear programming model.

6. The gait muscle strength dynamic estimation method based on muscle synergy according to claim 1, characterized in that, In step S5, a hierarchical hot start strategy is implemented based on static optimization priors, specifically including: Perform static optimization calculations with simplified cooperative constraints frame by frame for the gait cycle; When calculating the current frame, the automatic hierarchical degradation mechanism is executed sequentially: Prioritize performing fully constrained static optimization calculations that are subject to both joint torque balance and muscle physiological extremes; If the solution does not converge, a first-level degradation mechanism is triggered, which introduces relaxation variables to relax the boundary of muscle physiological extreme values ​​and performs static optimization calculation with relaxed one-way constraints. If convergence is still not achieved, a second-level degradation mechanism is triggered, the muscle physiological extreme value constraint is removed, and a static optimization calculation is performed that retains only the torque balance equation. After all gait frames have been calculated, the summarized static optimization discrete solutions are smoothed by spline interpolation and mapped onto the collocation grid to generate initial guess values ​​for hot start.

7. The gait muscle strength dynamic estimation method based on muscle synergy according to claim 1, characterized in that, Step S6 is as follows: The initial guess value for hot start generated in step S5, the boundary conditions extracted in step S1, and the nonlinear programming model constructed in step S4 are loaded into the nonlinear programming solver based on the interior point method for iterative optimization calculation. The nonlinear programming solver based on the interior point method uses the primal-dual interior point method to search for solutions that satisfy the Karush-Kuhn-Tucker optimality condition. The iteration stops when the descent gradient of the objective function and the constraint violation meet the preset convergence tolerance accuracy, and the dynamic muscle strength curves and activation sequences of each muscle in the lower limbs are output after eliminating non-physiological high-frequency oscillations.

Images