A multi-actuator biaxial attitude constraint predictive control method

Abstract

The present application belongs to the technical field of automatic attitude control, and discloses a multi-actuator two-axis attitude constraint predictive control method, which is used for slow-changing attitude extraction and control opportunity determination under variable-speed rotation disturbance. The two-axis slow-changing attitude state is obtained by resampling the two-axis inclination signal based on equal-phase pulses, synchronous averaging, zero-mean disturbance deduction, double-sided truncated averaging and slow feature analysis; the state is projected to multiple constraint directions, the prediction first arrival time is calculated according to the drift coefficient and diffusion coefficient of the projection increment; the candidate execution amount is generated under the constraints of actuator stroke, rate and common modal, the post-execution attitude is predicted by using the online control performance matrix, and the target execution amount is determined according to the attitude constraint margin, execution energy consumption and execution amount change. When the minimum prediction first arrival time is not greater than the sum of the execution time delay and the safety margin, the target execution amount is executed, and the control performance matrix is updated after the response is stable. The method is suitable for multi-actuator attitude stabilization control systems containing rotation disturbance sources.

Description

Technical Field

[0001] This invention belongs to the field of automatic control technology, specifically relating to a multi-actuator dual-axis attitude constraint predictive control method. Background Technology

[0002] Under the influence of load transfer, changes in support stiffness, slow foundation settlement, and periodic excitation, the tilt angles of the first and second attitude axes of a platform with rotating working components will slowly shift. Existing leveling control typically uses a threshold judgment for the current tilt angle and drives each support actuator according to a fixed ratio, without taking into account the time required for the state to reach the constraint boundary, the time required for the actuator to complete adjustment, and the cross-coupling between different actuators. The tilt angle signal also contains periodic disturbances locked to the rotation phase and random vibrations unrelated to the rotation phase; when the rotation speed changes, the time-domain frequency of the periodic disturbance changes accordingly and may overlap with the frequency band of random vibrations, making it difficult to stably retain the slow attitude shift using fixed cutoff frequency filtering. In addition, changes in load position, support contact state, and hydraulic parameters will change the dual-axis attitude response corresponding to a unit actuation quantity, and a fixed control performance relationship is prone to under-adjustment, over-adjustment, or mutual coupling between the two attitude axes. Therefore, it is necessary to establish an attitude state extraction method that does not rely on a constant rotation speed, determine the control timing based on the attitude constraint arrival probability and execution delay, generate control quantities under multi-actuator constraints, and use the actual response to correct the control performance relationship online. Summary of the Invention

[0003] The purpose of this invention is to provide a multi-actuator dual-axis attitude constraint predictive control method to solve the technical problems of slow-changing attitude state being difficult to extract stably when rotational synchronous disturbances and random vibrations are superimposed, current value triggering not taking into account execution delay, multi-directional attitude constraints not being handled uniformly, and fixed control performance relationships not being able to adapt to changes in operating state.

[0004] To address the aforementioned technical problems, this invention provides a multi-actuator dual-axis attitude constraint predictive control method, applied to an attitude control system including a dual-axis tilt sensor, a rotation phase sensor, an attitude controller, and at least three attitude adjustment actuators, comprising the following steps: Step 1: Using the isophase pulses output by the rotating phase sensor as a reference, the first attitude axis tilt angle signal and the second attitude axis tilt angle signal are resampled in the same phase. The zero-mean rotational synchronization disturbance is extracted according to the rotation period. The corrected sampled values ​​after deducting the disturbance are averaged by two-sided truncation to form a dual-axis imbalance characterization sequence. Step 2: Perform time-delay embedding, regularized whitening, and slow feature extraction on the biaxial imbalance representation sequence to obtain the biaxial slow-change attitude state; Step 3: Project the dual-axis slowly changing attitude state onto multiple attitude constraint directions, determine the drift coefficient and diffusion coefficient based on the increment of the projection sequence in each direction, and calculate the predicted first arrival time to the corresponding attitude constraint boundary. Step 4: Generate candidate execution quantities according to the stroke constraints, rate constraints and common modal constraints that do not cause the overall lifting and lowering of the controlled platform of the attitude adjustment actuator. Use the control efficiency matrix to predict the dual-axis attitude state after each candidate execution quantity is completed, and determine the target execution quantity based on the predicted attitude constraint margin, execution energy consumption and execution quantity change. Step 5: When the minimum predicted first arrival time is not greater than the sum of the execution delay of the target execution amount and the preset safety margin, or when the dual-axis slow-change attitude state reaches any attitude constraint boundary, execute the target execution amount; after the response stabilizes, update the control performance matrix according to the actual dual-axis attitude change amount, and return to step 1.

[0005] Furthermore, the attitude adjustment actuator is a first to a fourth electro-hydraulic support actuator installed at four support positions on the controlled bearing platform, wherein at least three of the four support positions are not collinear; each electro-hydraulic support actuator includes a servo valve, a hydraulic cylinder, a displacement sensor, and a pressure sensor, and the actuation amount is the hydraulic cylinder extension increment; the common modal constraint is that the sum of the extension increments of the four hydraulic cylinders is equal to 0.

[0006] Further, step 1 includes: recording the timestamps of each equal-phase pulse, performing linear interpolation on the dual-axis tilt signal according to the timestamps; dividing the data segment into 1024 sampling points per revolution of the rotating component, averaging the most recent 32 data segments at the same phase position to obtain the rotation synchronization component, subtracting the mean value of the rotation synchronization component point by point to obtain the zero-mean rotation synchronization disturbance; subtracting the zero-mean rotation synchronization disturbance point by point from the latest data segment, performing a two-sided truncated average on the obtained corrected sampling value to obtain the imbalance characterization value of the corresponding attitude axis; and performing sliding median filtering on the imbalance characterization values ​​of the current cycle and the previous 4 cycles.

[0007] Further, in step 2, taking each rotation cycle as one time step, the imbalance representation values ​​of the first and second attitude axes in the most recent 120 time steps are embedded with a 5th-order time delay to form a 6-dimensional sample vector; the sample vector is mean-removed and regularized whitened, the covariance matrix of the first-order difference between adjacent whitened sample vectors is calculated, the eigenvector corresponding to the smallest eigenvalue is taken as the slow direction vector, and the whitened sample vector is projected onto the slow direction vector to obtain the slow feature sequence of the corresponding attitude axis.

[0008] Furthermore, within the same sliding window, a least-squares linear fit is performed with the slow feature sequence as the independent variable and the imbalance representation value of the corresponding time step as the dependent variable. The slope and intercept obtained from the fit are used to convert the slow feature sequence into a slow-varying estimation sequence with the same dimensions as the tilt angle signal and the same sign as the tilt direction. The slow-varying estimation sequences of the two attitude axes constitute the dual-axis slow-varying attitude state.

[0009] Furthermore, multiple attitude constraint directions include the positive and negative directions of the first and second attitude axes, as well as four oblique directions formed by the normalized and equally weighted combination of the two attitude axes. For each attitude constraint direction, the dual-axis slowly changing attitude state is projected onto that direction, and the arithmetic mean and sample standard deviation of the increments of adjacent time steps in the projection sequence are taken as the original drift coefficient and diffusion coefficient, respectively. When the original drift coefficient is greater than 0 and less than a preset drift lower limit, the preset drift lower limit is used as the drift coefficient for prediction; otherwise, the original drift coefficient is used as the drift coefficient for prediction. When the diffusion coefficient is less than the preset lower diffusion limit and the predicted drift coefficient is greater than 0, the predicted first arrival time is determined by dividing the boundary distance by the predicted drift coefficient; when the diffusion coefficient is less than the preset lower diffusion limit and the predicted drift coefficient is not greater than 0, the predicted first arrival time is set to the preset maximum value; when the diffusion coefficient is not less than the preset lower diffusion limit and the predicted drift coefficient is greater than 0, the predicted first arrival time is determined according to the median of the inverse Gaussian distribution; when the diffusion coefficient is not less than the preset lower diffusion limit and the predicted drift coefficient is not greater than 0, the predicted first arrival time is determined according to the first arrival probability of the zero-drift Wiener process.

[0010] Furthermore, the control performance matrix is ​​a 2×2 matrix that maps the first attitude axis control mode and the second attitude axis control mode to the first attitude axis tilt angle change and the second attitude axis tilt angle change. During the calibration phase, the forward execution amount and the reverse execution amount of the first attitude axis control mode and the second attitude axis control mode are applied respectively, and the tilt angle change of the two axes before and after execution is recorded synchronously. The initial value of the control performance matrix is ​​determined by least squares fitting.

[0011] Further, step 4 includes: combining the first attitude axis control mode and the second attitude axis control mode according to the preset displacement step size, and converting them into four hydraulic cylinder extension increments that satisfy the common mode constraints; generating speed-limited displacement trajectories based on the extension increments of each hydraulic cylinder and the corresponding allowable speeds, determining the maximum value of the time required for each speed-limited displacement trajectory as the execution delay, extrapolating the unadjusted attitude within the execution delay according to the drift coefficient, and superimposing the attitude change amount obtained by the control efficiency matrix acting on the two attitude axis control modes to obtain the predicted attitude after execution; eliminating candidate execution quantities that cause any hydraulic cylinder to exceed the stroke limit, and determining the candidate execution quantity with the smallest control objective function as the target execution quantity, wherein the control objective function is composed of a weighted sum of the predicted attitude constraint margin penalty term, the execution energy consumption term, and the change term relative to the previous execution quantity.

[0012] Furthermore, half of the first attitude axis control mode is assigned with the same sign to the first and second electro-hydraulic support actuators, and the other half is assigned with the opposite sign to the third and fourth electro-hydraulic support actuators; half of the second attitude axis control mode is assigned with the same sign to the first and third electro-hydraulic support actuators, and the other half is assigned with the opposite sign to the second and fourth electro-hydraulic support actuators; the two control modes are algebraically added at each electro-hydraulic support actuator to obtain the hydraulic cylinder extension / retraction increment from the first to the fourth electro-hydraulic support actuators.

[0013] Furthermore, after the target execution is completed and the tilt rate of change of the two attitude axes is continuously within the preset stable range, the actual control response is obtained by subtracting the natural drift during the execution period from the dual-axis attitude change before and after execution. The regression vector is composed of the first attitude axis control mode and the second attitude axis control mode. The difference between the actual control response and the predicted response of the control performance matrix is ​​used as the update residual. The control performance matrix is ​​corrected by recursive least squares with a forgetting factor. When the norm of the regression vector is less than the preset excitation lower limit, the amplitude of the actual control response is less than 3 times the standard deviation of the tilt noise, or the tilt rate of change has not entered the preset stable range, the correction of the control performance matrix is ​​paused.

[0014] The present invention has the following beneficial effects: This invention establishes a phase reference unaffected by rotational speed changes using equal-phase pulses. Rotational synchronization components are extracted through multi-cycle synchronous averaging, and the phase mean is removed point-by-point to obtain a zero-mean rotational synchronization disturbance. The corrected sampled values ​​are then subjected to a two-sided truncated average to deduct the actual components involved in the imbalance characterization calculation. Slow feature analysis is used to extract a biaxial slowly varying attitude state with angular dimensions and a defined direction. Projecting this state onto the axial and oblique attitude constraint directions allows for handling single-axis and biaxial combined offsets within the same probabilistic framework. The segmented use of a zero-drift Wiener process and an inverse Gaussian first-arrival time model maintains numerical stability when drift approaches zero and accounts for the risk of out-of-bounds movement due to random diffusion. The 2×2 control effectiveness matrix simultaneously characterizes the main action and cross-coupling of the two attitude control modes. Candidate control solutions are constrained by travel, rate, and common mode, allowing the determination of the target execution quantity without causing overall lifting or lowering of the controlled platform. The actual attitude response after execution is used to correct the control performance matrix online. Excitation, signal-to-noise ratio, and steady-state criteria are used to limit the impact of unsteady response on model updates, thereby improving the stability and adaptability of multi-actuator dual-axis attitude control. Attached Figure Description

[0015] Figure 1 A schematic diagram of the layout of a multi-actuator dual-axis attitude control system provided in an embodiment of the present invention; Figure 2A schematic diagram of the phase domain synchronous averaging principle provided in an embodiment of the present invention; Figure 3 A schematic diagram showing the comparison between the imbalance characterization sequence and the periodically average original dip angle provided in an embodiment of the present invention; Figure 4 This is a schematic diagram of the extrapolation and control triggering criteria for the Wiener process with drift provided in an embodiment of the present invention; Figure 5 This is a schematic diagram of the cumulative distribution function of the inverse Gaussian distribution at the first arrival time provided in an embodiment of the present invention; Figure 6 This is a schematic diagram of the eight-direction attitude constraint allowable domain provided in an embodiment of the present invention; Figure 7 This is a schematic diagram of the dual-axis control mode and the extension / retraction increment allocation of the four actuators provided in an embodiment of the present invention. Detailed Implementation

[0016] The method of the present invention will be further described in detail below with reference to the accompanying drawings and embodiments.

[0017] The multi-actuator dual-axis attitude constraint predictive control method of this embodiment is applied to a heavy-duty controlled bearing platform with a rotating working component. The attitude control system includes a dual-axis tilt sensor, a rotating phase sensor that outputs equal-phase pulses synchronously with the rotating working component, an attitude controller, and first to fourth electro-hydraulic support actuators, and operates cyclically according to the following steps: Step 1: The attitude controller resamples the tilt signals of the first and second attitude axes with equal-phase pulses as a reference. The data segment is divided according to the rotation of the rotating working part. The zero-mean rotational synchronization disturbance is extracted by multi-cycle synchronous averaging. The corrected sampled value after deducting the disturbance is truncated and averaged by both sides. The imbalance characterization value of the two attitude axes is obtained cycle by cycle. Step 2: Perform time delay embedding, regularization whitening and slow direction projection on the imbalance representation sequences of the two attitude axes respectively to obtain slow feature sequences, and then convert them into dual-axis slow attitude states with the same dimensions as the tilt angle signal and the same sign as the tilt direction. Step 3: Project the dual-axis slow-change attitude state onto the axial and oblique attitude constraint directions, estimate the drift coefficient and diffusion coefficient based on the incremental projection sequence of each direction, and calculate the predicted first arrival time to the corresponding attitude constraint boundary. Step 4: Based on the control performance matrix, the dual-axis slowly varying attitude state, the drift coefficient, and the execution delay of the candidate execution quantities, predict the dual-axis attitude state after each candidate execution quantity is completed, calculate the control objective function, and determine the target execution quantity under the constraints of travel, rate, and common modal. Step 5: When the minimum predicted first arrival time is not greater than the sum of the execution delay of the target execution amount and the preset safety margin, or when the dual-axis slow-change attitude state has reached any attitude constraint boundary, execute the target execution amount; after the response stabilizes, update the control performance matrix according to the actual dual-axis attitude change amount, and continue to execute the above loop.

[0018] The phase domain descrambling, first arrival time prediction, and multi-actuator control process are explained below with reference to specific embodiments.

[0019] The tilt angle changes of the controlled platform include slow attitude shifts, rotational synchronous disturbances, and asynchronous random vibrations. Slow attitude shifts are caused by load transfer, changes in support contact state, slow foundation settlement, and hydraulic cylinder leakage, with a timescale typically ranging from tens of seconds to minutes; these are the state variables that the attitude controller needs to correct. Rotational synchronous disturbances are caused by mass eccentricity of rotating components, assembly deviations, and periodic cutting loads; their waveforms repeat with each revolution of the rotating working component. Asynchronous random vibrations originate from transient loads, excitation from nearby equipment, and structural vibrations; in the example system, the main energy distribution is between 0.05 Hz and 0.5 Hz. When the rotational speed varies from 8 r / min to 25 r / min, the fundamental frequency of the rotational synchronous disturbance is between 0.13 Hz and 0.42 Hz, overlapping with the frequency band of the asynchronous random vibrations. Time-domain fixed cutoff frequency filtering cannot simultaneously suppress both types of disturbances and preserve slow attitude shifts across the entire speed range. After resampling the tilt angle signal according to the rotation phase, the rotational synchronization disturbance remains coherent at the same phase position in each cycle, while the asynchronous random vibration does not remain coherent at the same phase position. Therefore, it can be separated by multi-cycle synchronous averaging.

[0020] This embodiment employs four electro-hydraulic support actuators to adjust the tilt angles of the first and second attitude axes of the controlled platform. Each electro-hydraulic support actuator includes a proportional servo valve, a double-acting hydraulic cylinder, a magnetostrictive displacement sensor, and a pressure sensor. The attitude controller outputs the hydraulic cylinder extension / retraction increment, and the displacement closed-loop controller adjusts the servo valve opening to bring the hydraulic cylinder to the target position. The sum of the extension / retraction increments of the four hydraulic cylinders is zero to suppress the common lifting / lowering component of the platform.

[0021] Figure 1The arrangement of a multi-actuator dual-axis attitude control system is shown. The first to fourth electro-hydraulic support actuators are positioned at four support locations on the controlled platform, with at least three of these locations being non-collinear. A rotating working component is located in the center of the platform. A rotating phase sensor is coaxially mounted with the rotating working component. A dual-axis tilt sensor is installed near the geometric center of the platform. The attitude controller is electrically connected to the rotating phase sensor, the dual-axis tilt sensor, and the four electro-hydraulic support actuators. The rotating phase sensor outputs 1024 isophase pulses per revolution, with a phase increment of approximately 0.35° between adjacent pulses. The dual-axis tilt sensor samples at 200Hz, has a range of ±5°, and a resolution better than 0.001°. The attitude controller records the timestamps of the isophase pulses with a timing accuracy of 1μs. Let the timestamp of one isophase pulse be... , Indicates the first The arrival time of each equal-phase pulse; this time is located at the adjacent sampling time of the tilt angle signal. and Between, the linear interpolation result is ,in and Respectively represent and Adjacent before and after sampling times, and These represent the dual-axis tilt sensor in and The output tilt angle value, The interpolated first value represents the first value obtained from the interpolation. The tilt angle values ​​at each equal-phase pulse moment. In the example system, the highest frequency of the vibration component to be separated is 0.5Hz, the sampling interval is 5ms, and a linear approximation can be used between adjacent sampling points; after pulse-by-pulse interpolation, a phase domain tilt angle sequence with rotational phase as the independent variable is obtained.

[0022] During installation, the geometric center of the controlled platform is used as the origin of the coordinate system. The two measuring axes of the dual-axis tilt sensor are aligned with the first and second attitude axes, respectively. Under the reference attitude confirmed by a standard level, the platform is kept stationary and there are no external load changes. Tilt angle data is continuously collected for 60 seconds, and the arithmetic mean of the two axes is used as the initial zero bias. Subsequently, known small-angle tilts are applied to the first and second attitude axes sequentially. The scaling factor and the non-orthogonality error of the two axes are corrected based on the linear fitting result between the standard tilt angle and the sensor output. The calibrated scaling factor deviation is no greater than 0.2%, and the cross sensitivity of the two axes is no greater than 0.5%.

[0023] The displacement sensors of the first to fourth electro-hydraulic support actuators are coaxially mounted with the hydraulic cylinder piston rods; each of the two working chambers of each hydraulic cylinder is connected to a pressure sensor, forming a total of 8 pressure measurement points. The attitude controller updates the displacement closed-loop in 50 milliseconds, adjusting the opening of the proportional servo valve according to the difference between the target displacement and the actual displacement. Pressure data is used to determine the actuator's load status and calculate execution energy consumption, but does not directly replace displacement closed-loop feedback. When the pressure difference between the two chambers of any hydraulic cylinder exceeds 90% of the rated pressure difference, the attitude controller reduces the allowable speed of that actuator; when the pressure difference reaches the rated value or the displacement feedback is interrupted, the current attitude adjustment is stopped.

[0024] After system startup, a self-test is performed on sensor zero bias, displacement range, and communication status. The sampling timestamps of the dual-axis tilt sensor must monotonically increase, and the interval between adjacent samplings must not exceed twice the rated interval. The number of pulses output by the rotating phase sensor in each complete rotation cycle must be within the range of 99.5% to 100.5% of the nominal pulse count. If the self-test fails, probability prediction and candidate execution quantity solving are not initiated; the attitude controller maintains each hydraulic cylinder in its current position and outputs a fault status. Only after the self-test passes is a phase domain buffer established to avoid introducing invalid data from the startup phase into the rotational synchronization component.

[0025] The attitude controller activates its control functions in three phases. After passing the self-test and forming the first complete rotation cycle, it activates the current attitude boundary protection, but does not perform probabilistic extrapolation. After accumulating 32 effective rotation cycles, the phase domain synchronously averages for a fixed length and enters rated operation. After accumulating 120 effective time steps and obtaining a reliable dual-axis slowly varying attitude state, it activates multi-directional first-arrival time prediction and candidate execution quantity solving. When any direction reaches the attitude constraint boundary during the activation phase, the control performance matrix saved during the calibration phase is directly called to generate the execution quantity within the constraints. Boundary protection is still performed before the statistical window is formed, and the drift coefficient or diffusion coefficient with insufficient samples is not used for probabilistic prediction.

[0026] The phase domain tilt sequence is divided into data segments according to each revolution of the rotating working component. Each data segment contains 1024 sampling points. Each sampling point corresponds to a fixed phase position. This represents the phase position number, with values ​​ranging from 1 to 1024. The attitude controller maintains a circular buffer with a capacity of 32 data segments and calculates the arithmetic mean of the sampled values ​​of the most recent data segments at the same phase position within the buffer. ,in This indicates the number of data segments involved in the averaging; in this example, it is set to 32. Indicates the data segment number, with a value ranging from 1 to... ; Indicates the most recent Each data segment at the phase position The sampled values ​​on; Indicates the phase position of the rotational synchronization component. The values ​​above. Rotational synchronous disturbances remain coherent at the same phase position and are retained after synchronous averaging; asynchronous random vibrations do not remain coherent at the same phase position, and their residual standard deviation is approximately calculated according to... attenuation. When the value is 32, the residual standard deviation is approximately 17.7% of the original value. Figure 2 The average result of 32-cycle synchronization and the mean of synchronization components are shown. The value of affects the ability to suppress random vibrations and the speed of tracking operating conditions; at 12 r / min, 32 cycles correspond to 160 s, which can cover multiple random vibration cycles while avoiding the mixing of rotational synchronization components with slow-changing states over an excessively long time range. The first 32 cycles after startup are averaged using all accumulated data segments. It equals the number of accumulated cycles.

[0027] The rotation synchronization component consists of a periodic waveform and a DC component formed by the attitude average over the most recent 32 cycles. If the rotation synchronization component is directly subtracted from the latest data segment, slow attitude drift will also be subtracted. First, calculate the phase average of the rotation synchronization component. ,in This represents the arithmetic mean of the rotation synchronization component at 1024 phase positions; then, the rotation synchronization component is subtracted point by point. The zero-mean rotating synchronization perturbation was obtained. ,in Indicates phase position Zero-mean rotating synchronization perturbation. Processed Only the periodic waveform is retained, and the DC attitude component is preserved in the original tilt angle signal.

[0028] The zero-mean rotating synchronization disturbance values ​​at the same phase position are subtracted point by point from the 1024 sampled values ​​of the latest data segment to obtain 1024 corrected sampled values. These corrected sampled values ​​are then sorted numerically, and the 51 largest and 51 smallest values ​​are removed. The remaining 922 corrected sampled values ​​are then arithmetically averaged, and this average is used as the imbalance characterization value for the current cycle. The two-sided truncated averaging method eliminates the direct involvement of phase domain disturbances in state estimation and weakens the impact of local shocks and anomalous sampled values ​​on the cycle characterization. The resulting imbalance characterization value retains the tilt angle dimension and characterizes the slowly varying attitude shift of the current cycle. The imbalance characterization sequence is calculated cycle by cycle. In the example data, the peak-to-peak fluctuation of the original tilt angle signal is approximately 0.4°, while the periodic fluctuation of the processed imbalance characterization sequence is less than 0.03°. Figure 3 The thin solid line in the figure represents the average original tilt angle for each period, while the thick solid line represents the imbalance characterization sequence.

[0029] In one implementation, a causal moving median filter of length 5 is applied to the imbalance characterization values ​​of the current cycle and the four cycles preceding it, and the median of the five values ​​is taken as the current output. This processing can suppress single-cycle anomalies caused by instantaneous impacts or sudden load changes on the workpiece, and does not use data after the current cycle. The number of isophase pulses output per revolution of the rotary phase sensor can be 512, 1024, or 2048; the sampling frequency of the dual-axis tilt sensor can be selected in the range of 100Hz to 500Hz, and cubic spline interpolation is used when the sampling frequency is below 100Hz. The number of data segments for synchronous averaging can be adjusted in the range of 16 to 64; the number of data segments is increased when random vibration is strong, and decreased when slow attitude offset changes rapidly. The point-by-point arithmetic mean used to construct the rotational synchronization component can also be replaced with a point-by-point truncated average, that is, after deleting the two largest and two smallest values ​​from the 32 sampled values, the average of the remaining 28 values ​​is calculated. The attitude controller uses a circular buffer to recursively update the synchronous averaging results. When the time interval between two adjacent equal-phase pulses exceeds 2 seconds, the rotating working component is determined to be paused, the data segment division and synchronous average update are frozen, and execution continues after the equal-phase pulses are restored.

[0030] When performing integrity checks on equal-phase pulses, the attitude controller simultaneously compares the median of the pulse interval in the current cycle with the median of the previous 8 valid cycles. When the difference exceeds 20%, the current cycle is marked as a speed transition cycle; equal-phase resampling can still be performed according to the timestamp during the transition cycle, but it is not immediately used to update the rotational synchronization component. After the difference in the median pulse intervals for three consecutive cycles is no greater than 5%, synchronization averaging updates are resumed. This process prevents transient periodic waveforms during acceleration or deceleration phases from contaminating the ring buffer.

[0031] When no more than three consecutive sampling points are missing from the tilt sampling point within a rotation cycle, linear interpolation is performed using the valid sampling values ​​on both sides of the missing point. When more than three consecutive sampling points are missing, or when the number of missing points within a week exceeds 1% of the total number of samples, the cycle is marked as an invalid data segment. Invalid data segments do not participate in synchronous averaging, sliding median filtering, and slow feature updates. The attitude controller uses the dual-axis slow-change attitude state from the previous valid cycle and reduces the available control time according to the most unfavorable drift coefficient during execution.

[0032] The imbalance characterization sequence obtained in step 1 still includes synchronous average residuals, tilt sensor quantization noise, and mechanical vibration, with periodic fluctuations of approximately 0.02° to 0.03°. Direct extrapolation based on adjacent period values ​​will cause the drift estimate to change sign with residual vibration. Slow feature analysis searches for the projection direction with the minimum difference variance under normalized variance conditions within the time-delay sample space, thereby extracting attitude states with a rate of change lower than that of residual vibration. This method does not preset a fixed time-domain cutoff frequency and is suitable for operating conditions where the frequency bands of rotating synchronous disturbances and asynchronous random vibrations overlap.

[0033] During implementation, each rotation cycle is considered as one time step, and the imbalance characteristics of the two attitude axes over the most recent 120 time steps are extracted to form a sliding window. When the rotation speed is 12 r / min, 120 time steps correspond to 600 s, which can cover multiple cycles of residual vibration; after a new operating state appears, the original window data is replaced within 120 rotation cycles.

[0034] For any attitude axis, starting from the 6th time step of the sliding window, the current imbalance representation value and the imbalance representation values ​​of the previous 5 time steps are arranged in chronological order to form a 6-dimensional sample vector. ,in Indicates the time step number within the sliding window. Indicates the first Imbalance representation values ​​at each time step, superscript This represents the vector transpose. 120 time steps form 115 6-dimensional sample vectors. A 5th-order time delay covers the residual vibration correlation lengths of 3 to 5 cycles in the example system, and ensures that the number of samples is greater than 15 times the vector dimension.

[0035] For each of the 115 6-dimensional sample vectors, the mean is removed dimension by dimension, and the denoted vector is... The mean of each dimension is , The dimension number is represented by a value from 1 to 6; the mean-removed sample vector is still denoted as... Calculate the covariance matrix of the sample vectors. , It is a 6×6 symmetric matrix, and its first... Line number The column element is the first The dimension and the first Covariance in each dimension Indicates the dimension index. (For...) Perform eigenvalue decomposition to obtain eigenvalues ​​arranged in descending order of numerical value. to and corresponding feature vectors to ; Indicates the first 1 eigenvalue, Indicates and The corresponding feature vector. Whitening results are calculated as follows: Calculation, where Represents the whitened sample vector The Each component. After dividing by the square root of the eigenvalue, the variance of each principal direction is normalized to the same order of magnitude, and the subsequent difference variance only reflects the rate of change over time.

[0036] When time delay components are highly correlated, the smallest eigenvalue of the covariance matrix may be close to 0. For example... for When it is one part per million, divide directly by This will amplify quantization noise. Set a lower bound for regularization. eigenvalues ​​less than At that time, with Instead of this feature value, it participates in whitening. This embodiment sets... equal of The regularization lower bound is multiplied by a factor of 1, which makes the regularization lower bound vary with the signal variance.

[0037] Calculate the first-order difference of the whitened sample vectors at adjacent time steps. This yields 114 difference vectors; Indicates the first The time step and the first The difference between the whitened sample vectors at each time step. Calculate the covariance matrix of the difference vectors. And perform feature decomposition, It is a 6×6 symmetric matrix. Take... The eigenvector corresponding to the smallest eigenvalue is used as the slow direction vector. In the example data, The six eigenvalues ​​are 2.61, 1.87, 1.42, 0.96, 0.39, and 0.011, respectively, with the smallest eigenvalue being significantly smaller than the second smallest. Projecting the whitened sample vector onto the slow direction vector yields... The obtained values ​​are arranged according to time steps to form a slow feature sequence. Indicates the first Slow feature values ​​at each time step.

[0038] The whitened slow features do not retain the original tilt angle dimensions, and the feature vectors and This represents the same characteristic direction. To restore the dimensions and direction, within the same sliding window, a least-squares linear fit is performed with the slow characteristic value as the independent variable and the imbalance representation value of the corresponding time step as the dependent variable to obtain the slope. and intercept and in accordance with Slow conversion features, among which Indicates the slope. Indicates the intercept. Indicates the first The slowly varying estimates at each time step. When the sign of the eigenvector changes, the sign of the fitted slope changes synchronously, so the transformed slowly varying estimates remain consistent with the slope direction. The example data has a slope of 0.047° and an intercept of 0.28°, and after transformation, the standard deviation of the periodic fluctuations decreases from 0.011° to 0.003°.

[0039] The sliding window length can be selected from 60 to 240 time steps, the time delay order can be selected from 3 to 9, and the number of effective samples should be no less than 15 times the dimension of the sample vector. When the ratio of the smallest eigenvalue to the second smallest eigenvalue of the difference covariance matrix is ​​greater than 0.5, the projection results of the two candidate slow directions are calculated separately, and the result with the larger absolute value of the correlation coefficient with the imbalanced representation sequence is selected. The lower bound of regularization can also be a fixed value. Square degree, or uniformly add it to the diagonal of the covariance matrix. Then, feature decomposition is performed. Feature decomposition and linear fitting are updated every 5 periods, and the most recent effective parameters are used in the remaining periods. When the absolute value of the fitted slope is less than 0.001°, the most recent effective slope and intercept are used until the correlation is recovered.

[0040] The attitude controller constructs a sliding window using the same effective period index for both attitude axes. If the tilt angle data for either attitude axis is invalid in a certain period, neither attitude axis will include that period in the time delay embedding, thus ensuring that the two components of the dual-axis slowly varying attitude state correspond to the same physical moment. When the number of effective time steps within the window is less than 60, the slow direction vector is not updated; when the number of effective time steps reaches 60 to 119 and the number of effective samples is not less than 15 times the dimension of the sample vector, a temporary slowly varying estimate is calculated using all currently valid data, and if the lower limit of the number of samples is not reached, the most recent valid result is used; after the number of effective time steps reaches 120, the system enters the rated sliding window operation.

[0041] Regularized whitening and eigenvalue decomposition of the difference covariance matrix are performed once every 5 rotation cycles. In other cycles, the most recently valid whitening matrix and slow direction vector are used, and only the time-delay samples and their projection values ​​are updated. When the variance of the imbalance representation within the sliding window changes by more than 25% compared to the previous eigenvalue decomposition, or when the ratio of the smallest eigenvalue to the second smallest eigenvalue changes by more than 0.2, eigenvalue decomposition is performed earlier. This update method reduces the computational load on the attitude controller while maintaining slow feature tracking capabilities.

[0042] When restoring the angular dimensions of slow features, in addition to calculating the linear fitting slope and intercept, the root mean square (RMS) of the fitting residuals is also calculated. If the RMS is no greater than 35% of the standard deviation of the corresponding attitude axis imbalance characterization, the current slow-varying estimation sequence is accepted; if it exceeds 35%, the previous effective fitting slope and intercept are retained, and the current sliding window is marked as low confidence. In the low confidence state, the attitude constraint boundaries are not amplified; instead, the preset safety margin is increased to 1.5 times the rated value to allow sufficient execution time.

[0043] The slowly varying estimates of the two attitude axes are combined to form a dual-axis slowly varying attitude state. The attitude allowable domain is described by eight directed boundaries, including positive and negative boundaries of the first and second attitude axes, and four oblique boundaries formed by a normalized, equally weighted combination of the two attitude axes. In the oblique projection, the absolute value of the weights of the two attitude axes is the reciprocal of the square root of 2. In this embodiment, the axial boundary is 0.5° and the oblique projection boundary is 0.6°. The oblique boundaries can reflect the combined attitude risk in advance when both attitude axes shift simultaneously but a single axis has not yet reached 0.5°. For each attitude constraint direction, the dual-axis slowly varying attitude state is projected onto that direction to obtain a one-dimensional direction projection sequence. The adjacent increments of the direction projection sequence consist of systematic drift and random diffusion. In this embodiment, a Wiener process with drift is used to calculate the first-arrival probability of the corresponding attitude constraint boundary, and the minimum predicted first-arrival time among the eight directions is used as the control triggering basis. Figure 6 The octagon in the diagram represents the biaxial attitude allowable domain, which is jointly defined by four axially directed boundaries and four obliquely directed boundaries. The dotted lines from the center to each boundary represent the corresponding attitude constraint directions.

[0044] For any attitude constraint direction, let the increment of adjacent time steps in the direction projection sequence be denoted as . ,in Indicates the first The time step relative to the first The increment of each time step Indicates the first The directional projection values ​​at each time step. 115 directional projection values ​​form 114 increments. Original drift coefficient. Take the arithmetic mean of 114 increments; diffusion coefficient We take the standard deviation of 114 increments, and use 113 as the denominator when calculating the sample variance. In the example data... It is 0.0012° / cycle. The mean drift is 0.004°; the standard error of the drift mean is 0.004° divided by the square root of 114, which is approximately 0.000375° / cycle.

[0045] For an attitude constraint direction, let the current direction projection value be... The corresponding original drift coefficient is Before calculating the first arrival time, if the original drift coefficient is positive and less than the preset drift lower limit, the preset drift lower limit is used to replace the original drift coefficient; otherwise, the original value is retained. The drift coefficient referred to below is the drift coefficient used for prediction. Current direction projection value. The difference to the attitude constraint boundary in that direction is denoted as .when When the value is greater than 0, the systematic drift points towards the boundary in that direction; when it is not greater than 0, the conservative first-arrival probability is still calculated based on the diffusion term. If the boundary has been reached or crossed, then Set it to 0, and directly set the predicted first arrival time in that direction to 0; otherwise, This represents the angular margin from the current attitude state to the corresponding directed boundary.

[0046] The predicted first arrival time is calculated in segments based on boundary distance, drift coefficient, and diffusion coefficient. Greater than 0 and When the diffusion threshold is less than the preset lower limit, Divide by Multiplying this by the current rotation period gives the deterministic first arrival time; when Not greater than 0 and When the diffusion threshold is less than the preset lower limit, the attitude state exhibits no tendency drift and no discernible diffusion, so the predicted first arrival time is set to 86400 s. Not less than the preset diffusion lower limit and When the value is not greater than 0, the zero-drift Wiener process is used as a conservative approximation, and its first-reach cumulative distribution function is: The cumulative probability is set to 0.5, and the median time step number is obtained using a bisection method. When the prediction drift coefficient is not greater than 0, the original drift is non-positive, and the zero drift approximation will not underestimate the hit probability. When the diffusion coefficient is not less than the preset diffusion lower limit and the prediction drift coefficient is greater than 0, Divide by Obtain the mean parameters of the inverse Gaussian distribution. ,Will Squared by The square of the shape parameter is obtained. . This represents the number of time steps required to reach the boundary based on the average drift. It reflects the concentration of the first arrival time; When the size increases, the first arrival time becomes more concentrated in .by and The cumulative distribution function of the inverse Gaussian distribution with parameter is: ,in Indicates the number of candidate time steps. This indicates that the Shouda incident occurred in The probability of occurrence within a time step The cumulative distribution function of the standard normal distribution. It is an exponential function. about Monotonically increasing, with 0.001 times is the initial lower bound of the binary search. The upper bound; because A value greater than 0.5 and 30 iterations yield the median number of time steps. Multiplying the median by the actual duration of the most recent complete rotation cycle gives the predicted first arrival time. The preset drift lower limit is 0.0001° / cycle, and the preset diffusion lower limit is 0.0001°.

[0047] For each directional attitude constraint direction, the drift coefficient and diffusion coefficient of the most recent 15 sliding windows are maintained separately. If the current diffusion coefficient exceeds three times the historical median and a single-cycle jump occurs in the imbalance characterization sequence, the window is identified as an impact window and is not used to update the long-term diffusion reference value. If the diffusion coefficient of three consecutive windows exceeds twice the historical median, it is considered a change in operating conditions rather than an instantaneous impact, and the first arrival time is calculated based on the new diffusion level.

[0048] The median of the inverse Gaussian distribution described above is solved using a binary search. The initial search lower bound is set to 0.001 times the mean parameter, and the search upper bound is also set to the mean parameter. If the cumulative probability at the initial lower bound is not less than 0.5, the lower bound is successively divided by 10 until the cumulative probability at the lower bound is less than 0.5. When the lower bound approaches 0, it is treated as having a cumulative probability of 0. The cumulative distribution function is calculated by taking the arithmetic mean of the upper and lower bounds each time. The upper bound is updated when the cumulative probability is not less than 0.5; otherwise, the lower bound is updated. The process ends when the difference between the upper and lower bounds is less than 0.01 time steps or the number of iterations reaches 30, and the average of the upper and lower bounds is used as the median time step. The second term of the inverse Gaussian cumulative distribution function is calculated using the logarithmic domain to avoid floating-point overflow or underflow caused by multiplying the exponential term with the normal tail probability. The zero-drift Wiener process uses the same cumulative probability and convergence criterion. The upper limit of the search is the number of time steps obtained by dividing the preset maximum prediction time by the current rotation cycle duration. If the cumulative probability at the upper limit of the search is still less than 0.5, the first prediction time is taken as the preset maximum prediction time.

[0049] Taking a boundary distance of 0.16°, a drift coefficient of 0.0035° / cycle, and a standard deviation of 0.004° per cycle as an example, the mean parameter of the inverse Gaussian distribution is 45.7143 time steps, the shape parameter is 1600 time steps, and the median of the cumulative distribution function is 45.0719 time steps. With a current rotation cycle length of 5 seconds, the predicted first arrival time is 225.3594 seconds, rounded to two decimal places, which is 225.36 seconds. The attached chart shows the rounded value as 225 seconds. The attitude controller retains the unrounded value for trigger comparison.

[0050] The predicted first arrival time determines the available control time, and the target execution quantity needs to be generated based on the actual response of the attitude adjustment actuator. The control effectiveness matrix is ​​a 2×2 matrix, with the two rows corresponding to the tilt angle changes of the first and second attitude axes, respectively, and the two columns corresponding to the unit changes of the first and second attitude axis control modes, respectively; off-diagonal elements represent the cross-coupling between the two control modes. During calibration, the rotating working part is kept stationary, and positive and negative execution quantities are applied to the first and second attitude axis control modes, respectively; the modal amplitude in each direction is successively taken as 2 mm, 4 mm, 8 mm, and 16 mm, and converted into four hydraulic cylinder commands with a sum of 0 extension and retraction increments according to the support position coordinates. After each execution is completed and the rate of change of the dual-axis tilt angle is less than 0.002° / s for 10 consecutive seconds, the change of the dual-axis tilt angle is recorded, and the initial value of the control effectiveness matrix is ​​obtained by least squares fitting. Positive and negative data are used together for fitting to reduce the influence of static zero bias and mechanical backlash on the matrix coefficients.

[0051] The attitude controller combines the first and second attitude axis control modes in 2mm increments within a range of -16mm to 16mm, forming 289 initial candidate modes. Based on the coordinate symbols of the four support positions, the two control modes are converted into four hydraulic cylinder extension / retraction increments, and candidates whose sum of the four extension / retraction increments is not zero or whose post-execution position exceeds the ±25mm working range are deleted. For the retained candidates, a speed-limited displacement trajectory is generated based on the extension / retraction increment of each hydraulic cylinder and the current allowable speed; the four electro-hydraulic support actuators move in parallel, and the candidate execution delay is the maximum value of the time required for each speed-limited displacement trajectory. The attitude controller extrapolates the unadjusted attitude within the execution delay according to the original drift coefficient without lower limit correction, and then superimposes the biaxial tilt angle change given by the control effectiveness matrix to obtain the predicted attitude after execution. The margin from the predicted attitude to the boundary is calculated for each of the eight attitude constraint directions; if the margin is less than the preset control margin, the square of the deficiency is penalized. The control objective function also includes the execution energy consumption calculated from the load measured by the pressure sensor, the effective area of ​​the hydraulic cylinder, and the extension / retraction increment, as well as the square of the change in the candidate execution quantity relative to the previous execution quantity. The candidate with the smallest control objective function is taken as the target execution quantity.

[0052] When generating candidate execution quantities, the selectable range of the two attitude axis control modes is first narrowed down based on the current displacement and allowable working range of each hydraulic cylinder, and then discretely combined. For each group of discrete control modes, common mode verification, stroke verification, generation of speed-limited displacement trajectory according to allowable speed, and pressure margin verification are performed sequentially. If any verification fails, the subsequent calculation of that candidate is immediately terminated; after all verifications pass, the duration of the longest speed-limited displacement trajectory is used as the execution delay, and the predicted attitude after execution is calculated. This sequence can reduce matrix operations and multi-directional constraint projections performed on infeasible candidates.

[0053] The predicted attitude constraint margin penalty term is calculated for each of the eight directed constraint margins. When the margin is positive, its ratio to the corresponding attitude constraint boundary is taken as the dimensionless safety margin; when the dimensionless safety margin is lower than the preset dimensionless control margin, a penalty is applied to the square of the shortfall. When the margin is 0 or negative, a high-weight penalty is applied to the square of the out-of-bounds amount. The execution energy consumption term is obtained by summing the products of the absolute values ​​of the pressure difference between the two chambers of each hydraulic cylinder, the effective piston area, and the extension displacement; the execution quantity change term is obtained by summing the squares of the differences between the current candidate control mode and the previously executed control mode. Each term is first normalized to the reference value during the calibration phase, and then weighted and summed to avoid any one term becoming dominant due to dimensional differences.

[0054] When multiple candidate values ​​with the same control objective function appear in the same solution, the candidate with the shortest execution delay, the smallest maximum single-cylinder displacement, and the smallest Euclidean distance from the previous execution value is selected in sequence. When all candidate execution values ​​are eliminated, the displacement step size is reduced from 2 mm to 1 mm and the solution is resolved; if there are still no feasible candidates, the current position of each hydraulic cylinder is maintained, and a stroke or rate constraint alarm is generated. Execution values ​​that exceed the constraints are not considered as alternative candidates.

[0055] The weights of the control objective function are determined during the calibration phase. First, the controlled platform is adjusted to near the center of the allowable domain. Then, each control mode is applied to obtain the reference distribution of the minimum attitude constraint margin, hydraulic energy consumption, and changes in the executed quantity after execution. Next, the three normalized values ​​are made to have similar orders of magnitude under rated operating conditions, and the weight of the attitude constraint margin penalty term is increased to ensure that candidate solutions preferentially satisfy all directed boundaries. After weight settings are completed, verification is performed using two load states not involved in the settings. If the predicted attitude after execution in any verification state exceeds the allowable domain, or the deviation between the actual attitude and the predicted attitude exceeds 10% of the corresponding axis boundary value, the control effectiveness matrix and weights are recalibrated, and the verification is not passed by relaxing the attitude constraint boundaries.

[0056] Control triggering is determined by the predicted first arrival time and the current boundary state. The target execution quantity is executed when the minimum predicted first arrival time in the eight attitude constraint directions is not greater than the sum of the execution delay of the target execution quantity and the preset safety margin; when the current projection value in any direction has reached the attitude constraint boundary, the target execution quantity is executed immediately without waiting for the probability model to be updated. Figure 4The current directional projection value for the shown condition is 0.34°, corresponding to a boundary value of 0.5° and a boundary distance of 0.16°. The drift coefficient is 0.0035° / cycle, and the standard deviation of the increment per cycle is 0.004°. This yields the inverse Gaussian distribution mean parameters over 45.7143 time steps, shape parameters over 1600 time steps, and median over 45.0719 time steps. With a cycle length of 5 seconds, the predicted first arrival time is 225.36 seconds. The execution delay of the target execution quantity is 120 seconds, and the preset safety margin is 120 seconds, totaling 240 seconds, which satisfies the prediction triggering condition. Figure 5 The corresponding first arrival time cumulative distribution function is shown, with the median time step rounded down to 45.

[0057] Half of the first attitude axis control mode is assigned with the same sign to the first and second electro-hydraulic support actuators, and the other half is assigned with the opposite sign to the third and fourth electro-hydraulic support actuators. Similarly, half of the second attitude axis control mode is assigned with the same sign to the first and third electro-hydraulic support actuators, and the other half is assigned with the opposite sign to the second and fourth electro-hydraulic support actuators. After algebraically adding the two control modes at each actuator, the sum of the extension / retraction increments of the four hydraulic cylinders is always 0. Figure 7 The diagram shows the distribution relationship of the extension and retraction increments of the two attitude axis control modes to the four hydraulic cylinders and the constraint that the common mode is 0.

[0058] Taking the first attitude axis control mode of 10mm and the second attitude axis control mode of 6mm as examples, the hydraulic cylinder extension increments of the first to fourth electro-hydraulic support actuators are 8mm, 2mm, -2mm, and -8mm respectively, and their sum is 0. The allowable working range of the example electro-hydraulic support actuator is -25mm to 25mm, and the current displacement is 0mm. All four target displacements meet the stroke constraints. When the allowable speed of the actuator is 0.0667mm / s, the time required for the parallel movement of the four actuators is determined by the actuator with the largest absolute extension of 8mm, and the calculated value is 119.94s. The controller includes 120s in the prediction. When the support spacing of the first attitude axis is 4m, the small angle change corresponding to the 10mm height difference between the front and rear supports is approximately 0.143°; when the support spacing of the second attitude axis is 3m, the angle change corresponding to the 6mm height difference between the two supports is approximately 0.115°, and the values ​​are consistent with the order of magnitude of the calibrated control performance matrix.

[0059] Candidate control modes are first quantized to integer multiples of twice the displacement sensor resolution during combination, and then... Figure 7The allocation relationship shown calculates the extension / retraction increments of the four hydraulic cylinders, ensuring that each increment is an integer multiple of the displacement resolution. After the conversion, the attitude controller recalculates the algebraic sum of the four extension / retraction increments. If the absolute value of the algebraic sum generated by floating-point arithmetic is no greater than 0.01 mm, the extension / retraction increment of the fourth electro-hydraulic support actuator is set to the negative of the sum of the first three increments; this adjustment must not exceed half of the displacement resolution, and the stroke and rate constraints are re-verified. If the absolute value of the algebraic sum exceeds 0.01 mm or the re-verification fails, the candidate is discarded. The final four extension / retraction increments satisfy the common modal constraints.

[0060] Before issuing the target execution quantity, a consistency check is performed on the latest data timestamps of the dual-axis tilt sensor, four displacement sensors, and eight pressure measuring points. If the difference between any timestamp and the current control time exceeds 200 milliseconds, the issuance is paused. After issuance, the four hydraulic cylinders start simultaneously and generate speed-limited displacement trajectories respectively, making the arrival times as consistent as possible. If the displacement tracking error of any actuator exceeds 0.5 mm for two consecutive seconds, or the deviation between the actual speed and the commanded speed exceeds 20% for two consecutive seconds, all actuators are stopped and the current position is retained; this response is not used to update the control performance matrix.

[0061] After the target execution is completed, the attitude controller calculates the tilt rate of change for the two attitude axes using the first-order difference of adjacent tilt angle sample values. If both rates of change are less than 0.002° / s for 10 consecutive seconds, the execution response is considered stable. The actual control response is obtained by subtracting the original drift coefficients (without lower bound correction) calculated from the slow-varying estimation sequences of the two attitude axes and the natural drift calculated based on the execution duration from the dual-axis attitude changes before and after execution. The control effectiveness matrix is ​​updated using recursive least squares with a forgetting factor of 0.995: a regression vector is formed by the first and second attitude axis control modes; the difference between the actual control response and the matrix-predicted response is used as the innovation; the gain is calculated based on the parameter covariance matrix and the regression vector from the previous cycle; and the control effectiveness matrix and parameter covariance matrix are then updated accordingly. The control effectiveness matrix is ​​not updated if the L2 norm of the regression vector is less than 2 mm, the actual control response amplitude is less than three times the standard deviation of the tilt noise, or the tilt rate of change has not yet entered a stable range. The absolute value of the deviation of each element of the control performance matrix from the initial calibration value must not exceed 60% of the absolute value of the corresponding initial calibration value, and the updated value must retain the same sign as the initial calibration value. Elements whose absolute value of the initial calibration value is less than the preset zero threshold are subject to a fixed symmetrical limiting interval centered at 0. After execution, the slow feature analysis sliding window is re-accumulated; if the window has not been formed, the imbalance characterization value continues to be calculated, and the current boundary state is used as the protection trigger condition. The above steps are repeated until the attitude control task ends.

[0062] Before recursive least squares update, the actual control response is checked for amplitude and direction. If the actual control response exceeds 2.5 times the predicted response of the control performance matrix on any attitude axis, or if the angle between the main response direction and the predicted direction exceeds 60°, the data is marked as an anomalous response. Anomalous responses are only used to generate support state check tasks and do not participate in matrix updates. If two consecutive anomalous responses in the same direction occur, small-amplitude forward and reverse calibrations must be re-executed to distinguish between changes in control performance and sensor anomalies.

[0063] The parameter covariance matrix used for recursive least squares is initialized as a diagonal matrix, with diagonal elements set according to the control mode amplitude during the calibration phase. After each update, the parameter covariance matrix is ​​symmetricized, and the eigenvalues ​​are limited to a preset positive range to prevent finite-precision calculations from compromising positive definiteness. When the updated value reaches the control performance matrix limit boundary, a parameter limiting event is recorded.

[0064] The attitude controller stores the trigger time of each control event, the current projection values ​​in eight directions, the original drift coefficient, the predicted drift coefficient, the diffusion coefficient, the predicted first arrival time, the target execution amount, the execution delay, the biaxial tilt angles before and after execution, and the matrix update results. Records are indexed and associated by rotation period, and the phase domain processing, probabilistic prediction, and execution feedback processes are replayed using the same time base. The operation log also retains invalid data segments, sensor timeouts, candidate solutions that are not found, and parameter limiting flags to determine the specific reasons why control was not executed or the model was not updated.

[0065] After a complete control event, the attitude controller verifies the execution results. Following execution, the dual-axis slow-change attitude state must be within the eight directed attitude constraint boundaries, the minimum attitude constraint margin must be greater than the preset control margin, the final displacement of the four hydraulic cylinders must be within the working range, and the common modal residual must not exceed 0.01 mm. When all the above conditions are met, the event state is recorded as complete. If the attitude constraint margin is insufficient but has not yet exceeded the limits, candidate execution quantities are regenerated using the latest state. If a limit violation occurs, the pressure reaches the rated value, displacement tracking fails, or the sensor times out, the event state is recorded as aborted, the current stable support position is maintained, and a manual inspection request is output. The response data for aborted events is not used for updating the control performance matrix.

[0066] The above embodiments are used to illustrate the technical solution of the present invention and do not constitute a limitation on the scope of protection. When there are 3 or more attitude adjustment actuators, a geometric allocation matrix is ​​constructed based on the support position coordinates of each actuator relative to the two attitude axes. Under the condition that the common mode is 0, the control modes of the two attitude axes are converted into the execution quantities of each actuator using the least second norm solution; the control performance matrix is ​​still used to describe the relationship between the control modes of the two attitude axes and the change in the biaxial tilt angle. Those skilled in the art can also make equivalent substitutions for the number of sampling points, window length, attitude adjustment actuator type, and candidate execution quantity solution method without departing from the technical concept of the present invention.

[0067] When there are three attitude control actuators, a 2×3 geometric allocation matrix is ​​first constructed using the coordinates of the three support positions on two attitude axes, and then solved simultaneously with the common modal constraints. The three support positions must not be collinear, and the rank of the geometric allocation matrix must be 2. When there are more than four attitude control actuators, the solution with the smallest weighted execution quantity L2 norm is selected from the solution set that satisfies the control modes of the two attitude axes and the common modal constraints. The weights can be determined based on the actuator's remaining stroke, allowable speed, and pressure margin. Actuators with smaller remaining stroke or lower pressure margins are given larger weights, thereby reducing their allocation quantity.

[0068] The attitude adjustment actuator can also be an electric cylinder, a servo screw, or other linear actuator capable of continuously adjusting the support height. After changing the actuator type, the methods for determining the stroke constraint, rate constraint, and execution energy consumption are adjusted accordingly, while the processing relationships for phase domain scrambling, slow-change attitude state extraction, multi-directional first-arrival time prediction, control performance matrix mapping, and post-execution updates remain unchanged. For electric actuators without hydraulic pressure parameters, the execution energy consumption can be determined by the integral of the drive current, bus voltage, and execution time.

Claims

1. A multi-actuator dual-axis attitude constraint predictive control method, applied to an attitude control system including a dual-axis tilt sensor, a rotation phase sensor, an attitude controller, and at least three attitude adjustment actuators, characterized in that, Includes the following steps: Step 1: Using the isophase pulses output by the rotating phase sensor as a reference, resample the first and second attitude axis tilt signals in the same phase. Extract the zero-mean rotational synchronization disturbance according to the rotation period. Perform a two-sided truncated average on the corrected sampled values ​​after deducting the disturbance to form a dual-axis imbalance characterization sequence. Step 2: Perform time-delay embedding, regularized whitening, and slow feature extraction on the dual-axis imbalance characterization sequence to obtain the dual-axis slow-changing attitude state. Step 3: Project the dual-axis slow-changing attitude state onto multiple attitude constraint directions. Determine the drift coefficient and diffusion coefficient based on the increment of the projection sequence in each direction, and calculate the distance to the corresponding attitude constraint boundary. Step 4: Generate candidate execution quantities based on the stroke constraints, rate constraints, and common modal constraints that do not cause the overall lifting and lowering of the controlled platform of the attitude adjustment actuator. Use the control performance matrix to predict the dual-axis attitude state after the candidate execution, and determine the target execution quantity based on the predicted attitude constraint margin, execution energy consumption, and execution quantity change. Step 5: Execute the target execution quantity when the minimum predicted first arrival time is not greater than the sum of the execution delay of the target execution quantity and the preset safety margin, or when the dual-axis slow-changing attitude state reaches any attitude constraint boundary. After the response stabilizes, update the control performance matrix based on the actual dual-axis attitude change and return to step 1.

2. The multi-actuator dual-axis attitude constraint prediction control method according to claim 1, characterized in that, The attitude adjustment actuators are the first to the fourth electro-hydraulic support actuators, which are installed at four support positions on the controlled bearing platform. At least three of the four support positions are not collinear. Each electro-hydraulic support actuator includes a servo valve, a hydraulic cylinder, a displacement sensor, and a pressure sensor. The actuation amount is the hydraulic cylinder extension increment. The common modal constraint is that the sum of the extension increments of the four hydraulic cylinders is equal to 0.

3. The multi-actuator dual-axis attitude constraint prediction control method according to claim 1, characterized in that, Step 1 includes: recording the timestamps of each equal-phase pulse; performing linear interpolation on the dual-axis tilt signal according to the timestamps; dividing the data segment into 1024 sampling points per revolution of the rotating component; averaging the most recent 32 data segments at the same phase position to obtain the rotation synchronization component; subtracting the mean value of the rotation synchronization component point by point to obtain the zero-mean rotation synchronization disturbance; subtracting the zero-mean rotation synchronization disturbance point by point from the latest data segment; performing a two-sided truncated average on the obtained corrected sampling value to obtain the imbalance characterization value of the corresponding attitude axis; and performing sliding median filtering on the imbalance characterization values ​​of the current cycle and the previous 4 cycles.

4. The multi-actuator dual-axis attitude constraint prediction control method according to claim 1, characterized in that, In step 2, each rotation cycle is taken as one time step. The imbalance representation values ​​of the first and second attitude axes in the most recent 120 time steps are embedded with a 5th-order time delay to form a 6-dimensional sample vector. The sample vector is mean-removed and regularized whitened. The covariance matrix of the first-order difference between adjacent whitened sample vectors is calculated. The eigenvector corresponding to the smallest eigenvalue is taken as the slow direction vector. The whitened sample vector is projected onto the slow direction vector to obtain the slow feature sequence of the corresponding attitude axis.

5. The multi-actuator dual-axis attitude constraint prediction control method according to claim 4, characterized in that, Within the same sliding window, a least-squares linear fit is performed with the slow feature sequence as the independent variable and the imbalance representation value of the corresponding time step as the dependent variable. The slope and intercept obtained from the fit are used to convert the slow feature sequence into a slow-varying estimation sequence with the same dimensions as the tilt angle signal and the same sign as the tilt direction. The slow-varying estimation sequences of the two attitude axes constitute the dual-axis slow-varying attitude state.

6. The multi-actuator dual-axis attitude constraint prediction control method according to claim 1, characterized in that, Multiple attitude constraint directions include the positive and negative directions of the first and second attitude axes, as well as four oblique directions formed by the normalized and equally weighted combination of the two attitude axes. For each attitude constraint direction, the dual-axis slowly changing attitude state is projected onto that direction, and the arithmetic mean and sample standard deviation of the increments of adjacent time steps in the projection sequence are taken as the original drift coefficient and diffusion coefficient, respectively. When the original drift coefficient is greater than 0 and less than a preset drift lower limit, the preset drift lower limit is used as the predicted drift coefficient; otherwise, the original drift coefficient is used as the predicted drift coefficient. When the diffusion coefficient is small... When the diffusion coefficient is below the preset lower limit and the predicted drift coefficient is greater than 0, the predicted first arrival time is determined by dividing the boundary distance by the predicted drift coefficient. When the diffusion coefficient is less than the preset lower limit and the predicted drift coefficient is not greater than 0, the predicted first arrival time is set to the preset maximum value. When the diffusion coefficient is not less than the preset lower limit and the predicted drift coefficient is greater than 0, the predicted first arrival time is determined according to the median of the inverse Gaussian distribution. When the diffusion coefficient is not less than the preset lower limit and the predicted drift coefficient is not greater than 0, the predicted first arrival time is determined according to the first arrival probability of the zero-drift Wiener process.

7. The multi-actuator dual-axis attitude constraint prediction control method according to claim 2, characterized in that, The control performance matrix is ​​a 2×2 matrix that maps the first attitude axis control mode and the second attitude axis control mode to the first attitude axis tilt angle change and the second attitude axis tilt angle change. During the calibration phase, the positive and negative execution quantities of the first attitude axis control mode and the second attitude axis control mode are applied respectively, and the tilt angle changes of the two axes before and after execution are recorded synchronously. The initial value of the control performance matrix is ​​determined by least squares fitting.

8. The multi-actuator dual-axis attitude constraint prediction control method according to claim 7, characterized in that, Step 4 includes: combining the first and second attitude axis control modes according to a preset displacement step size, and converting them into four hydraulic cylinder extension / retraction increments that satisfy common mode constraints; generating speed-limited displacement trajectories based on the extension / retraction increments of each hydraulic cylinder and the corresponding allowable speeds, determining the maximum time required for each speed-limited displacement trajectory as the execution delay, extrapolating the unadjusted attitude within the execution delay according to the drift coefficient, and superimposing the attitude change amount obtained by the control efficiency matrix acting on the two attitude axis control modes to obtain the predicted attitude after execution; eliminating candidate execution quantities that cause any hydraulic cylinder to exceed the stroke limit, and determining the candidate execution quantity with the smallest control objective function as the target execution quantity, wherein the control objective function is composed of a weighted sum of the predicted attitude constraint margin penalty term, the execution energy consumption term, and the change term relative to the previous execution quantity.

9. The multi-actuator dual-axis attitude constraint prediction control method according to claim 8, characterized in that, Half of the first attitude axis control mode is assigned to the first electro-hydraulic support actuator and the second electro-hydraulic support actuator with the same sign, and the other half is assigned to the third electro-hydraulic support actuator and the fourth electro-hydraulic support actuator with the opposite sign. Half of the second attitude axis control mode is assigned with the same sign to the first and third electro-hydraulic support actuators, and the other half is assigned with the opposite sign to the second and fourth electro-hydraulic support actuators; the two control modes are algebraically added at each electro-hydraulic support actuator to obtain the hydraulic cylinder extension / retraction increment from the first to the fourth electro-hydraulic support actuators.

10. The multi-actuator dual-axis attitude constraint prediction control method according to claim 8, characterized in that, After the target execution is completed and the tilt rate of change of the two attitude axes is continuously within the preset stable range, the actual control response is obtained by subtracting the natural drift during the execution period from the dual-axis attitude change before and after execution. The regression vector is composed of the first attitude axis control mode and the second attitude axis control mode. The difference between the actual control response and the predicted response of the control performance matrix is ​​used as the update residual. The control performance matrix is ​​corrected by recursive least squares with forgetting factor. When the norm of the regression vector is less than the preset excitation lower limit, the actual control response amplitude is less than 3 times the standard deviation of the tilt noise, or the tilt change rate does not enter the preset stable range, the correction of the control effectiveness matrix is ​​paused.

Images