Online optimal control method for nonlinear MIMO (Multiple Input Multiple Output) system based on approximate dynamic programming

Abstract

The invention relates to the technical field of on-line optimal control, and discloses a nonlinear MIMO system on-line optimal control method based on approximate dynamic programming, which comprises the following steps: step 1, obtaining and preprocessing states, instructions, errors, constraints, power and chassis postures of double mechanical arms; 2, generating a constraint indication sequence, and calculating a switching rate and a lag amount to obtain a constraint alternating judgment result; step 3, generating an anti-saturation judgment result and scheduling a weight; 4, obtaining a sensitivity and linearization coefficient matrix based on the input and state increment; step 5, combining the scheduling weight and the matrix construction stage cost and superposing future cost approximation, and solving an unlimited and actual instruction; step 6, using residual recursion least square to update future cost model parameters; and 7, when the judgment result is lower than the threshold value, recovering the weight to the reference value. According to the method, real-time optimal control and high-stability operation of the nonlinear MIMO system under the constraint condition are realized.

Description

Technical Field

[0001] This invention relates to the field of online optimal control, and more specifically, to an online optimal control method for nonlinear MIMO systems based on approximate dynamic programming. Background Technology

[0002] Nonlinear multiple-input multiple-output (MIMO) systems are widely used in industrial robots, collaborative robotic arms, and complex electromechanical equipment. For collaborative operation of two robotic arms, the system needs to process multiple inputs and state variables simultaneously, which often exhibit strong coupling and significant nonlinearity. During operation, the two robotic arms must achieve precise collaborative trajectory tracking while adhering to power limits and actuator constraints. However, due to the dynamic coupling between the two robotic arms, the phenomenon of alternating activation and deactivation of constraints frequently occurs. Furthermore, when one robotic arm experiences input saturation, it can affect the other robotic arm through load transfer, creating a risk of cross-arm saturation. These situations make it difficult to guarantee the stability and safety of the control process.

[0003] Dynamic programming theory can theoretically derive optimal control strategies, but its direct application faces limitations such as excessive computational cost and inability to meet real-time requirements due to the high dimensionality of the state and input spaces. To alleviate this problem, approximate dynamic programming methods have been proposed. The basic idea is to construct an optimization objective using an approximate model of stage costs and future costs, and then iteratively approach the optimal solution. However, existing methods mostly involve offline modeling or fixed parameter updates, making it difficult to adaptively adjust to changes in operating conditions. For dual-arm control, these methods fail to reflect alternating constraint changes in a timely manner and lack mechanisms to handle cross-arm saturation transmission. Summary of the Invention

[0004] This invention provides an online optimal control method for nonlinear MIMO systems based on approximate dynamic programming, which solves the technical problems in related technologies such as the excessive computational load of dynamic programming making it difficult to apply in real time, the inability of the future cost model to adaptively update with changes in operating conditions, and the difficulty in effectively dealing with the risks of constraint alternation and crossarm saturation.

[0005] This invention provides an online optimal control method for nonlinear MIMO systems based on approximate dynamic programming, comprising the following steps: Step 1: Obtain the system status, control commands, alignment error, constraint values, power data, and mobile chassis attitude of the dual flexible robotic arms, and perform preprocessing to obtain preprocessed data; Step 2: Generate constraint active indication sequences for the two robotic arms based on the preprocessed data, calculate the alternation switching rate and phase lag, and obtain the constraint alternation determination result; Step 3: Obtain the saturation shear amount based on the difference between the unlimited control command and the actual control command in the previous control cycle. Combined with the power increment of the other robotic arm, calculate the cross-arm load transfer coefficient and lag time to obtain the anti-saturation judgment result; schedule the constraint penalty weight, constraint bandwidth and power weight. Step 4: Within the sliding time window, obtain the local sensitivity matrix and linearization coefficient matrix based on the input increment and state increment of the preprocessed data; Step 5: Based on the scheduling constraint penalty weight, constraint bandwidth and power weight, as well as the local sensitivity matrix and linearization coefficient matrix, construct the stage cost function and superimpose the future cost approximation model. Under the premise of satisfying the input limiting and linearization constraint model, solve to obtain the unlimited control command and generate the actual control command. Step 6: Based on the current stage cost, the residual between the discounted approximate future cost of the next control period and the approximate future cost of the current control period, recursive least squares is used to update the parameters of the future cost approximation model online. Step 7: When the constraint alternation judgment result and the anti-saturation judgment result are lower than the corresponding threshold, restore the constraint penalty weight, constraint bandwidth and power weight to the baseline value.

[0006] Furthermore, the generation of the constraint alternation determination result includes: Step 11: For the two robotic arms, at the current sampling time, take the maximum value of each constraint value and compare it with the preset constraint judgment threshold. If the maximum value is greater than or equal to the preset constraint judgment threshold, generate a constraint activity indicator value of 1 for the corresponding robotic arm; otherwise, generate a constraint activity indicator value of 0. Step 12: Within the sliding time window, the constraint activity indicator values ​​of the continuously sampled robotic arm are arranged sequentially to form the constraint activity indicator sequence of the robotic arm. Step 13: Within the sliding time window, count the number of state transitions of the constraint active indication sequence of the two robotic arms, and use the ratio of the number of state transitions to the length of the sliding time window as the alternation rate. Step 14: Within the sliding time window, calculate the normalized cross-correlation function of the constraint active indication sequence of the two robotic arms, and determine the time shift corresponding to the maximum value of the cross-correlation function as the phase lag. Step 15: When the alternation switching rate is greater than the preset switching rate threshold and the phase lag is within the preset lag interval, the constraint alternation judgment result is 1; otherwise, it is 0, which respectively indicate that alternation exists and alternation does not exist.

[0007] Furthermore, the process of generating the anti-saturation determination result includes: Step 21: At the current sampling time, calculate the difference vector between the unlimited control command and the actual control command of the tested robotic arm, and take the Euclidean norm of the difference vector as the saturation shear quantity; calculate the power value difference of the other robotic arm at adjacent sampling times as the power increment, and perform unified timestamp alignment within the sliding time window. Step 22: Within the sliding time window, calculate the normalized cross-correlation function of the sequences corresponding to the saturation shear rate and the power increment, and determine the minimum positive time shift corresponding to the maximum value of the normalized cross-correlation function as the lag time. Step 23: Under the lag time, the least squares regression method is used to calculate the crossarm load transfer coefficient between the saturation shear and the power increment; Step 24: When the saturation shear is greater than the preset amplitude threshold, the crossarm load transfer coefficient is greater than the preset coefficient threshold, and the lag time is within the preset time shift interval, the anti-saturation judgment result is output as having crossarm saturation transfer risk; otherwise, the output is no risk.

[0008] Furthermore, scheduling mapping is performed on the constraint penalty weights, constraint bandwidth, and power weights, including: Step 31: When the constraint alternation judgment result is 1, the alternation switching rate is normalized according to the preset switching rate threshold to obtain the alternation risk intensity; when the anti-saturation judgment result is that there is a risk of cross-arm saturation transmission, the cross-arm load transmission coefficient is normalized according to the preset coefficient threshold to obtain the saturation risk intensity. Step 32: Based on alternating risk intensity and saturation risk intensity, use monotonic linear functions to adjust the constraint penalty weight, constraint bandwidth and power weight respectively, and limit them within the preset upper and lower bounds. Step 33: Update the adjusted constraint penalty weights, constraint bandwidth, and power weights using an exponential weighting method to obtain the scheduled constraint penalty weights, constraint bandwidth, and power weights.

[0009] Furthermore, within the sliding time window, based on the input increment and state increment of the preprocessed data, the local sensitivity matrix and linearization coefficient matrix are obtained, including: Step 41: Within the sliding time window, the difference between the current sampled input and the previous sampled input is used as the input increment, and the difference between the current sampled state and the previous sampled state is used as the state increment. The input increment sequence and the state increment sequence are then aligned one by one according to the unified timestamp. Step 42: Based on the concatenation of the input increment sequence and the state increment sequence to form a regression matrix, the least squares method is used to perform regression estimation to obtain the linearized coefficient matrix; Step 43: Limit the numerical magnitude of the linearized coefficient matrix, and trigger consistency verification and truncation when the matrix norm exceeds the threshold; Step 44: Determine the mapping coefficient matrix from input increment to state increment in the linearization coefficient matrix as the local sensitivity matrix, and update it using an exponential weighting method to obtain the local sensitivity matrix of the current control cycle.

[0010] Furthermore, the process of constructing the stage cost function and superimposing the future cost approximation model includes: Step 51: Obtain the state deviation, input increment, and constraint relaxation amount of the current control cycle. The state deviation is the difference between the current system state and the preset reference trajectory, the input increment is the difference between the current input and the input of the previous control cycle, and the constraint relaxation amount is the excess amount of the linearized constraint exceeding the preset bandwidth. Step 52: Construct a stage cost function based on the state deviation term, input increment term, and constraint relaxation term. The state deviation term is obtained by weighting the state deviation with a preset weight matrix. The input increment term is obtained by weighting the power weight with the input penalty weight. The constraint relaxation term is obtained by multiplying the square of the constraint relaxation amount with the constraint penalty weight. Step 53: The approximate future cost model determined in the previous control cycle is superimposed onto the stage cost function according to a preset discount factor to form the overall objective function of the current control cycle. The approximate future cost model is a parameterized function, and its parameters are updated online by the recursive least squares method.

[0011] Furthermore, the generation of unlimited control commands and actual control commands includes: Step 61: Based on the overall objective function, apply upper and lower limit constraints to the input, and use the linearization coefficient matrix and constraint bandwidth to perform a first-order expansion of the nonlinear constraint function at the current state and input point to form a linearized constraint model. Step 62: Based on the overall objective function, construct and solve the convex quadratic programming problem to obtain the optimal input increment for the current control cycle; Step 63: Add the optimal input increment to the input of the previous control cycle to obtain the unlimited control command; Step 64: Compare the unlimited control command with the input upper and lower limits. If it exceeds the upper limit, take the upper limit; if it is below the lower limit, take the lower limit to obtain the actual control command.

[0012] Furthermore, the parameters of the future cost approximation model are updated online during the control operation using a recursive least squares method. At the end of each control cycle, a consistency error signal is constructed based on the difference between the sum of the stage cost and the discounted future cost approximation of the next control cycle and the future cost approximation of the current control cycle. This error signal is then used as the driving input to recursively correct the parameter vector of the future cost approximation model, while simultaneously recursively updating the covariance matrix.

[0013] Furthermore, when the constraint alternation judgment result and the anti-saturation judgment result are lower than the corresponding thresholds, the constraint penalty weight, constraint bandwidth, and power weight are restored to the baseline values, including: Step 71: Detect alternating risk intensity and saturation risk intensity within the sliding time window. Trigger recovery operation when both are continuously lower than their respective preset thresholds for a preset number of times. Step 72: Call the preset baseline values ​​and upper and lower bound parameters of the system calibration stage, and read the baseline values ​​and corresponding safety boundaries of constraint penalty weight, constraint bandwidth and power weight respectively; Step 73: When recovery is triggered, perform a linear weighted average of the current constraint penalty weight, constraint bandwidth and power weight according to the preset regression rate, so that it regresses to the corresponding benchmark value.

[0014] The beneficial effects of this invention are as follows: It can simultaneously ensure real-time performance, constraint safety, and control performance under complex operating conditions; by introducing constraint alternation and anti-saturation judgment mechanisms, it effectively identifies the risks of frequent constraint state switching and cross-arm saturation transmission in dual-arm collaboration; through a risk intensity-driven weight scheduling and recovery mechanism, the constraint penalty weight, constraint bandwidth, and power weight can be dynamically adjusted according to the operating state, ensuring safety while avoiding excessive conservatism; the linearized coefficient matrix and local sensitivity matrix constructed based on input increment and state increment provide a computationally feasible approximate modeling method for nonlinear constraints; the future cost approximation model is updated online using the recursive least squares method, ensuring that the future cost estimate can be continuously corrected according to changes in operating conditions, thereby maintaining the foresight and accuracy of the control strategy. Overall, this invention, by constructing an objective function that superimposes stage costs and future costs and solving a convex quadratic programming problem, generates optimal control inputs that satisfy constraints in real time, achieving high-precision and high-stability control of nonlinear MIMO systems in actual operation. Attached Figure Description

[0015] Figure 1 This is a flowchart of the online optimal control method for nonlinear MIMO systems based on approximate dynamic programming, as described in this invention. Figure 2 This is a flowchart of the process for generating the constraint alternation determination result of the present invention; Figure 3 This is a flowchart of the process for generating the anti-saturation determination result according to the present invention. Detailed Implementation

[0016] The subject matter described herein will now be discussed with reference to exemplary embodiments. It should be understood that these embodiments are discussed only to enable those skilled in the art to better understand and implement the subject matter described herein, and changes may be made to the function and arrangement of the elements discussed without departing from the scope of this specification. Various processes or components may be omitted, substituted, or added as needed in the examples. Furthermore, features described in some examples may be combined in other examples.

[0017] It should be noted that, unless otherwise defined, the technical or scientific terms used in one or more embodiments of the present invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in one or more embodiments of the present invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" indicate that the element or object preceding the term encompasses the elements or objects listed following the term and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0018] like Figures 1-3 As shown, the online optimal control method for a nonlinear MIMO system based on approximate dynamic programming includes the following steps: Step 1: Obtain the system status, control commands, alignment error, constraint values, power data, and mobile chassis attitude of the dual flexible robotic arms, and perform preprocessing to obtain preprocessed data; Step 2: Generate constraint active indication sequences for the two robotic arms based on the preprocessed data, calculate the alternation switching rate and phase lag, and obtain the constraint alternation determination result; Step 3: Obtain the saturation shear amount based on the difference between the unlimited control command and the actual control command in the previous control cycle. Combined with the power increment of the other robotic arm, calculate the cross-arm load transfer coefficient and lag time to obtain the anti-saturation judgment result; schedule the constraint penalty weight, constraint bandwidth and power weight. Step 4: Within the sliding time window, obtain the local sensitivity matrix and linearization coefficient matrix based on the input increment and state increment of the preprocessed data; Step 5: Based on the scheduling constraint penalty weight, constraint bandwidth and power weight, as well as the local sensitivity matrix and linearization coefficient matrix, construct the stage cost function and superimpose the future cost approximation model. Under the premise of satisfying the input limiting and linearization constraint model, solve to obtain the unlimited control command and generate the actual control command. Step 6: Based on the current stage cost, the residual between the discounted approximate future cost of the next control period and the approximate future cost of the current control period, recursive least squares is used to update the parameters of the future cost approximation model online. Step 7: When the constraint alternation judgment result and the anti-saturation judgment result are lower than the corresponding threshold, restore the constraint penalty weight, constraint bandwidth and power weight to the baseline value.

[0019] In one embodiment of the present invention, system state refers to physical quantities reflecting the dynamic and static characteristics of the robotic arm, such as the joint angles, angular velocities, end-effector pose, and vibration mode information of the flexible arm itself; control commands refer to motion signals generated by manual input, including joint positions, end-effector trajectories, and torque or speed commands required for the operation task; alignment error refers to the deviation between the actual position of the robotic arm end-effector and the reference target position, used to characterize the deviation of the robotic arm from the target task's accuracy; constraint values ​​refer to the real-time values ​​of constraint functions involved in the system's operation, including joint angle limits, speed and acceleration limits, collision avoidance constraints, and the allowable power or torque range of the actuator, used to characterize whether the system is in a constraint-triggered state; power data refers to the power input and output values ​​of the robotic arm and its driver, including motor current, voltage, and corresponding mechanical power information, used to reflect the energy consumption and load transfer of the actuator; and mobile chassis attitude refers to the attitude information of the mobile platform mounted on the dual flexible robotic arms in the workspace, including position coordinates and attitude angles.

[0020] The preprocessing includes steps such as timestamp synchronization, filtering and denoising, coordinate system alignment, and data normalization. Timestamp synchronization ensures that data from multiple sensors are aligned within the same sampling period; filtering and denoising eliminates sensor jitter and environmental interference; coordinate system alignment unifies the reference frame between the robotic arm and the chassis; and normalization eliminates differences in the dimensions and orders of magnitude of different physical quantities. In this invention, the sampling period is consistent with the control period. Through these preprocessing steps, preprocessed data that is time-consistent, scale-uniform, and noise-suppressed can be obtained.

[0021] In one embodiment of the present invention, a constraint activity indication sequence for two robotic arms is generated based on preprocessed data, the alternation switching rate and phase lag are calculated, and a constraint alternation determination result is obtained, including: Step 11: For the two robotic arms, at the current sampling time, take the maximum value of each constraint value and compare it with the preset constraint judgment threshold. If the maximum value is greater than or equal to the preset constraint judgment threshold, generate a constraint activity indicator value of 1 for the corresponding robotic arm; otherwise, generate a constraint activity indicator value of 0. Step 12: Within the sliding time window, the constraint activity indicator values ​​of the continuously sampled robotic arm are arranged sequentially to form the constraint activity indicator sequence of the robotic arm. Step 13: Within the sliding time window, count the number of state transitions of the constraint active indication sequence of the two robotic arms respectively. The number of state transitions refers to the total number of changes from 0 to 1 and from 1 to 0. The ratio of the number of state transitions to the length of the sliding time window is used as the alternation rate. Step 14: Within the sliding time window, calculate the normalized cross-correlation function of the constraint active indication sequences of the two robotic arms, and determine the time shift corresponding to the maximum value of the cross-correlation function as the phase lag. When there are multiple maximum values, take the smallest positive time shift as the unique phase lag. The normalized cross-correlation function is used to measure the similarity between the constraint active indication sequences of the two robotic arms and their relative temporal relationship. Step 15: When the alternation switching rate is greater than the preset switching rate threshold and the phase lag is within the preset lag interval, the constraint alternation judgment result is 1; otherwise, it is 0, which respectively indicate that alternation exists and alternation does not exist.

[0022] Through the above processing, this embodiment realizes the automatic identification of the alternating constraint phenomenon of the two flexible robotic arms. The alternation switching rate quantifies the frequency of the active constraint state and can reflect the trend of system instability. The phase lag quantifies the time difference characteristics of the constraint triggering of the two robotic arms and can identify whether there is an alternation phenomenon. The combination of the two can significantly improve the robustness and reliability of the results, thereby providing accurate input basis for subsequent control weight scheduling.

[0023] In one embodiment of the present invention, the process of generating the anti-saturation determination result includes: Step 21: At the current sampling time, calculate the difference vector between the unlimited control command and the actual control command of the tested robotic arm, and take the Euclidean norm of the difference vector as the saturation shear quantity; calculate the power value difference of the other robotic arm at adjacent sampling times as the power increment, and perform unified timestamp alignment within the sliding time window to avoid timing deviation; wherein, the saturation shear quantity is used to quantify the degree of reduction caused by the limited amplitude constraint of the control command. Step 22: Within the sliding time window, calculate the normalized cross-correlation function of the sequences corresponding to the saturation shear rate and the power increment, and determine the minimum positive time shift corresponding to the maximum value of the normalized cross-correlation function as the lag time. When there are multiple maximum values, take the minimum positive time shift as the unique lag time. Step 23: Under the lag time, the least squares regression method is used to calculate the crossarm load transfer coefficient between the saturation shear and the power increment, and a non-negative constraint is applied to the crossarm load transfer coefficient to ensure that only physically reasonable positive values ​​are retained. Step 24: When the saturation shear is greater than the preset amplitude threshold, the crossarm load transfer coefficient is greater than the preset coefficient threshold, and the lag time is within the preset time shift interval, the anti-saturation judgment result is output as having crossarm saturation transfer risk; otherwise, the output is no risk.

[0024] Through the above steps, this embodiment can accurately identify cross-arm load transfer risks during dual-arm collaborative operation. The saturation shear rate reflects the degree of constraint of the tested robotic arm, the power increment and lag time reveal the dynamic coupling relationship between the other robotic arm and the tested robotic arm, and the cross-arm load transfer coefficient provides a quantitative indicator. The joint determination of these three factors significantly improves the accuracy and robustness of risk identification, providing a reliable basis for subsequent control weight scheduling.

[0025] In one embodiment of the present invention, in order to dynamically balance constraint satisfaction and performance indicators in the online control of a dual flexible robotic arm, a scheduling mapping is performed on constraint penalty weights, constraint bandwidth, and power weights, including: Step 31: When the constraint alternation judgment result is 1, the alternation switching rate is normalized according to the preset switching rate threshold to obtain the alternation risk intensity; when the anti-saturation judgment result is that there is a risk of cross-arm saturation transmission, the cross-arm load transmission coefficient is normalized according to the preset coefficient threshold to obtain the saturation risk intensity; wherein, the saturation risk intensity is used to measure the degree to which the tested robotic arm may cause load transmission risk to another robotic arm due to saturation shear. Step 32: Based on alternating risk intensity and saturation risk intensity, adjust the constraint penalty weight, constraint bandwidth, and power weight using monotonic linear functions, and limit them within preset upper and lower bounds; specifically, the constraint penalty weight increases monotonically with the increase of risk intensity; the constraint bandwidth decreases monotonically with the increase of alternating risk intensity; and the power weight increases monotonically with the increase of saturation risk intensity. Step 33: Update the adjusted constraint penalty weights, constraint bandwidth, and power weights using an exponential weighting method to obtain the scheduled constraint penalty weights, constraint bandwidth, and power weights. Specifically, the update formula is: wt = β⋅wt' + (1−β)⋅wt−1, where wt represents the updated weight parameters, including constraint penalty weights, constraint bandwidth, and power weights; wt' represents the weight parameters adjusted in the current control cycle in step 32; wt−1 represents the updated weight parameters from the previous control cycle; and β represents the smoothing coefficient, with a value ranging from 0 to 1. Through the above steps, this embodiment realizes the dynamic scheduling of constraint penalty weights, constraint bandwidth, and power weights; the alternating risk intensity and saturation risk intensity provide risk quantification indicators, the monotonic linear mapping ensures the interpretability and monotonicity of weight adjustment, and the exponential weighted update ensures the stability of the scheduling process.

[0026] In one embodiment of the present invention, within a sliding time window, a local sensitivity matrix and a linearization coefficient matrix are obtained based on the input increment and state increment of the preprocessed data, including: Step 41: Within the sliding time window, the difference between the current sampled input and the previous sampled input is used as the input increment, and the difference between the current sampled state and the previous sampled state is used as the state increment. The input increment sequence and the state increment sequence are then aligned one by one according to a unified timestamp. Specifically, the current sampled input represents the control command vector, and the current sampled state represents the system state vector, including joint angles, angular velocities, and related attitude variables. Step 42: Based on the concatenation of the input increment sequence and the state increment sequence to form a regression matrix, the least squares method is used to perform regression estimation to obtain a linearized coefficient matrix, which is used to reflect the approximate linear mapping relationship between the input increment and the state increment. Step 43: Limit the numerical amplitude of the linearized coefficient matrix, and trigger consistency check and truncation when the matrix norm exceeds the threshold, so as to avoid numerical divergence caused by data anomalies or noise and ensure the numerical stability and physical rationality of the results. Step 44: The mapping coefficient matrix from input increment to state increment in the linearization coefficient matrix is ​​determined as the local sensitivity matrix, and updated using an exponential weighting method to obtain the local sensitivity matrix for the current control cycle. The local sensitivity matrix is ​​used to quantify the system state's response strength to input disturbances.

[0027] Through the above processing, this embodiment ensures the consistency of input increment and state increment in the time domain and the effectiveness of linear fitting; through numerical amplitude limitation and exponential weighted update mechanism, the stability and robustness of local sensitivity matrix and linearization coefficient matrix in the iteration process are enhanced, providing a reliable basis for subsequent constraint model construction and control command optimization.

[0028] In one embodiment of the present invention, the process of constructing a stage cost function and superimposing a future cost approximation model includes: Step 51: Obtain the state deviation, input increment, and constraint relaxation amount of the current control cycle. The state deviation is the difference between the current system state and the preset reference trajectory, represented by a vector. The input increment is the difference between the current input and the input of the previous control cycle, also represented by a vector. The constraint relaxation amount is the excess amount of the linearized constraint exceeding the preset bandwidth, used to quantify the compensation variable when the constraint is slightly violated. Under multiple constraints, the constraint relaxation amount is represented by a vector. Step 52: Construct a stage cost function based on the state deviation term, input increment term, and constraint relaxation term. The state deviation term is obtained by weighting the state deviation with a preset weight matrix. The input increment term is obtained by weighting the power weight with the input penalty weight. The constraint relaxation term is obtained by multiplying the square of the constraint relaxation amount with the constraint penalty weight. Step 53: The approximate future cost model determined in the previous control cycle is superimposed onto the stage cost function according to a preset discount factor to form the overall objective function for this control cycle. The approximate future cost model is a parameterized function, and its parameters are updated online using a recursive least squares method. Specifically, the formula for calculating the overall objective function is: Jtotal(k)=ex(k)TQex(k)+Δu(k)T(γpI+γuI)Δu(k)+γc‖ξ(k)‖22+λ⋅Jfuture(x(k+1)); Where Jtotal(k) represents the value of the total objective function, k represents the control cycle index, ex(k) represents the state deviation, x represents the system state, T represents the transpose operation, Q represents the state weight matrix, Δu(k) represents the input increment, γp and γu represent the power weight and input penalty weight, respectively, reflecting the degree of penalty for input energy consumption and the degree of suppression of input mutations. The combination of the two can ensure that the input is both economical and stable. γc represents the constraint penalty weight, which is used to adjust the strength of the constraint relaxation in the cost function. The larger the value, the lower the system's tolerance for constraint violation. I represents the identity matrix, ξ(k) represents the constraint relaxation, λ represents the discount factor, and Jfuture(x(k+1)) represents the future cost approximation model. It takes the predicted state of the next control cycle as input and outputs the future cost approximation. Preferably, the future cost approximation model is a quadratic function with the state deviation as the independent variable: Jfuture(x(k+1))=ex(k+1)TP ex(k+1)+qTex(k+1)+c, where x(k+1) represents the system state in control period k+1, represented by a vector; ex(k+1) represents the state deviation in the next control period; P represents a symmetric positive semi-definite matrix used to quantify the importance of different state components in the future cost approximation; q represents the weight vector of the first-order term; and c represents the constant term. Preferably, q is a zero vector and c is a non-negative constant. The above parameters are updated online using a recursive least squares method, combined with a forgetting factor to adapt to the dynamic changes in operating conditions.

[0029] In this embodiment, the stage cost function ensures that the state control and constraints of the current control cycle are satisfied, and the future cost approximation model provides a prediction of the performance of subsequent control cycles. The combination of the two can realize the idea of ​​approximate dynamic programming in the finite time domain. This design not only improves the foresight of the control strategy, but also enhances the robustness and stability of the method in complex nonlinear systems.

[0030] In one embodiment of the present invention, the generation of the unlimited control command and the actual control command includes: Step 61: Based on the overall objective function, apply upper and lower bound constraints on the input, and use the linearization coefficient matrix and constraint bandwidth to perform a first-order expansion of the nonlinear constraint function at the current state and input point to form a linearized constraint model. The nonlinear constraint function is a function composed of the system state and control input, and has non-convexity and computational complexity. By performing a first-order expansion at the current state and input point, and adding constraint bandwidth as a tightening margin on the right side, the complex nonlinear constraints can be effectively transformed into a linear constraint form that can be handled by the convex optimization problem, while maintaining a certain degree of conservatism, thereby enhancing the stability and feasibility of the solution. Step 62: Based on the overall objective function, construct and solve the convex quadratic programming problem to obtain the optimal input increment for the current control cycle. The convex quadratic programming problem refers to an optimization problem with a quadratic objective function and linear constraints, which has a globally unique solution and can guarantee the efficiency and stability of real-time calculation. Step 63: Add the optimal input increment to the input of the previous control cycle to obtain the unlimited control command. This unlimited control command reflects the desired control action of the system under the condition of not considering the input saturation limit. Step 64 compares the unlimited control command with the upper and lower limits of the input. If it exceeds the upper limit, the upper limit is used; if it is below the lower limit, the lower limit is used, thus obtaining the actual control command. This step ensures that the control command meets the physical constraints of the actuator, avoiding system instability or hardware damage due to control input exceeding limits.

[0031] The processing procedure in this embodiment realizes a complete closed loop from optimizing the objective function to generating executable control instructions. The linearized constraint model ensures the computability and safety margin of complex constraints; the quadratic programming solution ensures that the control input achieves the optimal balance between performance and constraints; and the final limiting operation ensures the implementability of the instructions and hardware security.

[0032] In one embodiment of the present invention, the parameters of the future cost approximation model are updated online during the control operation using a recursive least squares method. At the end of each control cycle, a consistency error signal is constructed based on the difference between the sum of the stage cost and the discounted future cost approximation of the next control cycle and the future cost approximation of the current control cycle. This error signal is then used as a driving input to recursively correct the parameter vector of the future cost approximation model, while simultaneously recursively updating the covariance matrix.

[0033] Specifically, the consistency error signal reflects the deviation between the future cost approximation model and the actual cost recursion relationship. The future cost approximation model adopts a linear parameterized form with a fixed eigenmap, and its input is an eigenvector composed of system state components and their quadratic terms, thus simultaneously capturing the linear and nonlinear effects of the state on the future cost approximation. Using this eigenvector as a regressor, combined with the consistency error signal, the parameter vector is recursively corrected using the recursive least squares method, and the corresponding covariance matrix is ​​updated synchronously. The covariance matrix is ​​initialized as a diagonal positive definite matrix and maintains positive definiteness during the update process to ensure the numerical stability of the parameter estimation process.

[0034] Through the above process, the future cost approximation model can continuously and adaptively update during system operation, thereby gradually approximating the theoretical future cost function. Among them, the introduction of the consistency error signal ensures that the update direction is consistent with the real dynamic programming equation; the recursive least squares algorithm provides efficient online estimation capability, while the dynamic update of the covariance matrix ensures the stability of parameter convergence.

[0035] In one embodiment of the present invention, when the constraint alternation determination result and the anti-saturation determination result are lower than the corresponding thresholds, the constraint penalty weight, constraint bandwidth, and power weight are restored to the baseline values, including: Step 71: Detect alternating risk intensity and saturation risk intensity within the sliding time window. Trigger recovery operation when both are continuously lower than their respective preset thresholds for a preset number of times. Step 72: Call the preset baseline values ​​and upper and lower bound parameters of the system calibration stage, and read the baseline values ​​and corresponding safety boundaries of constraint penalty weight, constraint bandwidth and power weight respectively; Step 73: When recovery is triggered, the current constraint penalty weight, constraint bandwidth and power weight are linearly weighted and averaged according to the preset regression rate to return to the corresponding benchmark value. After each update, the upper and lower bounds are pruned to ensure that the updated values ​​are always within the allowable range and to avoid affecting control performance or causing instability due to numerical overshoot.

[0036] Through the above process, this embodiment can proactively remove the conservatism brought about by weight scheduling when the system risk is significantly reduced, so that the constraint penalty weight, constraint bandwidth and power weight are restored to the calibration state. It not only ensures the adaptive adjustment of control parameters under high risk conditions, but also restores them to the baseline settings in a timely manner after the risk subsides, thereby achieving a dynamic balance between safety and performance and improving the flexibility and robustness of the overall control method.

[0037] It should be noted that the interval and threshold sizes are set for ease of comparison. The size of the threshold depends on the amount of sample data and the base number set by those skilled in the art for each set of sample data, as long as it does not affect the proportional relationship between the parameter and the quantized value. Furthermore, the above formulas are all dimensionless calculations, and the formulas are derived from software simulations using a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0038] The embodiments of the present invention have been described above, but the embodiments are not limited to the specific implementation methods described above. The specific implementation methods described above are merely illustrative and not restrictive. Those skilled in the art can make many other forms under the guidance of the embodiments described above, all of which are within the protection scope of the embodiments described above.

Claims

1. An online optimal control method for nonlinear MIMO systems based on approximate dynamic programming, characterized in that, Includes the following steps: Step 1: Obtain the system status, control commands, alignment error, constraint values, power data, and mobile chassis attitude of the dual flexible robotic arms, and perform preprocessing to obtain preprocessed data; Step 2: Generate constraint active indication sequences for the two robotic arms based on the preprocessed data, calculate the alternation switching rate and phase lag, and obtain the constraint alternation determination result; Step 3: Obtain the saturation shear amount based on the difference between the unlimited control command and the actual control command in the previous control cycle. Combined with the power increment of the other robotic arm, calculate the cross-arm load transfer coefficient and lag time to obtain the anti-saturation judgment result; schedule the constraint penalty weight, constraint bandwidth and power weight. Step 4: Based on the input increment and state increment of the preprocessed data, obtain the local sensitivity matrix and linearization coefficient matrix; Step 5: Based on the scheduling constraint penalty weight, constraint bandwidth and power weight, as well as the local sensitivity matrix and linearization coefficient matrix, construct the stage cost function and superimpose the future cost approximation model. Under the premise of satisfying the input limiting and linearization constraint model, solve to obtain the unlimited control command and generate the actual control command. Step 6: Based on the current stage cost, the residual between the discounted approximate future cost of the next control period and the approximate future cost of the current control period, update the parameters of the future cost approximation model using recursive least squares. Step 7: When the constraint alternation judgment result and the anti-saturation judgment result are lower than the corresponding threshold, restore the constraint penalty weight, constraint bandwidth and power weight to the baseline value.

2. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, The generation of the constraint alternation determination result includes: Step 11: For the two robotic arms, at the current sampling time, take the maximum value of each constraint value and compare it with the preset constraint judgment threshold. If the maximum value is greater than or equal to the preset constraint judgment threshold, generate a constraint activity indicator value of 1 for the corresponding robotic arm; otherwise, generate a constraint activity indicator value of 0. Step 12: Within the sliding time window, the constraint activity indicator values ​​of the continuously sampled robotic arm are arranged sequentially to form the constraint activity indicator sequence of the robotic arm. Step 13: Within the sliding time window, count the number of state transitions of the constraint active indication sequence of the two robotic arms, and use the ratio of the number of state transitions to the length of the sliding time window as the alternation rate. Step 14: Within the sliding time window, calculate the normalized cross-correlation function of the constraint active indication sequence of the two robotic arms, and determine the time shift corresponding to the maximum value of the cross-correlation function as the phase lag. Step 15: When the alternation switching rate is greater than the preset switching rate threshold and the phase lag is within the preset lag interval, the constraint alternation judgment result is 1; otherwise, it is 0, which respectively indicate that alternation exists and alternation does not exist.

3. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, The process of generating the anti-saturation determination result includes: Step 21: At the current sampling time, calculate the difference vector between the unlimited control command and the actual control command of the tested robotic arm, and take the Euclidean norm of the difference vector as the saturation shear quantity; calculate the power value difference of the other robotic arm at adjacent sampling times as the power increment, and perform unified timestamp alignment within the sliding time window. Step 22: Within the sliding time window, calculate the normalized cross-correlation function of the sequences corresponding to the saturation shear rate and the power increment, and determine the minimum positive time shift corresponding to the maximum value of the normalized cross-correlation function as the lag time. Step 23: Under the lag time, the least squares regression method is used to calculate the crossarm load transfer coefficient between the saturation shear and the power increment; Step 24: When the saturation shear is greater than the preset amplitude threshold, the crossarm load transfer coefficient is greater than the preset coefficient threshold, and the lag time is within the preset time shift interval, the anti-saturation judgment result is output as having crossarm saturation transfer risk; otherwise, the output is no risk.

4. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, The scheduling mapping for constraint penalty weights, constraint bandwidth, and power weights includes: Step 31: When the constraint alternation judgment result is 1, the alternation switching rate is normalized according to the preset switching rate threshold to obtain the alternation risk intensity; when the anti-saturation judgment result is that there is a risk of cross-arm saturation transmission, the cross-arm load transmission coefficient is normalized according to the preset coefficient threshold to obtain the saturation risk intensity. Step 32: Based on alternating risk intensity and saturation risk intensity, use monotonic linear functions to adjust the constraint penalty weight, constraint bandwidth and power weight respectively, and limit them within the preset upper and lower bounds. Step 33: Update the adjusted constraint penalty weights, constraint bandwidth, and power weights using an exponential weighting method to obtain the scheduled constraint penalty weights, constraint bandwidth, and power weights.

5. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, Based on the input increment and state increment of the preprocessed data, the local sensitivity matrix and linearization coefficient matrix are obtained, including: Step 41: Within the sliding time window, the difference between the current sampled input and the previous sampled input is used as the input increment, and the difference between the current sampled state and the previous sampled state is used as the state increment. The input increment sequence and the state increment sequence are then aligned one by one according to the unified timestamp. Step 42: Based on the concatenation of the input increment sequence and the state increment sequence to form a regression matrix, the least squares method is used to perform regression estimation to obtain the linearized coefficient matrix; Step 43: Limit the numerical magnitude of the linearized coefficient matrix, and trigger consistency verification and truncation when the matrix norm exceeds the threshold; Step 44: Determine the mapping coefficient matrix from input increment to state increment in the linearization coefficient matrix as the local sensitivity matrix, and update it using an exponential weighting method to obtain the local sensitivity matrix of the current control cycle.

6. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, The process of constructing the stage cost function and superimposing the future cost approximation model includes: Step 51: Obtain the state deviation, input increment, and constraint relaxation amount of the current control cycle. The state deviation is the difference between the current system state and the preset reference trajectory, the input increment is the difference between the current input and the input of the previous control cycle, and the constraint relaxation amount is the excess amount of the linearized constraint exceeding the preset bandwidth. Step 52: Construct a stage cost function based on the state deviation term, input increment term, and constraint relaxation term. The state deviation term is obtained by weighting the state deviation with a preset weight matrix. The input increment term is obtained by weighting the power weight with the input penalty weight. The constraint relaxation term is obtained by multiplying the square of the constraint relaxation amount with the constraint penalty weight. Step 53: The approximate future cost model determined in the previous control cycle is superimposed onto the stage cost function according to a preset discount factor to form the overall objective function of the current control cycle. The approximate future cost model is a parameterized function, and its parameters are updated online by the recursive least squares method.

7. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 6, characterized in that, The generation of unlimited control commands and actual control commands includes: Step 61: Based on the overall objective function, apply upper and lower limit constraints to the input, and use the linearization coefficient matrix and constraint bandwidth to perform a first-order expansion of the nonlinear constraint function at the current state and input point to form a linearized constraint model. Step 62: Based on the overall objective function, construct and solve the convex quadratic programming problem to obtain the optimal input increment for the current control cycle; Step 63: Add the optimal input increment to the input of the previous control cycle to obtain the unlimited control command; Step 64: Compare the unlimited control command with the input upper and lower limits. If it exceeds the upper limit, take the upper limit; if it is below the lower limit, take the lower limit to obtain the actual control command.

8. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, The parameters of the future cost approximation model are updated online during the control operation using a recursive least squares method. At the end of each control cycle, a consistency error signal is constructed based on the difference between the sum of the stage cost and the discounted future cost approximation of the next control cycle and the future cost approximation of the current control cycle. This error signal is then used as the driving input to recursively correct the parameter vector of the future cost approximation model, while simultaneously recursively updating the covariance matrix.

9. The online optimal control method for nonlinear MIMO systems based on approximate dynamic programming according to claim 1, characterized in that, When the constraint alternation judgment result and the anti-saturation judgment result are lower than the corresponding threshold, the constraint penalty weight, constraint bandwidth, and power weight will be restored to the baseline value, including: Step 71: Detect alternating risk intensity and saturation risk intensity within the sliding time window. Trigger recovery operation when both are continuously lower than their respective preset thresholds for a preset number of times. Step 72: Call the preset baseline values ​​and upper and lower bound parameters of the system calibration stage, and read the baseline values ​​and corresponding safety boundaries of constraint penalty weight, constraint bandwidth and power weight respectively; Step 73: When recovery is triggered, perform a linear weighted average of the current constraint penalty weight, constraint bandwidth and power weight according to the preset regression rate, so that it regresses to the corresponding benchmark value.

Images