Synergistic control method, system and apparatus for pumping tests

Abstract

The application provides a kind of pumping test synergic control method, system and device, the method includes collecting main pump control instruction and pipeline pressure, flow signal, reconstructs phase space and calculates local recurrence rate;When local recurrence rate is lower than short-term average and deviation exceeds fluctuation tolerance, it is judged to enter unstable precursor state.Subspace identification is carried out on the data before triggering, and a high-order model is obtained;Combined with Hankel singular value and eigenvalue, the minimum order reduction q is determined to meet the energy threshold and retain all unstable modes.Modal decomposition is carried out on the high-order model, and slow trend mode and key mode are identified according to wavelet sub-band energy proportion, and combined with unstable mode to obtain q-order reduced model.Based on the model, model predictive control is carried out, and the predicted output is decomposed into low-frequency trend component and high-frequency disturbance component to generate main pump reference control sequence and bypass fine adjustment valve dynamic compensation sequence respectively, to realize synergic control.

Description

Technical Field

[0001] This application belongs to the field of collaborative control, and in particular relates to a collaborative control method, system and device for pumpability testing. Background Technology

[0002] In actual pumping operations or pumpability tests, the system is susceptible to external environmental disturbances, changes in pipeline structure, and variations in the fluid medium's viscosity, gas content, particle distribution, and other physical properties, causing irregular fluctuations in pipeline pressure and flow signals. These fluctuations not only reduce the stability and reliability of pumpability test results but may also induce local or global instability under conditions of high load, long-distance transport, or sudden changes in medium state. Existing technologies typically rely on fixed threshold alarms or conventional linear state monitoring methods, which struggle to characterize the system's dynamic evolution from multi-dimensional time-series signals such as pressure and flow. This leads to a lag in identifying instability precursors, often resulting in reactive adjustments only after an instability trend has become apparent or after instability has occurred, posing response delays and operational safety risks.

[0003] Meanwhile, although pumping systems can be identified through mathematical models, the directly obtained high-order models are usually computationally intensive and difficult to meet the requirements of online prediction and real-time control. Conventional model reduction methods may weaken or lose key modal information related to system stability. On the other hand, existing control strategies mostly rely on a single actuator, the main pump, which struggles to balance long-term baseline adjustment and rapid disturbance compensation when faced with complex operating conditions involving both macroscopic low-frequency trend changes and high-frequency transient disturbances. Therefore, a collaborative control method for pumpability testing is urgently needed, capable of identifying early signs of instability, preserving key stability information while ensuring real-time model prediction, and improving operational smoothness and control reliability during pumpability testing through the synergistic effect of the main pump and auxiliary fine-tuning actuators. Summary of the Invention

[0004] To address the problems of lagging identification of instability precursors in existing pumping systems, difficulty in meeting real-time predictive control requirements of high-order models, easy loss of stability modes during order reduction, and difficulty of a single actuator in simultaneously handling low-frequency trends and high-frequency disturbances, this invention proposes a collaborative control method, system, and device for pumpability testing.

[0005] In a first aspect, the present invention proposes a cooperative control method for pumpability testing, comprising: Real-time acquisition of control commands from the main pump motor and pressure and flow signals from the pipeline; reconstruction of the phase space based on the pressure and flow signals; calculation of the local recurrence rate of the state trajectory; when the local recurrence rate is lower than the short-term mean and the absolute value of the difference from the short-term mean exceeds the fluctuation tolerance, it is determined that an unstable precursor state has been entered and subsequent steps are triggered. A high-order model is obtained by subspace identification of the pre-trigger data; the Hankel singular values ​​and eigenvalues ​​of the model are calculated; the order q of the reduced model is determined, wherein the order q is the smallest positive integer that satisfies the condition that the cumulative singular value energy is not lower than a threshold and retains all unstable modes; modal decomposition is performed on the high-order model, and the wavelet approximate subband and detail subband coefficients of the impulse response of each stable mode are calculated; the slowly varying trend mode is identified based on the energy proportion of the approximate subband coefficients, and the key mode is identified based on the energy proportion of the detail subband coefficients among the remaining stable modes not identified as slowly varying trend modes; the slowly varying trend mode, the key mode, and all unstable modes are combined to obtain the q-order reduced model; The q-order reduced model is used to generate output prediction for model predictive control and decomposed into low-frequency trend components and high-frequency disturbance components. The main pump reference control sequence is calculated based on the low-frequency trend components. The high-frequency disturbance components are used as feedforward signals to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control.

[0006] On the other hand, the present invention also proposes a cooperative control system for pumpability testing, comprising: The calculation module is used to collect control commands from the main pump motor and pressure and flow signals from the pipeline in real time; reconstruct the phase space based on the pressure and flow signals; calculate the local recurrence rate of the state trajectory; when the local recurrence rate is lower than the short-term mean and the absolute value of the difference between the local recurrence rate and the short-term mean exceeds the fluctuation tolerance, it is determined that the system has entered an unstable precursor state and triggers subsequent steps. The combination module is used to perform subspace identification on the pre-trigger data to obtain a higher-order model; calculate the Hankel singular values ​​and eigenvalues ​​of the model; determine the order q of the reduced-order model, where the reduced-order model order q is the smallest positive integer that satisfies the condition that the accumulated singular value energy is not lower than a threshold and retains all unstable modes; perform mode decomposition on the higher-order model, calculate the wavelet approximate subband and detail subband coefficients of the impulse response of each stable mode; identify the slow-varying trend mode based on the energy proportion of the approximate subband coefficients, and identify the critical mode among the remaining stable modes not identified as slow-varying trend modes based on the energy proportion of the detail subband coefficients; and combine the slow-varying trend mode, the critical mode, and all unstable modes to obtain the q-order reduced-order model. The control module is used to generate output predictions by using the q-order reduced model for model predictive control, and decomposes it into low-frequency trend components and high-frequency disturbance components; calculates the main pump reference control sequence based on the low-frequency trend components; and uses the high-frequency disturbance components as feedforward signals to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control.

[0007] This invention identifies potential instability precursors in pumping systems by reconstructing the phase space of pressure and flow signals and combining the deviation relationship between local recurrence rate and short-term mean, thereby improving the ability to detect abnormal trends in advance. After the precursor is detected, a higher-order model is obtained using subspace identification, and the order of the reduced model is determined by combining Hankel singular values ​​and eigenvalues. While retaining unstable modes, the main stable modes are selected, reducing the burden of online prediction. Based on the energy proportions of wavelet approximation subbands and detail subbands, slowly varying trend modes and key disturbance modes are distinguished, ensuring that the reduced-order model takes into account both low-frequency evolution and high-frequency response. Subsequently, the predicted output is decomposed into low-frequency trend components and high-frequency disturbance components, generating the main pump reference control sequence and bypass fine-tuning valve compensation sequence, respectively, improving the control stability and operational reliability of the pumping test process. Attached Figure Description

[0008] Figure 1 A flowchart of a coordinated control method for pumpability testing; Figure 2 For mutual information computation graph; Figure 3 This is a distribution chart of singular values. Detailed Implementation

[0009] The technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of this application, and not all embodiments. Based on the embodiments of this application, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of this application.

[0010] In the first embodiment, the present invention proposes a cooperative control method for pumpability testing, such as... Figure 1 As shown, it includes: S1 collects control commands from the main pump motor and pressure and flow signals from the pipeline in real time.

[0011] An industrial-grade data acquisition card is used, communicating with a programmable logic controller (PLC) via an EtherCAT bus. It reads the pulse width modulation (PWM) control commands from the main pump motor. Simultaneously, analog voltage signals are acquired through a piezoresistive pressure sensor and an electromagnetic flowmeter installed on the pumping pipeline, and converted into digital signals via an analog-to-digital converter. A pre-allocated fixed-length annular buffer stores the time-series data. The sampling frequency is set to 1000Hz, and a windowed moving average filtering algorithm is applied to reduce high-frequency electromagnetic noise from the sensors.

[0012] S2, reconstruct the phase space based on the pressure and flow signals.

[0013] For the collected pressure and flow time series, the mutual information method is used to calculate the first minimum value to determine the optimal delay time. A spurious nearest neighbor algorithm is used to calculate the dimension at which the proportion of neighbor points is lower than a preset neighbor threshold, which is then used as the optimal embedding dimension.

[0014] Based on the embedding theorem, a one-dimensional time series is reconstructed into a multi-dimensional state space vector matrix, which is used to characterize the evolution relationship of the system's dynamic characteristics in phase space.

[0015] As an optional implementation, the reconstruction of the phase space based on the pressure and flow signals includes: The collected pressure and flow signals are subtracted from their corresponding time-series average values, and the subtraction results are divided by the corresponding standard deviations to obtain the normalized pressure sequence and the normalized flow sequence. The delay time of the phase space is determined by calculating the first minimum point using the mutual information method, and the embedding dimension of the phase space is determined by the false nearest neighbor method when the proportion of false neighbors drops below the preset neighbor threshold. Using the determined delay time and embedding dimension, the normalized pressure sequence and normalized flow sequence are reconstructed in a multivariate phase space. The resulting state vectors each contain a time series phase space matrix of pressure and flow data corresponding to the number of dimensions specified by the embedding dimension.

[0016] When performing Z-score standardization on the collected signals, if the sampling frequency is 1000Hz and the data is extracted for a continuous acquisition time of T=10 seconds, the pressure sequence and flow sequence with a length of N=10000 are obtained. The sample mean and standard deviation of the two sets of data are calculated respectively, and zero mean and unit variance transformation is performed.

[0017] The mutual information method is used to calculate the delay time τ of normalized data. It iterates through the delay steps (e.g., 1 to 100 steps) and calculates the mutual information I(τ). The first τ value that reaches a local minimum is then found; for example, τ = 12, corresponding to a 12ms delay. The curve showing the change in mutual information with delay time is shown below. Figure 2 As shown.

[0018] After determining the delay time, the embedding dimension m is calculated using the spurious nearest neighbor method. Initially, m=2. The relative change in the nearest neighbor distance of all phase space points is calculated when the dimension is increased to m+1. If the distance increase ratio exceeds a preset threshold, for example, set to... =15, or the distance itself exceeds the tolerance limit, for example, A=2, then it is judged as a false neighbor; gradually increase m until the proportion of false neighbors decreases and stabilizes below the preset neighbor threshold, such as 1%, and determine that the embedding dimension m=4.

[0019] Based on the obtained τ=12 and m=4, the normalized pressure sequence and flow sequence are transformed into a multivariate reconstructed phase space matrix X, where each state vector has a dimension of 2m=8, that is, 4-dimensional pressure features and 4-dimensional flow features are alternately spliced. The size of this reconstruction matrix is ​​8×(N-(m-1)τ), thereby realizing the mapping of the system evolution trajectory in the high-dimensional phase space.

[0020] S3, calculate the local recurrence rate of the state trajectory.

[0021] Using a recursive quantitative analysis method, the Euclidean distance between state vectors in phase space is calculated, and a distance matrix is ​​constructed. A fixed spatial threshold is set, and a step decision function is used to transform the distance matrix into a binary recursive graph.

[0022] Within a preset sliding time window, the density of non-zero elements is statistically analyzed. By calculating the ratio of the total number of recurrence points to the total number of computable paired state points within this window, the local recurrence rate sequence is obtained.

[0023] As an optional implementation, the calculation of the local recurrence rate of the state trajectory includes: Extract the phase space moving window where the current time step is located, traverse any two non-identical state vectors in the phase space moving window, and calculate the Euclidean distance between the two state vectors. When the Euclidean distance is less than a preset threshold radius, a relapse state point is recorded. The total number of recurring state points within the phase space moving window is counted, and the total number of recurring state points is divided by the total number of calculable paired state point combinations within the phase space moving window to obtain the local recurrence rate of the state trajectory.

[0024] For the constructed time series phase space matrix, a sliding window mechanism is used to extract the instantaneous changes of the system dynamic characteristics. The length of the phase space moving window is set to W=500 sampling points, and it slides with a step size Δw=1 as time progresses.

[0025] Within each fixed window, iterate through any two non-identical state vectors within the window. and Calculate the Euclidean distance between the two. .

[0026] During the stable period of normal operation, extract the reference phase space matrix, such as the phase space constructed from the first 5000 data points. Iterate through all state vector pairs to calculate the Euclidean distance and obtain the overall average distance. If the average distance is calculated to be 3.5, multiply it by 10% to set the threshold radius ϵ=0.35.

[0027] In the moving window, a comparison is performed using a decision function: if... If the value is less than 0.35, the corresponding element in the state matrix is ​​assigned a value of 1, which is recorded as a recurrence state point; otherwise, it is assigned a value of 0.

[0028] After traversal, sum all recurrence state points with a value of 1 extracted in the current window, excluding self-matching cases (i.e., not calculating the distance at i=j on the diagonal). The total number of recurrence points is then calculated. =15000, and then divide it by the total number of possible ordered non-self-matching state points within the window, 249500, to calculate the local recurrence rate RR of the state trajectory corresponding to the window, RR=0.0601, or 6.01%. Based on this, the spatial clustering of the system state trajectory within the current local time period can be quantified.

[0029] S4. When the local recurrence rate is lower than the short-term mean and the absolute value of the difference between the local recurrence rate and the short-term mean exceeds the fluctuation tolerance, it is determined that the state of instability has entered the precursor state and subsequent steps are triggered.

[0030] In the controller's backend, a fixed-length exponential moving average filter is continuously run to calculate the short-term mean of the local recurrence rate in real time. The absolute value of the difference between the current local recurrence rate and the short-term mean is evaluated.

[0031] When the absolute value is greater than the fluctuation tolerance calculated based on three standard deviations of historical steady-state data, and the current local recurrence rate is lower than the short-term average, the interrupt function based on the event triggering mechanism is called to determine whether the fluid is undergoing a transition from laminar to turbulent flow or cavitation precursors, a Boolean trigger flag is generated, and the order reduction identification module is activated.

[0032] S5, subspace identification is performed on the data before triggering to obtain a higher-order model.

[0033] The main pump motor control command, generated within a specified window length prior to the trigger flag, is extracted and used as the control input data for system identification. The pressure sequence within the same time period is extracted and used as the output response data for system identification, constructing corresponding input and output data matrices. The flow rate and pressure sequences are used together for the aforementioned phase-space state monitoring and instability precursor triggering. A numerical subspace state-space system identification method is employed to construct block Hankel matrices corresponding to the input and output data, and these block Hankel matrices are divided into past and subsequent data blocks.

[0034] By performing singular value decomposition on the projection result matrix used for subspace identification, the system extended observability matrix or observable subspace is estimated, and the system state matrix, input matrix, output matrix and direct transfer matrix of the discrete-time high-order state-space model are obtained by using the least squares method.

[0035] As an optional implementation, the step of obtaining a higher-order model by subspace identification of the pre-trigger data includes: Extract the control command sequence of the main pump motor of a fixed length before the trigger moment, construct the input block Hankel matrix according to the time sliding sampling method, and divide it into past input data matrix and subsequent input data matrix; extract the pressure sequence of a fixed length before the trigger moment, construct the output block Hankel matrix according to the same time sliding sampling method, and divide it into past output data matrix and subsequent output data matrix; Based on the past input data matrix, past output data matrix, and subsequent input data matrix, the subsequent output data matrix is ​​projected to obtain a projection result matrix for subspace identification. Singular value decomposition is performed on the projection result matrix to obtain a left singular vector matrix, a right singular vector matrix, and a diagonal matrix with singular values ​​on the diagonal. The observability matrix of the system is reconstructed by retaining the first-order maximum singular value in the diagonal matrix and the corresponding left and right singular vectors. The transition matrix equations are then solved to obtain a high-order model that includes the state matrix, input control matrix, output observation matrix and feedforward matrix.

[0036] When the system is determined to be in an unstable precursor state, backtracking is performed to extract control and response data from a fixed period of 2 seconds prior to the trigger point, for example, data containing M=2000 sampling points. The input block Hankel matrix is ​​constructed using the main pump motor control command sequence u(t), with row block parameters i=20 and column numbers j=M-2i+1=1961; for the k-th column, the past input data matrix... It consists of i consecutive sampling points from u(k) to u(k+i-1), and the subsequent input data matrix It consists of i consecutive sampling points from u(k+i) to u(k+2i-1), where k slides sequentially from 1 to j, thus forming and Similarly, the past output data matrix is ​​constructed using the response pressure sequence y(t) according to the same sliding sampling rule. and subsequent output data matrix .

[0037] Past input data matrix Compared with past output data matrix Concatenate row by row to form a joint input / output data matrix. ; then input the data matrix As the control input influence terms to be eliminated, the subsequent output data matrix Orthogonalization is performed, and the processed subsequent output data matrix is ​​projected onto the past input-output joint data matrix. The row space is used to obtain the projection result matrix O for subspace identification.

[0038] Performing singular value decomposition on the projection matrix O yields a left singular vector matrix U, a right singular vector matrix V, and a singular value diagonal matrix S. Since the system's energy is typically concentrated in a few principal states, the singular values ​​along the diagonal of S exhibit a steep descent characteristic; for example, the first 25 singular values ​​are all on the order of 10. 1 The above decreases to 10 starting from the 26th. -3 A truncated singular value diagonal matrix is ​​constructed by extracting the 25 largest singular values ​​set beforehand. Extract the corresponding left singular vector matrix. This allows us to obtain the extended observability matrix. Singular value distribution as follows Figure 3 As shown.

[0039] use By utilizing the displacement-invariant property, the first i-1 blocks and the last i-1 blocks of the matrix are extracted, and the system state matrix is ​​obtained by solving using the least squares method. With output observation matrix Simultaneously, the input control matrix is ​​reconstructed by solving the residual equations. and feedforward matrix Thus, a 25th-order linear discrete high-order state-space model that can approximately characterize the dynamic features of this transient condition is obtained.

[0040] S6, calculate the Hankel singular values ​​and eigenvalues ​​of the model.

[0041] First, calculate all the eigenvalues ​​of the state matrix A in the high-order discrete state-space model to obtain the eigenvalue set. Calculate the magnitude of each eigenvalue one by one. The modes corresponding to eigenvalues ​​with amplitudes not less than 1 are marked as unstable modes, and the modes corresponding to eigenvalues ​​with amplitudes less than 1 are marked as stable modes. For paired complex conjugate eigenvalues, their corresponding modes are retained or divided as a group.

[0042] Construct the mode transformation matrix T obtained from eigenvectors or real Schur decomposition, and transform the original high-order model to the modal coordinate system to obtain... , , Based on the indices of unstable and stable modes, , and By dividing the system into blocks, an unstable subsystem is obtained. , , ) and stable subsystem ( , , In this case, all modes of the unstable subsystem are not truncated and are directly retained as the fixed part of the subsequent reduced-order model.

[0043] For a stable subsystem, since Since all eigenvalues ​​have amplitudes less than 1, the discrete Lyapunov equation can be solved to obtain the controllable Gauram matrix. Solving the discrete Lyapunov equations yields the observable Gryllium matrix. Then the matrix is ​​calculated. The eigenvalues ​​of the stable subsystems are calculated, and the square roots of the non-negative eigenvalues ​​are taken to obtain the Hankel singular value sequence. These Hankel singular values ​​are then sorted in descending order of numerical value for subsequent energy ranking and determination of the order of the reduced stable modes. The unstable subsystems then participate in the subsequent order-reduction model reconstruction together with the selected stable subsystems.

[0044] S7, determine the order q of the reduced model, wherein the order q of the reduced model is the smallest positive integer that satisfies the condition that the singular value cumulative energy is not lower than the threshold and retains all unstable modes.

[0045] Arrange the Hankel singular values ​​of the stable subsystem in descending order, and calculate the ratio of the sum of squares of the first few Hankel singular values ​​to the sum of squares of all Hankel singular values ​​to obtain the cumulative energy ratio. Extract unstable poles with a modulus greater than or equal to 1 from the eigenvalues ​​of the state matrix, and count the number of such unstable modes.

[0046] An iterative search algorithm is used to find the smallest integer q that satisfies the condition that the cumulative energy ratio is not lower than a preset energy threshold and contains all unstable poles. This integer is used as the order q of the reduced-order model.

[0047] As an optional implementation, determining the order q of the reduced-order model, wherein the order q is the smallest positive integer that satisfies the condition that the singular value cumulative energy is not lower than a threshold and retains all unstable modes, includes: The state matrix of the higher-order model is decomposed to obtain all eigenvalues. Eigenvalues ​​with an amplitude of not less than 1 are extracted as unstable poles, and the number of poles is recorded as the unstable mode number constant. The higher-order model is decomposed into an unstable subsystem and a stable subsystem, while retaining all modes of the unstable subsystem; Calculate the Hankel singular values ​​for the stable subsystem and sort them in descending order of their numerical values; Starting with the first singular value, calculate the ratio of the sum of squares of the progressively increasing preceding singular values ​​to the sum of squares of all singular values ​​in the stable subsystem; Find the smallest positive integer that ensures the ratio of the sum of squares is not lower than a preset energy threshold, and assign the sum of the smallest positive integer and the constant of the number of unstable modes as the order q of the reduced model.

[0048] For the 25th order high-order model obtained by subspace identification, its system state matrix A is extracted and eigenvalue decomposition is performed using the characteristic polynomial to obtain the complex magnitudes of the 25 discrete poles.

[0049] Assuming that a pole is found through testing. =1.05+0.1i and =1.05-0.1i, with an amplitude of |1.05±0.1i|≈1.0547>1, is marked as an unstable pole, and the number constant of unstable modes is statistically obtained as k=2. The original model is decoupled using a block diagonal transformation matrix, separating the unstable subsystem containing these two poles and the stable subsystem containing the remaining 23 poles. For the 23rd-order stable subsystem, the Hankel singular value sequence is extracted by solving the controllability and observability Lyapunov equations. Let i = 1, 2, ..., 23, and sort the values ​​in descending order from largest to smallest, such as... =120.5、 =85.3、 =45.1.

[0050] Entering the cumulative energy ratio traversal loop, the threshold is set to 0.90, or 90%. When successively adding accumulation terms, assuming the initial p=5 terms, the energy percentage is calculated. =0.932>0.90, while the proportion is only 0.865 when p=4. Therefore, it can be determined that the smallest positive integer that satisfies the truncation criterion is p=5.

[0051] A merging calculation of the reduced-order dimensions is performed, algebraically summing the p=5 order modes to be retained in the stable subsystem with the k=2 order modes of the fully retained unstable part, establishing the reduced-order model order q=7. This step, while eliminating redundant modes with low energy proportions, helps to preserve the divergent evolution trend and energy information of the dominant operating characteristics of the original system.

[0052] S8. Perform modal decomposition on the higher-order model and calculate the wavelet approximation subband and detail subband coefficients of each stable mode impulse response.

[0053] The state matrix of the high-order model is decomposed into eigenvalues ​​or real Schur decomposition to obtain the eigenvalues ​​and mode transformation matrices corresponding to each mode, and the model is transformed into a modal coordinate system. The amplitude of each modal eigenvalue is calculated one by one, and the modes with amplitudes less than 1 are identified as stable modes. For paired modes corresponding to complex conjugate eigenvalues, they are treated as the same mode group. For each stable mode or stable mode group, a unit discrete pulse signal with an amplitude of 1 and lasting for one sampling period is applied, and the time-domain impulse response sequence of the mode within a preset sampling length is obtained by step-by-step recursion according to the discrete state update equation. Orthogonal wavelet basis functions are selected, and multi-level discrete wavelet decomposition is performed on the time-domain impulse response sequence. The last layer approximate sub-band coefficients are extracted as low-frequency feature coefficients, and the detail sub-band coefficients of each layer are extracted as high-frequency feature coefficients.

[0054] S9 identifies slowly varying trend modes based on the approximate subband coefficient energy ratio, and identifies key modes based on the detailed subband coefficient energy ratio among the remaining stable modes that are not identified as slowly varying trend modes.

[0055] Calculate the sum of squares of the approximate subband coefficients and the sum of squares of the detail subband coefficients for each stable mode, and calculate the proportion of the approximate subband energy to the total wavelet subband energy of that mode, as the low-frequency energy proportion. Set a low-frequency proportion threshold, and filter stable modes that exceed this threshold according to the low-frequency energy proportion from high to low, and determine the slow-varying trend mode set by combining the total number of stable modes to be filtered.

[0056] In the remaining stable modes, the proportion of detail subband energy to the total wavelet subband energy of the mode is calculated as the high-frequency energy proportion. The modes are sorted in descending order of high-frequency energy proportion, and the top few modes are extracted and saved to the key mode set. Redundant high-frequency modes with low remaining energy proportions are removed.

[0057] As an optional implementation, the step of identifying slowly varying trend modes based on the approximate subband coefficient energy ratio, and identifying key modes based on the detailed subband coefficient energy ratio among the remaining stable modes not identified as slowly varying trend modes, includes: Unit impulse excitation is applied to each stable mode of the higher-order model to obtain the impulse response time series corresponding to each stable mode; Multi-resolution orthogonal wavelet decomposition is performed on the impulse response time series of each stable mode to extract the detailed subband coefficients of the corresponding high frequency and the approximate subband coefficients of the corresponding low frequency. Calculate the approximate subband energy of each stable mode and its proportion of the total wavelet subband energy to obtain the low-frequency energy proportion; among the stable modes whose low-frequency energy proportion is greater than the set proportion threshold, select modes whose proportion is no more than the total number of stable modes to be screened minus 1, and confirm them as slow-changing trend modes, wherein the total number of stable modes to be screened is the order q of the reduced model minus the total number of unstable modes. Calculate the detail subband energy of the remaining stable modes and its proportion of the total wavelet subband energy to obtain the high-frequency energy ratio; select a target number of modes according to the high-frequency energy ratio from high to low to confirm them as key modes, wherein the target number is the total number of stable modes to be screened minus the total number of confirmed slow-changing trend modes.

[0058] For the 23 identified stable modes, ideal discrete pulse signals with an amplitude of 1 and a duration of one sampling period are applied to the input terminal. Under interference-free conditions, the pulse response output discrete sequence with a length of L=1024 points is obtained by solving the problem. m=1,2,...,23. Perform a five-level discrete wavelet transform on each response sequence to separate the fifth-level low-frequency approximate subband coefficient sequence cA5 in the 0-15.6Hz frequency band, and the first to fifth level detail subband coefficient sequence set [cD1,cD2,...,cD5] covering the broadband high-frequency fluctuation characteristics.

[0059] Calculate the proportion of the approximate subband energy of a single mode to the total wavelet subband energy of that mode, and set a threshold for the proportion of low-frequency energy. =0.60. If the low-frequency energy proportions of certain modes are 0.72, 0.68, and 0.63 respectively, then they are identified as slow-changing trend modes. At this time, the boundary conditions are checked: since the total number of stable modes to be retained is p=5, the upper limit is 4. The current number of 3 meets the rule, so these 3 modes are determined to be slow-changing trend modes of the dominant macroscopic fluctuations.

[0060] Extract all level detail subband coefficients from the remaining 20 stable modes that were not identified as slowly varying trend modes, calculate the proportion of detail subband energy to the total wavelet subband energy of that mode, and sort them according to the proportion of high-frequency energy. For example, the two highest high-frequency energy proportions are 0.55 and 0.49, respectively. Based on the formula, the target number to be added is 2. Accordingly, the two stable modes with the highest proportions are selected and defined as modes responding to high-frequency disturbances.

[0061] The above-mentioned operation, which combines wavelet band isolation and modal feature localization, enables the preserved 7th-order system structure to take into account unstable poles, low-frequency trend features, and high-frequency disturbance response features.

[0062] S10, combine the slow-changing trend mode, the key mode and all unstable modes to obtain a q-order reduced-order model.

[0063] The selected slow-changing trend modes, key modes, and all previously retained unstable modes are used as a retained mode set. In the modal coordinate system, the state sub-matrix, input coupling sub-matrix, output coupling sub-matrix, and direct transfer matrix corresponding to the retained mode set are extracted and reorganized according to a diagonal block structure. An inverse coordinate transformation is performed based on the transformation matrix corresponding to the retained modes to obtain a q-order reduced state-space model. This q-order reduced state-space model includes a reduced state matrix, a reduced input matrix, a reduced output matrix, and a reduced direct transfer matrix, used to predict the state at the next time step according to the discrete state update equation and to generate a pressure prediction output according to the output equation.

[0064] S11, the q-order reduced model is used for model predictive control to generate output prediction, and decomposed into low-frequency trend components and high-frequency disturbance components.

[0065] The q-order reduced model is configured as an internal prediction model, and the prediction and control time domains are defined. A quadratic objective function containing tracking error penalties and control increment penalties is constructed. The sequential quadratic programming solver is invoked to continuously calculate the output prediction sequence in the future time domain.

[0066] The predicted output sequence is segmented by calling a low-pass filter, a sliding smoothing filter, or a frequency decomposition algorithm based on the predicted sequence, separating the smooth low-frequency trend component and the high-frequency disturbance component containing transient fluctuations.

[0067] S12, calculate the main pump reference control sequence based on the low-frequency trend component.

[0068] The separated low-frequency trend component is input into the model predictive control optimization framework. Combined with the reference trajectory, control increment penalty term and physical amplitude constraint, a reference control sequence for adjusting the speed or displacement of the main pump motor is generated and sent to the main pump driver through the fieldbus interface to perform steady-state tracking control.

[0069] As an optional implementation, the step of calculating the main pump reference control sequence based on the low-frequency trend component includes: Establish a quadratic cost function that includes the sum of a weight function containing the square of the first-order difference term of the control quantity and a weight function containing the square of the low-frequency trend prediction error term. The low-frequency trend components of each moment within the predicted line of sight in the future are extracted to form a trend prediction output column vector. The test target reference data of the pumping system within the same time period are extracted to form a reference trajectory column vector. The trend error column vector is obtained by subtracting the trend prediction output column vector from the reference trajectory column vector. Based on the system's dynamic constraints and physical amplitude hard constraints, a sequential quadratic programming solver is used to minimize the quadratic cost function, solve for the main pump motor reference control command at the first future moment, and output it for execution.

[0070] Based on the evolution path of the baseline state variables obtained by real-time deduction of the recombined reduced-order model in the model predictive controller, the control time domain length is set. =5 and predicted sight distance parameters =15. Establish a quadratic cost function. In this setting, the directional correction parameter is assigned as the low-frequency error weight constant matrix. =diag([10,...,10]), and the damping stabilization penalty coefficient. =diag([0.5,...,0.5]). The first summation term iterates through the trend error within the prediction line of sight, and the second summation term iterates through the control increment within the control time domain. Both Δu and Δu are normalized variables.

[0071] A 15-dimensional prediction output column vector is constructed using the low-frequency trend component sequences corresponding to the time window obtained from the above decomposition. Simultaneously, a target reference column vector is constructed using the ideal pressure ramp-up planning data of the pumping equipment within the same step size under this test scenario. Calculate the error residual column vector .

[0072] Within this constrained space, hard limitations are applied using system physical barriers, such as setting the main pump motor speed output range as follows: =800rpm and =2200rpm, and set a hard limit constraint on the inverter speed for its first-order differential control increment, i.e. -100rpm≤Δu≤100rpm.

[0073] This cost function, along with the inequalities and equality boundary conditions, is loaded into the system's built-in sequential quadratic programming solver engine. Partial derivatives are iteratively calculated and the Hessian matrix is ​​updated until the function converges to find a local minimum in the global control space. After convergence, the obtained value is extracted... The first value of the reference command column vector composed of discrete elements, for example, if the speed increment obtained at this time is Δu=45rpm, then the current basic setting value is the speed at the previous moment plus this increment, which is 1545rpm. This value is latched and output as the formal main pump motor reference control command for controlling the cycle update step, thereby improving the fit between the system's basic pressure establishment process and the macroscopic target.

[0074] S13 uses the high-frequency disturbance component as a feedforward signal to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control.

[0075] High-frequency disturbance components are extracted and converted into valve compensation control commands using a dynamic feedforward compensator. These commands are then sent to a faster-responding bypass fine-tuning electromagnetic proportional valve via a digital-to-analog converter or fieldbus to execute dynamic compensation actions.

[0076] This allows for the maintenance of stable macroscopic flow rate during pumping, while simultaneously reducing high-frequency pulsation amplitude through the coordinated action of the main pump and the fine-tuning valve's dual execution channels.

[0077] As an optional implementation, the step of using high-frequency disturbance components as feedforward signals to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control includes: A dynamic feedforward compensator calculation model based on the discrete time domain is established, which has independent calculation channels for proportional, integral, and differential calculations. The first control quantity is generated by multiplying the input value of the high-frequency disturbance component generated in real time by a preset proportional constant. The second control quantity is generated by discretely summing the high-frequency disturbance components within a fixed window and multiplying them by a preset integral constant. The high-frequency disturbance component of the current sampling step is subtracted from the high-frequency disturbance component of the previous sampling step by a first-order backward difference operation, and then multiplied by a preset differential constant to generate the third control quantity. The valve compensation control quantity is obtained by algebraically summing the first control quantity, the second control quantity, and the third control quantity, generating a dynamic compensation control sequence and sending it to the bypass fine-tuning valve servo drive.

[0078] The high-frequency disturbance component separated from the model predictive control output, i.e., the response residual signal. As a feedforward excitation source, a feedforward dynamic compensation model is established using a discrete PID structured operation function block pre-integrated into the industrial controller. The preset proportional constant, integral constant, and derivative constant are all calibrated dimensionless conversion gains, which convert the first control quantity, the second control quantity, and the third control quantity into a unified valve compensation control quantity.

[0079] When implementing the operation in each processing channel, the PID tuning parameter group is first set using empirical parameter tuning methods or relay self-tuning methods, such as configuring the proportional constant. =2.5, Integral constant =0.8, differential constant =0.15.

[0080] In the proportional calculation channel, the real-time high-frequency disturbance value at the k-th discrete time point is obtained as follows: =0.4MPa, the first control value of the proportional channel is obtained by multiplying. =1.0, because Valve opening compensation dimension calibration has been completed. This indicates the valve compensation control value. Simultaneously, a historical data cache queue is started, and the retrieval time is set. The high-frequency component set within a cumulative sliding window of 50 steps is discretely accumulated and summed. Assuming the sum is 1.2, it is multiplied by a pre-set integration constant to reduce long-term steady-state error, and the second control variable is obtained. =0.96, because It has been calibrated according to the sampling period and valve opening compensation dimensions. This represents the valve compensation control quantity. The high-frequency spike ripple rate of pressure is detected in the differential independent channel, and the current value is taken. =0.4MPa and retrieve the previous clock value from the k-1 buffer register. =0.2MPa, perform a first-order backward differential operation, which yields 0.2, then multiply it with the differential gain to generate the third control quantity. =0.03. Because It has been calibrated according to the sampling period and valve opening compensation dimensions. It also indicates the valve compensation control quantity.

[0081] The valve compensation control quantities obtained from the three channels are linearly superimposed to obtain... =1.99. This value, after being converted by limiting and calibrating the ratio, forms the opening compensation command of the bypass fine-tuning valve, and is encapsulated into a high-speed dynamic compensation control sequence. The control frequency is, for example, 100Hz, and is sent to the high-speed electromagnetic servo driver of the bypass fine-tuning valve based on the EtherCAT communication bus to reduce the high-frequency oscillation amplitude of the system in a parallel intervention manner.

[0082] In a second embodiment, the present invention also proposes a cooperative control system for pumpability testing, comprising: The calculation module is used to collect control commands from the main pump motor and pressure and flow signals from the pipeline in real time; reconstruct the phase space based on the pressure and flow signals; calculate the local recurrence rate of the state trajectory; when the local recurrence rate is lower than the short-term mean and the absolute value of the difference between the local recurrence rate and the short-term mean exceeds the fluctuation tolerance, it is determined that the system has entered an unstable precursor state and triggers subsequent steps. The combination module is used to perform subspace identification on the pre-trigger data to obtain a higher-order model; calculate the Hankel singular values ​​and eigenvalues ​​of the model; determine the order q of the reduced-order model, where the reduced-order model order q is the smallest positive integer that satisfies the condition that the accumulated singular value energy is not lower than a threshold and retains all unstable modes; perform mode decomposition on the higher-order model, calculate the wavelet approximate subband and detail subband coefficients of the impulse response of each stable mode; identify the slow-varying trend mode based on the energy proportion of the approximate subband coefficients, and identify the critical mode among the remaining stable modes not identified as slow-varying trend modes based on the energy proportion of the detail subband coefficients; and combine the slow-varying trend mode, the critical mode, and all unstable modes to obtain the q-order reduced-order model. The control module is used to generate output predictions by using the q-order reduced model for model predictive control, and decomposes it into low-frequency trend components and high-frequency disturbance components; calculates the main pump reference control sequence based on the low-frequency trend components; and uses the high-frequency disturbance components as feedforward signals to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control.

[0083] As an optional implementation, the reconstruction of the phase space based on the pressure and flow signals includes: The collected pressure and flow signals are subtracted from their corresponding time-series average values, and the subtraction results are divided by the corresponding standard deviations to obtain the normalized pressure sequence and the normalized flow sequence. The delay time of the phase space is determined by calculating the first minimum point using the mutual information method, and the embedding dimension of the phase space is determined by the false nearest neighbor method when the proportion of false neighbors drops below the preset neighbor threshold. Using the determined delay time and embedding dimension, the normalized pressure sequence and normalized flow sequence are reconstructed in a multivariate phase space. The resulting state vectors each contain a time series phase space matrix of pressure and flow data corresponding to the number of dimensions specified by the embedding dimension.

[0084] As an optional implementation, the calculation of the local recurrence rate of the state trajectory includes: Extract the phase space moving window where the current time step is located, traverse any two non-identical state vectors in the phase space moving window, and calculate the Euclidean distance between the two state vectors. When the Euclidean distance is less than a preset threshold radius, a relapse state point is recorded. The total number of recurring state points within the phase space moving window is counted, and the total number of recurring state points is divided by the total number of calculable paired state point combinations within the phase space moving window to obtain the local recurrence rate of the state trajectory.

[0085] As an optional implementation, the step of obtaining a higher-order model by subspace identification of the pre-trigger data includes: Extract the control command sequence of the main pump motor of a fixed length before the trigger moment, construct the input block Hankel matrix according to the time sliding sampling method, and divide it into past input data matrix and subsequent input data matrix; extract the pressure sequence of a fixed length before the trigger moment, construct the output block Hankel matrix according to the same time sliding sampling method, and divide it into past output data matrix and subsequent output data matrix; Based on the past input data matrix, past output data matrix, and subsequent input data matrix, the subsequent output data matrix is ​​projected to obtain a projection result matrix for subspace identification. Singular value decomposition is performed on the projection result matrix to obtain a left singular vector matrix, a right singular vector matrix, and a diagonal matrix with singular values ​​on the diagonal. The observability matrix of the system is reconstructed by retaining the first-order maximum singular value in the diagonal matrix and the corresponding left and right singular vectors. The transition matrix equations are then solved to obtain a high-order model that includes the state matrix, input control matrix, output observation matrix and feedforward matrix.

[0086] As an optional implementation, determining the order q of the reduced-order model, wherein the order q is the smallest positive integer that satisfies the condition that the singular value cumulative energy is not lower than a threshold and retains all unstable modes, includes: The state matrix of the higher-order model is decomposed to obtain all eigenvalues. Eigenvalues ​​with an amplitude of not less than 1 are extracted as unstable poles, and the number of poles is recorded as the unstable mode number constant. The higher-order model is decomposed into an unstable subsystem and a stable subsystem, while retaining all modes of the unstable subsystem; Calculate the Hankel singular values ​​for the stable subsystem and sort them in descending order of their numerical values; Starting with the first singular value, calculate the ratio of the sum of squares of the progressively increasing preceding singular values ​​to the sum of squares of all singular values ​​in the stable subsystem; Find the smallest positive integer that ensures the ratio of the sum of squares is not lower than a preset energy threshold, and assign the sum of the smallest positive integer and the constant of the number of unstable modes as the order q of the reduced model.

[0087] As an optional implementation, the step of identifying slowly varying trend modes based on the approximate subband coefficient energy ratio, and identifying key modes based on the detailed subband coefficient energy ratio among the remaining stable modes not identified as slowly varying trend modes, includes: Unit impulse excitation is applied to each stable mode of the higher-order model to obtain the impulse response time series corresponding to each stable mode; Multi-resolution orthogonal wavelet decomposition is performed on the impulse response time series of each stable mode to extract the detailed subband coefficients of the corresponding high frequency and the approximate subband coefficients of the corresponding low frequency. Calculate the approximate subband energy of each stable mode and its proportion of the total wavelet subband energy to obtain the low-frequency energy proportion; among the stable modes whose low-frequency energy proportion is greater than the set proportion threshold, select modes whose proportion is no more than the total number of stable modes to be screened minus 1, and confirm them as slow-changing trend modes, wherein the total number of stable modes to be screened is the order q of the reduced model minus the total number of unstable modes. Calculate the detail subband energy of the remaining stable modes and its proportion of the total wavelet subband energy to obtain the high-frequency energy ratio; select a target number of modes according to the high-frequency energy ratio from high to low to confirm them as key modes, wherein the target number is the total number of stable modes to be screened minus the total number of confirmed slow-changing trend modes.

[0088] As an optional implementation, the step of calculating the main pump reference control sequence based on the low-frequency trend component includes: Establish a quadratic cost function that includes the sum of a weight function containing the square of the first-order difference term of the control quantity and a weight function containing the square of the low-frequency trend prediction error term. The low-frequency trend components of each moment within the predicted line of sight in the future are extracted to form a trend prediction output column vector. The test target reference data of the pumping system within the same time period are extracted to form a reference trajectory column vector. The trend error column vector is obtained by subtracting the trend prediction output column vector from the reference trajectory column vector. Based on the system's dynamic constraints and physical amplitude hard constraints, a sequential quadratic programming solver is used to minimize the quadratic cost function, solve for the main pump motor reference control command at the first future moment, and output it for execution.

[0089] As an optional implementation, the step of using high-frequency disturbance components as feedforward signals to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control includes: A dynamic feedforward compensator calculation model based on the discrete time domain is established, which has independent calculation channels for proportional, integral, and differential calculations. The first control quantity is generated by multiplying the input value of the high-frequency disturbance component generated in real time by a preset proportional constant. The second control quantity is generated by discretely summing the high-frequency disturbance components within a fixed window and multiplying them by a preset integral constant. The high-frequency disturbance component of the current sampling step is subtracted from the high-frequency disturbance component of the previous sampling step by a first-order backward difference operation, and then multiplied by a preset differential constant to generate the third control quantity. The valve compensation control quantity is obtained by algebraically summing the first control quantity, the second control quantity, and the third control quantity, generating a dynamic compensation control sequence and sending it to the bypass fine-tuning valve servo drive.

[0090] In this specification, relational terms such as "first" and "second" are used merely to distinguish one entity or operation from another, and do not necessarily require or imply any such actual relationship or order between these entities or operations. Furthermore, the terms "comprising," "including," or any other variations thereof are intended to cover non-exclusive inclusion, such that a process, method, article, or apparatus that comprises a list of elements includes not only those elements but also other elements not expressly listed, or elements inherent to such a process, method, article, or apparatus. Unless otherwise limited, an element defined by the phrase "comprising one..." does not exclude the presence of other identical elements in the process, method, article, or apparatus that includes said element. In this document, "a," "an," "the," "the," and "its" may also include plural forms unless the context clearly indicates otherwise. "Multiple" refers to at least two, such as 2, 3, 5, or 8, etc. "And / or" includes any and all combinations of the associated listed items.

[0091] The various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. The various embodiments can be combined as needed, and the same or similar parts can be referred to each other.

[0092] The above description of the disclosed embodiments enables those skilled in the art to make or use this application. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of this application. Therefore, this application is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.

Claims

1. A collaborative control method for pumpability testing, characterized in that, include: Real-time acquisition of control commands from the main pump motor and pressure and flow signals from the pipeline; reconstruction of the phase space based on the pressure and flow signals; Calculate the local recurrence rate of the state trajectory; When the local recurrence rate is lower than the short-term mean and the absolute value of the difference from the short-term mean exceeds the fluctuation tolerance, it is determined that the state of instability precursor is entered and subsequent steps are triggered. Subspace identification is performed on the data before triggering to obtain a higher-order model; the Hankel singularity value and eigenvalue of the model are calculated; The order q of the reduced-order model is determined, which is the smallest positive integer that satisfies the condition that the singular value cumulative energy is not lower than a threshold and retains all unstable modes. Modal decomposition is performed on the higher-order model, and the wavelet approximate subband and detail subband coefficients of the impulse response of each stable mode are calculated. Slow-trend modes are identified based on the energy proportion of the approximate subband coefficients, and key modes are identified based on the energy proportion of the detail subband coefficients among the remaining stable modes that are not identified as slow-trend modes. The slow-changing trend mode, the key mode, and all unstable modes are combined to obtain a q-order reduced model; The q-order reduced model is used to generate output predictions for model predictive control, and is decomposed into low-frequency trend components and high-frequency disturbance components. The main pump reference control sequence is calculated based on the low-frequency trend component; the dynamic compensation control sequence of the bypass fine-tuning valve is calculated using the high-frequency disturbance component as a feedforward signal for coordinated control.

2. The method according to claim 1, characterized in that, The reconstruction of the phase space based on the pressure and flow signals includes: The collected pressure and flow signals are subtracted from their corresponding time-series average values, and the subtraction results are divided by the corresponding standard deviations to obtain the normalized pressure sequence and the normalized flow sequence. The delay time of the phase space is determined by calculating the first minimum point using the mutual information method, and the embedding dimension of the phase space is determined by the false nearest neighbor method when the proportion of false neighbors drops below the preset neighbor threshold. Using the determined delay time and embedding dimension, the normalized pressure sequence and normalized flow sequence are reconstructed in a multivariate phase space, and the resulting state vectors are combined to generate a time series phase space matrix containing a number of dimensions of pressure and flow data corresponding to the embedding dimension.

3. The method according to claim 1, characterized in that, The calculation of the local recurrence rate of the state trajectory includes: Extract the phase space moving window where the current time step is located, traverse any two non-identical state vectors in the phase space moving window, and calculate the Euclidean distance between the two state vectors. When the Euclidean distance is less than a preset threshold radius, a relapse state point is recorded. The total number of recurring state points within the phase space moving window is counted, and the total number of recurring state points is divided by the total number of calculable paired state point combinations within the phase space moving window to obtain the local recurrence rate of the state trajectory.

4. The method according to claim 2, characterized in that, The step of obtaining a higher-order model by subspace identification of the data before triggering includes: Extract the control command sequence of the main pump motor of a fixed length before the trigger moment, construct the input block Hankel matrix according to the time sliding sampling method, and divide it into past input data matrix and subsequent input data matrix; extract the pressure sequence of a fixed length before the trigger moment, construct the output block Hankel matrix according to the same time sliding sampling method, and divide it into past output data matrix and subsequent output data matrix; Based on the past input data matrix, past output data matrix, and subsequent input data matrix, the subsequent output data matrix is ​​projected to obtain a projection result matrix for subspace identification. Singular value decomposition is performed on the projection result matrix to obtain a left singular vector matrix, a right singular vector matrix, and a diagonal matrix with singular values ​​on the diagonal. The observability matrix of the system is reconstructed by retaining the first-order maximum singular value in the diagonal matrix and the corresponding left and right singular vectors. The transition matrix equations are then solved to obtain a high-order model that includes the state matrix, input control matrix, output observation matrix and feedforward matrix.

5. The method according to claim 1, characterized in that, The determination of the order q of the reduced-order model, wherein the order q of the reduced-order model is the smallest positive integer that satisfies the condition that the singular value cumulative energy is not lower than a threshold and retains all unstable modes, includes: The state matrix of the higher-order model is decomposed to obtain all eigenvalues. Eigenvalues ​​with an amplitude of not less than 1 are extracted as unstable poles, and the number of poles is recorded as the unstable mode number constant. The higher-order model is decomposed into an unstable subsystem and a stable subsystem, while retaining all modes of the unstable subsystem; Calculate the Hankel singular values ​​for the stable subsystem and sort them in descending order of their numerical values; Starting with the first singular value, calculate the ratio of the sum of squares of the progressively increasing preceding singular values ​​to the sum of squares of all singular values ​​in the stable subsystem; Find the smallest positive integer that ensures the ratio of the sum of squares is not lower than a preset energy threshold, and assign the sum of the smallest positive integer and the constant of the number of unstable modes as the order q of the reduced model.

6. The method according to claim 1, characterized in that, The step of identifying slowly varying trend modes based on the approximate subband coefficient energy ratio, and identifying key modes based on the detailed subband coefficient energy ratio among the remaining stable modes not identified as slowly varying trend modes, includes: Unit impulse excitation is applied to each stable mode of the higher-order model to obtain the impulse response time series corresponding to each stable mode; Multi-resolution orthogonal wavelet decomposition is performed on the impulse response time series of each stable mode to extract the detailed subband coefficients of the corresponding high frequency and the approximate subband coefficients of the corresponding low frequency. Calculate the approximate subband energy of each stable mode and its proportion of the total wavelet subband energy to obtain the low-frequency energy proportion; among the stable modes whose low-frequency energy proportion is greater than the set proportion threshold, select modes whose proportion is no more than the total number of stable modes to be screened minus 1, and confirm them as slow-changing trend modes, wherein the total number of stable modes to be screened is the order q of the reduced model minus the total number of unstable modes. Calculate the detail subband energy of the remaining stable modes and its proportion of the total wavelet subband energy to obtain the high-frequency energy ratio; select a target number of modes according to the high-frequency energy ratio from high to low to confirm them as key modes, wherein the target number is the total number of stable modes to be screened minus the total number of confirmed slow-changing trend modes.

7. The method according to claim 1, characterized in that, The calculation of the main pump reference control sequence based on the low-frequency trend component includes: Establish a quadratic cost function that includes the sum of a weight function containing the square of the first-order difference term of the control quantity and a weight function containing the square of the low-frequency trend prediction error term. The low-frequency trend components of each moment within the predicted line of sight in the future are extracted to form a trend prediction output column vector. The test target reference data of the pumping system within the same time period are extracted to form a reference trajectory column vector. The trend error column vector is obtained by subtracting the trend prediction output column vector from the reference trajectory column vector. Based on the system's dynamic constraints and physical amplitude hard constraints, a sequential quadratic programming solver is used to minimize the quadratic cost function, solve for the main pump motor reference control command at the first future moment, and output it for execution.

8. The method according to claim 1, characterized in that, The method of using high-frequency disturbance components as feedforward signals to calculate the dynamic compensation control sequence of the bypass fine-tuning valve for coordinated control includes: A dynamic feedforward compensator calculation model based on the discrete time domain is established, which has independent calculation channels for proportional, integral, and differential calculations. The first control quantity is generated by multiplying the input value of the high-frequency disturbance component generated in real time by a preset proportional constant. The second control quantity is generated by discretely summing the high-frequency disturbance components within a fixed window and multiplying them by a preset integral constant. The high-frequency disturbance component of the current sampling step is subtracted from the high-frequency disturbance component of the previous sampling step by a first-order backward difference operation, and then multiplied by a preset differential constant to generate the third control quantity. The valve compensation control quantity is obtained by algebraically summing the first control quantity, the second control quantity, and the third control quantity, generating a dynamic compensation control sequence and sending it to the bypass fine-tuning valve servo drive.

9. A collaborative control system for pumpability testing, characterized in that, include: The calculation module is used to acquire control commands from the main pump motor and pressure and flow signals from the pipeline in real time; and to reconstruct the phase space based on the pressure and flow signals. Calculate the local recurrence rate of the state trajectory; when the local recurrence rate is lower than the short-term mean and the absolute value of the difference between the local recurrence rate and the short-term mean exceeds the fluctuation tolerance, it is determined that the state has entered an unstable precursor state and subsequent steps are triggered. The combination module is used to perform subspace identification on the pre-trigger data to obtain a higher-order model; and to calculate the Hankel singular value and eigenvalue of the model. The order q of the reduced-order model is determined, which is the smallest positive integer that satisfies the condition that the singular value cumulative energy is not lower than a threshold and retains all unstable modes. Modal decomposition is performed on the higher-order model, and the wavelet approximate subband and detail subband coefficients of the impulse response of each stable mode are calculated. Slow-trend modes are identified based on the energy proportion of the approximate subband coefficients, and key modes are identified based on the energy proportion of the detail subband coefficients among the remaining stable modes that are not identified as slow-trend modes. The slow-changing trend mode, the key mode, and all unstable modes are combined to obtain a q-order reduced model; The control module is used to use the q-order reduced model for model predictive control to generate output predictions and decompose them into low-frequency trend components and high-frequency disturbance components. The main pump reference control sequence is calculated based on the low-frequency trend component; the dynamic compensation control sequence of the bypass fine-tuning valve is calculated using the high-frequency disturbance component as a feedforward signal for coordinated control.

10. A cooperative control device for pumpability testing, comprising a memory, a processor, and a computer program stored in the memory and running on the processor, characterized in that, When the processor executes the computer program, it implements the method as described in any one of claims 1 to 8.

Images