Adaptive control method and system for a filling installation

Abstract

The application discloses a kind of self-adapting control method and system of filling equipment, and the accurate compensation of multiple temperature point influence is realized by the hybrid compensation model combining physical model and neural network, the filling error is separated into random error, periodic error and trend drift three components and is carried out predictive compensation using empirical mode decomposition technique, adaptive setting of controller parameters is realized based on online system identification and multi-objective optimization, system reliable operation is guaranteed through multi-level fault detection and fault-tolerant control strategy, and finally global collaborative optimization is realized through hierarchical optimization strategy.The application solves the technical problems such as simple temperature compensation mechanism of traditional filling control system, lack of adaptive ability of control parameters, limited fault diagnosis capability, insufficient multi-objective collaborative optimization, etc., improves the filling precision from ±2% to ±0.3%, and significantly improves the precision, stability and adaptability of filling equipment.

Description

Technical Field

[0001] This invention relates to the field of industrial automation control technology, specifically to precision control technology for filling equipment, and particularly to an adaptive control method and system for filling equipment. Background Technology

[0002] Industrial filling technology, as an important component of modern manufacturing, is widely used in food, pharmaceutical, and chemical industries. The precision, stability, and adaptability of filling equipment directly affect product quality and production efficiency, especially the precise control required under conditions of varying temperatures, product batches, and environmental interference.

[0003] Currently, common filling technologies mainly fall into two categories: volumetric filling and flow metering filling. Volumetric filling controls the quantity of filling by pre-setting a volume space. It has a simple structure, but its accuracy is limited by the precision of mechanical components. Flow metering filling, on the other hand, achieves quantitative control by monitoring flow parameters in real time and controlling the filling time. Theoretically, it can achieve higher accuracy, but it is easily affected by external factors such as temperature and pressure fluctuations.

[0004] Existing high-precision filling control systems typically employ PID control algorithms combined with simple temperature compensation mechanisms. The control system monitors filling process parameters using pressure sensors and flow meters, adjusting valve opening and filling time based on a preset control model. This type of control system can provide a certain level of filling accuracy under stable conditions, and some advanced systems also employ simple feedback correction mechanisms to address batch-to-batch variations.

[0005] However, these traditional systems have obvious shortcomings:

[0006] First, the temperature compensation mechanism is too simplistic. Traditional systems typically monitor only a single temperature point (such as liquid temperature) and use a linear compensation model, which cannot handle the complex coupling effects of multiple temperature points. In actual filling processes, liquid temperature, gas temperature, sensor body temperature, and ambient temperature all affect filling accuracy through different mechanisms. These effects are interdependent and exhibit nonlinear characteristics, making simple single-point linear compensation insufficient to meet high-precision requirements.

[0007] Second, the control parameters lack adaptive capability. Traditional systems typically use PID parameters that are manually tuned during the commissioning phase and then fixed, making it difficult to cope with dynamic factors such as changes in the characteristics of different batches of products, equipment aging, and changes in environmental conditions. When the physical properties of the liquid, such as viscosity and density, change, the fixed control parameters cannot maintain optimal performance, leading to fluctuations in filling accuracy.

[0008] Third, the fault diagnosis capability is limited. Traditional systems mainly rely on simple threshold alarms, and cannot adjust the control strategy in time when abnormalities occur. They lack the ability to identify gradual faults such as sensor drift and actuator performance degradation in the early stages, and are often only discovered when the fault seriously affects product quality.

[0009] Fourth, the system optimization objective is singular. Traditional control systems mainly focus on the single indicator of filling accuracy, failing to comprehensively consider multiple performance objectives such as response speed, control smoothness, and energy consumption. They cannot achieve multi-objective collaborative optimization, and may sacrifice other important performance in the pursuit of ultimate accuracy.

[0010] These issues lead to significant fluctuations in filling accuracy during actual production, typically achieving only ±2% accuracy. Especially in environments with large temperature variations and frequent product batch changes, filling accuracy and stability are difficult to guarantee, affecting product quality consistency and production efficiency. Summary of the Invention

[0011] The purpose of this invention is to provide an adaptive control method and apparatus for filling equipment, so as to solve the technical problems in the prior art, such as simple temperature compensation mechanism, lack of adaptive control parameters, limited fault diagnosis capability, and single system optimization objective leading to large fluctuations in filling accuracy.

[0012] To achieve the above objectives, the present invention adopts the following technical solution:

[0013] An adaptive control method for filling equipment includes:

[0014] Real-time temperature data of liquid temperature, gas temperature, sensor body temperature and ambient temperature of filling equipment are collected. The influence coefficient of each temperature change on the filling volume is calculated by a hybrid compensation model that combines physical model and neural network, and the corrected filling control parameters are output.

[0015] Filling is performed according to the modified filling control parameters and filling error data is collected. Empirical mode decomposition is performed on the filling error data to separate the filling error into random error components, periodic error components and trend drift components. Error trend prediction is performed based on each error component and time compensation value is calculated.

[0016] Based on the corrected filling control parameters and the time compensation value, the flow coefficient, time constant and valve response delay parameters of the filling equipment are obtained, and multi-objective optimization calculations are performed on the flow coefficient, time constant and valve response delay parameters to output the optimal controller parameters;

[0017] The filling equipment is controlled to operate according to the optimal controller parameters. The operating status data of the filling equipment is monitored and analyzed in real time through multi-level fault detection to obtain fault diagnosis results. Based on the fault diagnosis results, the corresponding fault-tolerant control strategy is activated.

[0018] Based on the fault-tolerant control strategy and the operating status data of the filling equipment, the control parameters of the filling equipment are collaboratively optimized and the filling is controlled through a hierarchical optimization strategy.

[0019] Furthermore, the hybrid compensation model, which combines a physical model with a neural network, calculates the influence coefficient of temperature changes on the filling volume and outputs corrected filling control parameters, including:

[0020] The real-time temperature data is filtered, and the filtered temperature data is input into the corresponding component of the physical model to calculate the individual influence coefficient of each temperature point on the filling volume.

[0021] The real-time temperature data is input into the neural network for nonlinear compensation calculation. The neural network takes the liquid temperature, gas temperature, sensor body temperature and ambient temperature as inputs and outputs nonlinear coupling compensation coefficients.

[0022] Based on the individual influence coefficients of each temperature point and the nonlinear coupling compensation coefficient, a comprehensive compensation coefficient is generated by fusion calculation according to the adaptive weights of each temperature point. The filling control time is then corrected based on the comprehensive compensation coefficient, and the corrected filling control parameters are output.

[0023] Further, the step of performing empirical mode decomposition on the filling error data to separate the filling error into random error components, periodic error components, and trend drift components, and then predicting the error trend and calculating the time compensation value based on each error component, includes:

[0024] The filling error data is subjected to iterative screening. The intrinsic mode function is extracted by identifying extreme points and the mean of the envelope. The filling error data is decomposed into multiple intrinsic mode function components and residual components. According to the frequency characteristics of each component, the intrinsic mode function components and residual components are classified as the random error component, the periodic error component, and the trend drift component.

[0025] The trend drift component is calculated using a prediction algorithm that combines linear extrapolation and curvature correction. The periodic error component is analyzed in the frequency domain to extract the main frequency and amplitude parameters and calculate the periodic prediction value. The trend prediction value and the periodic prediction value are integrated to determine the comprehensive error prediction value.

[0026] Based on the comprehensive error prediction value, the filling control adjustment amount is calculated by combining the feedback compensation component and the feedforward compensation component, and the filling control adjustment amount is converted into the time compensation value.

[0027] Furthermore, the process of acquiring the flow coefficient, time constant, and valve response delay parameters of the filling equipment, performing multi-objective optimization calculations on the flow coefficient, time constant, and valve response delay parameters, and outputting optimal controller parameters includes:

[0028] The current working state of the filling system is determined based on the corrected filling control parameters and the time compensation value. A preset excitation signal is applied to the filling equipment, and the input and output response data of the filling equipment to the excitation signal are collected.

[0029] The input and output response data are evaluated using a recursive least squares algorithm to estimate parameters and identify the flow coefficient, time constant, and valve response delay parameters of the filling equipment.

[0030] The first controller parameters are calculated based on the identified flow coefficient, time constant and valve response delay parameters. The first controller parameters are then subjected to rate of change limitation and absolute boundary constraint processing to generate the second controller parameters.

[0031] Using filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives, a multi-objective genetic algorithm is used to optimize the parameters of the second controller, and the optimal controller parameters are output.

[0032] Furthermore, the real-time monitoring and analysis of the operating status data of the filling equipment through multi-level fault detection to obtain fault diagnosis results, and the activation of corresponding fault-tolerant control strategies based on the fault diagnosis results, includes:

[0033] The operating status data of the filling equipment are sequentially subjected to signal-level detection, feature-level detection, and decision-level detection. Signal-level detection is used to perform threshold judgment and trend analysis on the raw sensor signals, feature-level detection is used to extract signal feature parameters and identify anomalies, and decision-level detection is used to fuse and judge multi-sensor data, and output the fault type and fault severity level.

[0034] The fault diagnosis result is determined by matching the fault type and fault severity level with a preset fault mode library.

[0035] Based on the fault diagnosis results, the fault-tolerant control strategy is generated by switching to a backup sensor or state estimation mode for sensor faults, activating redundant actuators for actuator faults, and switching to a degraded control mode for control algorithm faults.

[0036] Furthermore, the step of collaboratively optimizing the control parameters of the filling equipment and controlling the filling process using a hierarchical optimization strategy includes:

[0037] The current set of optimization variables and constraints are determined based on the fault-tolerant control strategy and the operating status data of the filling equipment.

[0038] The Bayesian optimization algorithm was used to perform offline global parameter optimization on historical filling data to determine the first optimization parameter;

[0039] Based on the first optimized parameters, the gradient descent method is used to fine-tune the parameters online in real time according to the preset filling cycle to generate the second optimized parameters;

[0040] Based on measurable disturbance data, the second optimization parameter is corrected through feedforward compensation calculation, a globally optimal filling control scheme is output, and filling control is executed according to the globally optimal filling control scheme.

[0041] Furthermore, the parameter estimation of the input / output response data using a recursive least squares algorithm, identifying the flow coefficient, time constant, and valve response delay parameters of the filling equipment, and the optimization of the second controller parameters using a multi-objective genetic algorithm with filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives, outputting the optimal controller parameters, including:

[0042] A Gaussian process regression model is constructed, with the control input and observable state variables in the input-output response data as the input of the Gaussian process regression model, and the filling volume measurement value as the output of the Gaussian process regression model. The posterior mean and posterior variance of the flow coefficient, time constant and valve response delay parameter are calculated by Bayesian inference.

[0043] The confidence index of each parameter is calculated based on the posterior variance. The spectrum of the excitation signal corresponding to the parameter whose confidence is lower than the preset threshold is adjusted to increase the signal energy of the sensitive frequency band of the parameter. The Gaussian process regression model is updated after collecting supplementary response data.

[0044] Monte Carlo sampling is performed on the posterior mean and the posterior variance to generate multiple sets of parameter samples. The stability margin and dynamic performance index of the closed-loop system are calculated for each set of parameter samples to determine the expected value and variance of each optimization objective function.

[0045] A weight matrix is ​​constructed based on the variance values ​​of each objective function. A first update step size is assigned to the objective function with a large variance value, and a second update step size is assigned to the objective function with a small variance value. The candidate solutions generated by the multi-objective genetic algorithm are locally refined using weighted Jacobi iteration. The second update step size is larger than the first update step size.

[0046] For the locally refined candidate solution set, calculate the performance confidence interval of each solution. Filter the candidate solutions by interval dominance relationship, retain robust solutions whose lower bound of the performance confidence interval dominates the upper bound of the confidence intervals of other solutions, and select the solution that meets the current production requirements from the robust solutions as the optimal controller parameters.

[0047] Furthermore, the step of collaboratively optimizing the control parameters of the filling equipment and controlling the filling process using a hierarchical optimization strategy includes:

[0048] The control parameters of the filling equipment are organized into a fourth-order control parameter tensor according to the dimensions of control type, parameter type, operating condition type, and optimization objective.

[0049] The fourth-order control parameter tensor is decomposed into a core tensor and four factor matrices. The core tensor has a lower dimension than the fourth-order control parameter tensor, and each factor matrix represents the main direction of change of the corresponding dimension.

[0050] The dimensions of the core tensor are determined based on the reconstruction error threshold and the cumulative energy ratio. The column vectors of each factor matrix are analyzed to extract the control mode features of each dimension.

[0051] In the low-dimensional space corresponding to the core tensor, the global control parameters are optimized offline using the Bayesian optimization algorithm to determine the low-dimensional optimal solution. The low-dimensional optimal solution is then reconstructed into the optimized parameters of the original parameter space through the factor matrix.

[0052] Based on the optimized parameters, the filling equipment is fine-tuned online and the feedforward compensation is corrected to output a globally optimal filling control scheme, and the filling control is executed according to the globally optimal filling control scheme.

[0053] Furthermore, after performing Tucker decomposition on the fourth-order control parameter tensor, the process further includes:

[0054] The control decision of the filling equipment is divided into four time scales: rapid control layer, single filling layer, batch adjustment layer and long-term optimization layer. Control decision tensors corresponding to each time scale are constructed respectively.

[0055] Define a scaling function between control decision tensors of adjacent time scales, and establish synchronization constraints between control decision tensors of each time scale through the scaling function.

[0056] The optimization problem with synchronization constraints is decomposed into independent subproblems at each time scale by using the alternating direction multiplier method. The controller at each time scale independently optimizes its own objective function, and the coordination between time scales is maintained through synchronization update steps and Lagrange multiplier update steps.

[0057] The first interaction mode and the second interaction mode are identified based on the numerical values ​​of each element in the core tensor. The synchronization constraints corresponding to the second interaction mode are relaxed, while the synchronization constraints corresponding to the first interaction mode are kept in strength.

[0058] For the missing elements corresponding to unobserved operating conditions in the fourth-order control parameter tensor, tensor completion calculation is performed based on the low-rank structure and physical feasibility constraints of the tensor to infer the control parameters under unobserved operating conditions. The completed control parameters are then used for the operating condition switching and filling control of the filling equipment.

[0059] The present invention also provides an adaptive control device for a filling equipment, comprising:

[0060] The temperature compensation module is used to collect real-time temperature data of the liquid temperature, gas temperature, sensor body temperature and ambient temperature of the filling equipment, calculate the influence coefficient of each temperature change on the filling volume through a hybrid compensation model that combines physical model and neural network, and output the corrected filling control parameters.

[0061] The dynamic time compensation module is used to perform filling and collect filling error data according to the corrected filling control parameters, perform empirical mode decomposition on the filling error data, separate the filling error into random error components, periodic error components and trend drift components, predict the error trend based on each error component and calculate the time compensation value.

[0062] The parameter self-tuning module is used to obtain the flow coefficient, time constant and valve response delay parameters of the filling equipment based on the corrected filling control parameters and the time compensation value, perform multi-objective optimization calculations on the flow coefficient, time constant and valve response delay parameters, and output the optimal controller parameters.

[0063] The fault diagnosis and fault tolerance module is used to control the operation of the filling equipment according to the optimal controller parameters, monitor and analyze the operating status data of the filling equipment in real time through multi-level fault detection, obtain fault diagnosis results, and start the corresponding fault tolerance control strategy based on the fault diagnosis results.

[0064] The collaborative optimization module is used to perform collaborative optimization and filling control on the control parameters of the filling equipment based on the fault-tolerant control strategy and the operating status data of the filling equipment through a hierarchical optimization strategy.

[0065] The beneficial effects of this invention are as follows:

[0066] This invention uses a hybrid compensation model that combines a physical model with a neural network to accurately calculate the influence coefficients of multiple temperature points, such as liquid temperature, gas temperature, sensor body temperature, and ambient temperature, on the filling volume. This solves the problem that the temperature compensation mechanism of traditional control systems is simple and cannot handle the complex coupling effects of multiple temperature points, enabling the system to maintain high-precision filling control within a wide temperature range of 10-40℃.

[0067] This invention uses empirical mode decomposition technology to separate filling errors into random error components, periodic error components, and trend drift components. Based on each error component, targeted error trend prediction and advance compensation are performed, which solves the problem that traditional control systems cannot cope with changes in the characteristics of different batches of products, and realizes fine adjustment within batches and rapid adaptation between batches.

[0068] This invention solves the problem of the lack of adaptive capability of control parameters in traditional control systems by automatically tuning controller parameters through online system parameter identification and multi-objective optimization algorithms, ensuring that the filling control system maintains optimal performance under various operating conditions. Compared with empirical tuning, it can reduce control error by 40%.

[0069] This invention solves the problem of limited fault diagnosis capability of traditional control systems by using multi-level fault detection (signal level, feature level, decision level) and adaptive fault-tolerant control strategy, so that the system can still maintain basic functions in the event of a fault, and reduce production losses caused by sudden failures by 90%.

[0070] This invention uses a hierarchical optimization strategy (offline optimization, online optimization, and feedforward compensation) to collaboratively optimize key system parameters, solving the problem of single optimization objectives in traditional control systems. It achieves the best balance between multiple objectives such as filling accuracy, response speed, control smoothness, and energy consumption, improving filling accuracy from ±2% to ±0.3%.

[0071] This invention provides uncertainty quantification for parameter estimation through Gaussian process regression and combines it with weighted Jacobian iteration for robust multi-objective optimization, enabling the system to maintain good performance under parameter uncertainty. Compared with traditional methods, it can reduce the risk of parameter adjustment by more than 95%.

[0072] This invention achieves dimensional reduction and pattern extraction of control parameters through Tucker decomposition, and combines four-focus tensor synchronization technology to achieve coordinated control across multiple time scales. This enables the system to systematically handle high-dimensional multivariable coupled problems, reducing computation time by 5-10 times while improving the quality of optimization results. Attached Figure Description

[0073] Figure 1This is a schematic diagram of the overall process of an adaptive control method and system for a filling equipment according to the present invention;

[0074] Figure 2 This is a schematic diagram of the adaptive control method and system for a filling equipment according to the present invention. Detailed Implementation

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

[0076] Example 1

[0077] like Figure 1 As shown, this embodiment provides an adaptive control method for a filling equipment, including the following steps:

[0078] Step S1: Collect real-time temperature data of liquid temperature, gas temperature, sensor body temperature and ambient temperature of the filling equipment, calculate the influence coefficient of each temperature change on the filling volume through a hybrid compensation model that combines physical model and neural network, and output the corrected filling control parameters.

[0079] In this embodiment, the temperature acquisition system deploys four high-precision temperature sensors to monitor the liquid temperature (using a Pt100 platinum resistance thermometer with an accuracy of ±0.1℃), the gas system temperature (using a K-type thermocouple), the sensor body temperature (integrated inside the pressure sensor), and the ambient temperature (using a DHT22 digital temperature and humidity sensor). Data acquisition from the four temperature sensors uses a unified time reference, with a sampling period set to 100ms to ensure the capture of dynamic temperature changes.

[0080] The hybrid compensation model consists of a physical model and a neural network. The physical model is based on fundamental principles of fluid mechanics. According to the Hagen-Poiseuille law, the relationship between the flow rate and pressure difference, pipe size, and liquid viscosity of the fluid in the middle layer of the pipe can be expressed as:

[0081]

[0082] in, Indicates flow rate. Indicates the pipe radius. Indicates pressure difference, Indicates the viscosity of the liquid. Indicates the length of the pipe.

[0083] The relationship between liquid viscosity and temperature is expressed using the modified Arrhenius equation:

[0084]

[0085] in, and For liquid-specific constants, This refers to absolute temperature.

[0086] The neural network part adopts a 3-layer backpropagation structure, which includes an input layer (4 nodes, corresponding to four temperature measurement points), a hidden layer (10 nodes, using the hyperbolic tangent activation function), and an output layer (1 node, representing the compensation coefficient).

[0087] The process of building and training a neural network is as follows:

[0088] (1) Data preprocessing: The four temperature inputs are normalized and mapped to the interval [-1, 1]. The normalization formula is as follows:

[0089]

[0090] in, This is the original temperature value. and These represent the minimum and maximum values ​​of the temperature measurement range, respectively.

[0091] Network initialization: The weights of the hidden and output layers are initialized using the Xavier method, with the initial weights following a uniform distribution.

[0092]

[0093] in Input the number of nodes. This is the number of output nodes. The bias term is initialized to 0.

[0094] (2) Forward Propagation: The output of the hidden layer is calculated as follows:

[0095]

[0096] in, This is the output of the j-th hidden node. The weights from the input layer to the hidden layer. For the i-th normalized temperature value, This is the bias for the hidden layer.

[0097] The output layer is calculated as follows:

[0098]

[0099] in, For compensation coefficient, The weights from the hidden layer to the output layer. This is the output layer bias.

[0100] (3) Loss function: The mean squared error loss function plus an L2 regularization term is used.

[0101]

[0102] in, The number of training samples, For the first Calculate the compensation coefficient for each sample. For the first The target compensation coefficient for each sample This is the regularization coefficient (value 0.001). The weights from the input layer to the hidden layer. These are the weights from the hidden layer to the output layer.

[0103] (4) Backpropagation and parameter update: The Adam optimizer is used, with the initial learning rate set to 0.01 and the momentum parameter... , The training process employs an early stopping method, using a validation set to monitor the training process. Training is stopped when the validation error fails to decrease for five consecutive epochs.

[0104] (5) Incremental learning: The initial weights of the model are preset based on the theoretical relationships of the physical model. For example, the initial weights corresponding to liquid temperature are set to a large positive value (reflecting the relationship that the viscosity decreases and the flow rate increases with the increase of temperature). Then, they are fine-tuned with actual data. The model is automatically updated every 8 hours, using the most recent 2000 sets of data. Each update only allows weight changes within ±10% to ensure smooth model evolution.

[0105] The output of the hybrid compensation model is a comprehensive compensation coefficient. Used to correct filling control time:

[0106]

[0107] in, This is the corrected time value. This is the original time value. This is the final compensation coefficient.

[0108] Step S2: Perform filling according to the corrected filling control parameters and collect filling error data. Perform empirical mode decomposition on the filling error data to separate the filling error into random error components, periodic error components and trend drift components. Based on each error component, predict the error trend and calculate the time compensation value.

[0109] In this embodiment, the filling volume is measured in real time using a high-precision electronic scale (resolution 0.01g, accuracy 0.02%) to measure the net weight of each filling container. The electronic scale communicates with the main control system via an RS-485 interface to transmit weight data in real time. Based on the liquid's density information, the weight data is converted into volume data, compared with the target filling volume, and the relative error is calculated.

[0110]

[0111] in, This represents the actual filling volume. For the target filling volume, This represents the relative error in filling volume.

[0112] The specific steps for processing error time series using Empirical Mode Decomposition (EMD) are as follows:

[0113] (1) Identify error sequences All local maxima and minima;

[0114] (2) Connect all maxima using cubic spline interpolation to form the upper envelope. Connect all the minimum values ​​to form the lower envelope. ;

[0115] (3) Calculate the average value of the upper and lower envelopes:

[0116]

[0117] (4) Extracting detail components:

[0118]

[0119] (5) Inspection Does it satisfy the Intrinsic Modulus Function (IMF) condition: the difference between the number of extreme points and the number of zero-crossing points does not exceed 1, and the mean of the upper and lower envelopes is zero at any time? If it satisfies this condition, then... For an IMF component If not satisfied, As a new Repeat steps (1)-(4);

[0120] (6) Subtract the obtained IMF from the original signal: ,Will As a new Repeat the above process until the residual component is reached. It is a monotonic function or its amplitude is less than a preset threshold.

[0121] Through EMD analysis, filling errors were separated into:

[0122]

[0123] Among them, the first 1-2 IMF components The errors are classified as high-frequency random errors, the middle IMF components are classified as periodic errors, and the last few IMFs and residual terms are classified as periodic errors. It is classified as a trend drift.

[0124] For the trend drift component, a weighted least squares method is used to fit a linear trend model:

[0125]

[0126] in, This is the filling sequence number. and These are the model parameters. To capture nonlinear variations, a curvature correction term is introduced:

[0127]

[0128] in, The curvature coefficient, Used as a reference point.

[0129] For periodic error components, the main frequency components are identified using Fast Fourier Transform (FFT). and its corresponding amplitude and phase Construct a periodic error prediction model:

[0130]

[0131] in, For the first Trend drift error during the second filling For the first Periodic errors in each filling process The number of main frequency components identified.

[0132] The overall error prediction value is:

[0133]

[0134] The time compensation value is calculated based on the prediction error. A hybrid feedback-feedforward control architecture is adopted, with the feedback component based on a PID control algorithm.

[0135]

[0136] in, Due to historical error, , , These are the proportional, integral, and derivative control parameters, respectively.

[0137] Feedforward components are based on prediction error:

[0138] in, This is the feedforward compensation coefficient (value range: 0.8-0.95).

[0139] The total compensation control amount is:

[0140]

[0141] Convert to time compensation value:

[0142]

[0143] in, The baseline filling time is used.

[0144] Step S3: Based on the corrected filling control parameters and the time compensation value, obtain the flow coefficient, time constant and valve response delay parameters of the filling equipment, perform multi-objective optimization calculation on the flow coefficient, time constant and valve response delay parameters, and output the optimal controller parameters.

[0145] In this embodiment, parameter self-tuning includes four steps: excitation signal design, system parameter identification, controller parameter calculation, and multi-objective optimization.

[0146] First, a pseudo-random binary sequence (PRBS) is designed as the excitation signal for system identification. PRBS generation is based on a shift register and XOR feedback:

[0147]

[0148] in, This represents the XOR operation. For register length, This is the tap position. The number of iterations in the sequence. For the first The PRBS sequence value of bits, For the first The PRBS sequence value of bits, For the first PRBS sequence value of bits; generation period is The sequence. The PRBS amplitude is controlled within ±5% of the rated pressure, and the clock cycle is set to 3-5 filling cycles.

[0149] The recursive least squares (RLS) algorithm is used to identify system parameters. The dynamic behavior of the filling system is described by difference equations:

[0150]

[0151] in, for The filling volume at any given time, for The filling volume at any given time, for The filling volume at any given time, for The amount of control at any given moment for The amount of control at any given moment for The amount of control at any given moment For discrete delay steps, These are the parameters to be identified.

[0152] Define parameter vector and regression vector The recursive formula for the RLS algorithm is:

[0153]

[0154]

[0155]

[0156] in, for The parameter vector estimate at time t. for The parameter vector estimate at time t. for Gain matrix at time step For regression vectors transpose, for The parameter covariance matrix at time t, for The parameter covariance matrix at time t, It is the identity matrix. The forgetting factor (values ​​0.95-0.99).

[0157] The identified difference equation parameters are converted into physical parameters:

[0158] Flow coefficient:

[0159] Time constant: ,in The system time constant, The sampling period is It is the natural logarithm function. It represents the principal pole of the characteristic polynomial of the difference equation.

[0160] Valve response delay: ,in This refers to the valve response delay time.

[0161] Based on the identified system model, the preliminary PID parameters are calculated using the improved Ziegler-Nichols method (CHR modified version):

[0162]

[0163]

[0164]

[0165] in, These are the initial values ​​for the proportional parameters of the PID controller. K'' is the system's master time constant, and K'' is the flow coefficient. This is the initial value of the integral time constant of the PID controller. This is the initial value of the derivative time constant of the PID controller.

[0166] Apply safety constraints to the calculated PID parameters:

[0167]

[0168]

[0169]

[0170] in, For the updated PID proportional parameters, The PID proportional parameters before the update. This is the updated PID integral time constant. The PID integral time constant before the update. The updated PID derivative time constant. The PID derivative time constant before the update. This is for absolute value operations.

[0171] A multi-objective optimization problem is constructed, with filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives:

[0172] Filling accuracy indicators:

[0173] Response speed metrics:

[0174] Control smoothness index:

[0175] Energy consumption indicators:

[0176] in, This refers to the filling accuracy index value. For the first Error value of each filling, This represents the sum of squares of all filling error values, where n'' is the total number of fillings. For response speed index value, For the first Response time for each filling cycle For the first Response time for each filling cycle It is the sum of the absolute values ​​of the differences between all two consecutive filling response times; To control the smoothness index value, For the first The change in the control quantity. It is the sum of the squares of all changes in control variables; Energy consumption index value, For the first Power value of the next filling process. For the first Duration of each filling cycle This represents the total energy consumption of all filling processes.

[0177] The NSGA-II multi-objective genetic algorithm is used to solve the problem.

[0178] minimize

[0179] Subject to: Parameter boundary constraints and stability constraints.

[0180] Among the decision variables Includes basic PID parameters and controller structure parameters. For multi-objective optimization, the objective function vector, Decision variables The corresponding filling accuracy index value, Decision variables The corresponding response speed index value, Decision variables The corresponding control smoothness index value, Decision variables The corresponding energy consumption index value; These are the controller structure parameters.

[0181] NSGA-II algorithm parameter settings: population size 100, number of iterations 100-200, crossover probability 0.9, mutation probability 0.1, using simulated binary crossover (SBX) and polynomial mutation.

[0182] After the algorithm generates the Pareto front, it selects the most suitable solution as the optimal controller parameters based on the current production requirements (such as high-precision mode, high-speed mode, and energy-saving mode).

[0183] Step S4: Control the operation of the filling equipment according to the optimal controller parameters, monitor and analyze the operating status data of the filling equipment in real time through multi-level fault detection, obtain the fault diagnosis results, and start the corresponding fault-tolerant control strategy based on the fault diagnosis results.

[0184] In this embodiment, multi-level fault detection includes three levels: signal-level detection, feature-level detection, and decision-level detection.

[0185] Signal-level detection layer: Signal-level detection performs threshold judgment and trend analysis on the raw sensor signal. This is applied to pressure sensor signals. Set upper and lower thresholds and ,when or A signal-level alarm is triggered. Simultaneously, the signal change rate is calculated.

[0186]

[0187] in, The rate of change of the pressure signal. for The pressure sensor signal value at any given time. for The pressure sensor signal value at any given time. The time step for calculating the rate of change; when the rate of change exceeds a preset threshold, it is determined to be an abnormal fluctuation.

[0188] Feature-level detection layer: Feature-level detection extracts signal feature parameters for anomaly identification. Statistical features are calculated for the pressure signal.

[0189] Mean:

[0190] Standard deviation:

[0191] Kuroshi:

[0192] in, The mean of the pressure signal samples. The standard deviation of the pressure signal samples. For pressure signal kurtosis, The total number of pressure signal samples. For the first One pressure signal sample value, It is the sum of all pressure signal sample values. This is the sum of squares of the deviations of all pressure signal samples from the mean. It is the fourth power sum of the deviations of all pressure signal samples from the mean; by comparing with normal benchmark values, abnormal characteristic parameters are identified.

[0193] Decision-level detection layer: Decision-level detection fuses and judges data from multiple sensors. It employs the DS evidence theory to fuse fault judgment results from different sensors. Let the confidence level of sensor 1 in judging fault A be... The confidence level of sensor 2 in determining fault A is 1. The confidence level after fusion is:

[0194]

[0195] in, The confidence level of fault A after fusion. For sensor 1, a pair of propositions The basic probability distribution, For sensor 2 pairs of propositions The basic probability distribution, For all intersections equal to The sum of the probability products of the propositional combinations. It is the sum of the probability products of all propositional combinations whose intersection is an empty set (i.e., conflict terms). It is an empty set; based on the fused confidence level and the preset threshold, the fault type (sensor fault, actuator fault, control algorithm fault, etc.) and fault severity level (prompt, attention, warning, alarm) are output.

[0196] Based on the fault diagnosis results, the corresponding fault-tolerant control strategy is activated:

[0197] (1) For sensor failure: If a pressure sensor failure is detected, automatically switch to the backup pressure sensor; if there is no backup sensor, activate the state estimation mode to estimate the pressure value based on data from other sensors (such as flow sensors and temperature sensors) and the system model.

[0198]

[0199] in, for Pressure estimate at time, Let be the system state estimation function. for The flow signal value at time. for Temperature signal value at time, This is the parameter vector of the state estimation model.

[0200] (2) For actuator failure: If the main filling valve is detected to be stuck, activate the backup filling valve; if there is no backup valve, adjust the control strategy to adapt to degraded operation, such as reducing the filling speed and increasing the filling time to compensate for insufficient valve response.

[0201] (3) In case of control algorithm failure: If the MPC solution times out or the parameter identification diverges, the system will automatically switch to the simpler and more reliable PID control mode to ensure the basic functions of the system.

[0202] Step S5: Based on the fault-tolerant control strategy and the operating status data of the filling equipment, the control parameters of the filling equipment are collaboratively optimized and the filling is controlled through a hierarchical optimization strategy.

[0203] In this embodiment, the hierarchical optimization strategy includes an offline optimization layer, an online optimization layer, and a feedforward compensation layer.

[0204] The offline optimization layer, based on historical big data, employs a Bayesian optimization algorithm to optimize global parameters. Bayesian optimization uses a Gaussian process (GP) to establish a surrogate model of the objective function.

[0205]

[0206] in, Let x'' be the objective function (such as comprehensive indicators like filling accuracy and energy consumption), and let x'' be the decision variable to be optimized (such as PID parameters, feedforward coefficients, etc.). Let be the mean function of the Gaussian process, describing the prior expectation of the objective function. The covariance function (kernel function) describes the correlation between the objective function values ​​at different sampling points. For any comparison sampling point, This indicates a Gaussian process distribution. The expected improved (EI) acquisition function is used to guide sampling.

[0207]

[0208] in, For the desired improvement value, For the expected operation, This represents the optimal objective function value among the currently evaluated samples. The operation involves taking the maximum value and zero; the point with the largest EI value is selected as the next evaluation point, and the optimization is iterated until convergence.

[0209] The online optimization layer, based on the offline optimization results, uses gradient descent to fine-tune the parameters in real time. It is updated every 10 filling cycles, and the gradient calculation uses a finite difference approximation.

[0210]

[0211] in, For the objective function For decision variables gradient, Indicates approximate equality. This is the step size for small perturbations (used for numerical approximation of gradients). The objective function value after adding perturbations to the decision variables. This represents the objective function value corresponding to the current decision variable. Parameter update formula:

[0212]

[0213] in, For the updated decision variables, For the decision variables before the update, The learning rate (ranging from 0.01 to 0.05) controls the step size for updating the parameters.

[0214] The feedforward compensation layer is designed to address measurable disturbances (such as temperature and liquid level changes) using a feedforward compensation mechanism. Control parameters are rapidly adjusted through lookup tables and interpolation calculations. A disturbance-compensation mapping table is established for temperature disturbances. The corresponding pressure compensation can be obtained by looking up the table. and time compensation :

[0215]

[0216]

[0217] in, To compensate for the final filling pressure, Based on the filling pressure, This is the pressure compensation amount; The final filling time after compensation. Based on the filling time, For time compensation amount, This is the difference between the current temperature and the reference temperature.

[0218] The three-layer optimization works in synergy: offline optimization provides the global optimal operating point, online optimization quickly responds to real-time changes, feedforward compensation actively suppresses known disturbances, and finally outputs the global optimal filling control scheme and executes the filling control.

[0219] Through the synergistic effect of the above steps S1-S5, this embodiment achieves a technological breakthrough by improving the filling accuracy from the traditional ±2% to ±0.3%, and maintains high-precision filling for a long time under various conditions, including a temperature range of 10-40℃, different batches of products, and various working conditions.

[0220] Example 2

[0221] This embodiment refines step S1 based on embodiment 1. The step of calculating the influence coefficient of temperature changes on filling volume using a hybrid compensation model combining a physical model and a neural network, and outputting the corrected filling control parameters, specifically includes the following sub-steps:

[0222] Step S1.1: Filter the real-time temperature data, input the filtered temperature data into the corresponding components of the physical model, and calculate the individual influence coefficient of each temperature point on the filling volume.

[0223] In this embodiment, a low-pass filter is used to filter the four temperature data streams, with a cutoff frequency set to 0.5Hz to eliminate high-frequency interference. The filter is a second-order Butterworth digital filter with the following transfer function:

[0224]

[0225] in, Let be the z-domain transfer function of the discrete-time system. For unit delay operators, These are the filter numerator coefficients. These are the denominator coefficients of the filter, calculated based on the cutoff frequency of 0.5Hz and the system sampling frequency.

[0226] The filtered temperature data are then input into the corresponding components of the physical model. For liquid temperature... Based on the viscosity-temperature relationship and the flow equation, the influence coefficient on the filling volume is calculated:

[0227]

[0228] in, This is the coefficient representing the influence of liquid temperature on filling flow rate. For traffic For liquid temperature The partial derivatives; Pi The radius of the filling pipeline, The pressure difference between the two ends of the pipeline, This refers to the length of the pipeline. For the dynamic viscosity of the liquid, This is the partial derivative of the reciprocal of viscosity with respect to the liquid temperature; , To correct the Arrhenius viscosity equation The empirical constant in This is the natural index term.

[0229] Similarly, calculate the gas temperature Sensor temperature and ambient temperature The influence coefficient.

[0230] Step S1.2: Input the real-time temperature data into the neural network for nonlinear compensation calculation. The neural network takes the liquid temperature, gas temperature, sensor body temperature and ambient temperature as inputs and outputs nonlinear coupling compensation coefficients.

[0231] In this embodiment, the normalized four-channel temperature data are:

[0232] Input the pre-trained neural network. The forward propagation process of the neural network is as described in Example 1, and the output is the nonlinear coupling compensation coefficient. .

[0233] This compensation coefficient captures the coupling effect between nonlinear effects and temperature that cannot be fully described by the physical model, for example:

[0234] Synergistic effect when liquid temperature and gas temperature rise simultaneously;

[0235] The interaction between sensor temperature drift and changes in ambient temperature;

[0236] Local flow field changes caused by temperature gradient.

[0237] Step S1.3: Based on the individual influence coefficients of each temperature point and the nonlinear coupling compensation coefficient, perform fusion calculation according to the adaptive weights of each temperature point to generate a comprehensive compensation coefficient, correct the filling control time based on the comprehensive compensation coefficient, and output the corrected filling control parameters.

[0238] In this embodiment, the adaptive weights are dynamically adjusted based on the temperature change rate and historical compensation effects:

[0239]

[0240] in, For the first Road temperature ( Corresponding liquid ,gas ,sensor ,environment Adaptive weights, For the first The basic weighting coefficient for road temperature (preset value, used to determine the initial weighting percentage of each temperature component). For the first Rate of change of road temperature For the first Road temperature at time intervals The change within, The time step for calculating the rate of temperature change is... For the first Historical compensation effectiveness evaluation factor for road temperature (adjusted based on feedback from past compensation accuracy; better compensation effect indicates better performance). The larger the value).

[0241] Weight normalization:

[0242]

[0243] in, For the first The normalized value of the road temperature adaptive weight (ensuring that the sum of the weights of all temperature components is 1). For the first The original adaptive weights for road temperature. For all road temperatures ( correspond , , , The sum of the original adaptive weights.

[0244] Calculation of comprehensive compensation coefficient:

[0245]

[0246] in, The overall compensation coefficient is the one without amplitude limiting processing. This is a weighted sum of the compensation coefficients corresponding to the four temperature streams (liquid, gas, sensor, and ambient). For the first Normalized weights for road temperature, For the first The compensation coefficient for the impact of road temperature. This is the weight for nonlinear compensation (usually a value of 0.2-0.3, used to balance the proportion of linear and nonlinear compensation). It is a nonlinear compensation coefficient (used to compensate for the nonlinear relationship between temperature and filling volume).

[0247] To ensure the stability of compensation, a limit is imposed on the comprehensive compensation coefficient:

[0248]

[0249] in, This is the final comprehensive compensation coefficient after amplitude limiting. This is the comprehensive compensation coefficient without amplitude limit. This is the comprehensive compensation coefficient used in the last application. To perform the maximum value operation, To minimize the value, this formula ensures that the adjustment range of the comprehensive compensation coefficient does not exceed ±5% of the previous adjustment, thus avoiding sudden changes in the compensation amount that could lead to system instability.

[0250] Revised filling control time:

[0251]

[0252] in, The revised final filling control time. The original filling control time before the correction. This is the comprehensive compensation coefficient after amplitude limiting.

[0253] Simultaneously output the corrected control parameters, such as the target pressure:

[0254]

[0255] in, The revised target filling pressure. The original filling target pressure before correction. This is a comprehensive compensation coefficient that has undergone amplitude limiting to ensure that the target pressure and filling time are corrected synchronously, thus maintaining stable filling accuracy.

[0256] These revised control parameters serve as inputs to step S2, used to perform filling and collect error data.

[0257] Example 3

[0258] This embodiment refines step S2 based on embodiment 1. The step of performing empirical mode decomposition on the filling error data, separating the filling error into random error components, periodic error components, and trend drift components, and then predicting the error trend and calculating the time compensation value based on each error component, specifically includes the following sub-steps:

[0259] Step S2.1: Perform iterative sieving processing on the filling error data. Extract the intrinsic modulus function by identifying extreme points and the mean of the envelope. Decompose the filling error data into multiple intrinsic modulus function components and residual components. According to the frequency characteristics of each component, classify the intrinsic modulus function components and residual components into the random error component, the periodic error component, and the trend drift component.

[0260] In this embodiment, the detailed iterative screening process of empirical mode decomposition is as follows:

[0261] (1) Initialization: Let the original error sequence be... Number of iterations ;

[0262] in, The initial error sequence for empirical mode decomposition. For time variables, This refers to the original error signal of the filling system (such as the relative error of the filling volume). This indicates that the iteration starts from step 0 (the initial state).

[0263] (2) Identification All local maxima and local minimum points ;

[0264] in, For the first Error sequence at the next iteration Let this be the set of times corresponding to all local maxima in the sequence. The index of the maximum point; The set of times corresponding to all local minima. This is the index of the minimum point.

[0265] (3) Construct the upper envelope using cubic spline interpolation. The lower envelope passes through all maxima. Pass through all local minimum points;

[0266] in, To pass through all maximum points The constructed upper envelope is used to wrap the upper edge fluctuations of the error sequence; To pass through all minimum points The constructed lower envelope is used to wrap the lower edge fluctuations of the error sequence; cubic spline interpolation is used to ensure that the envelope is smooth and fits the extreme value distribution, avoiding overshoot or distortion.

[0267] (4) Calculate the mean of the envelope:

[0268]

[0269] in, For the first The mean function of the upper and lower envelopes at the next iteration. The upper envelope, The lower envelope is used; this formula obtains the trend component of the error sequence by averaging, which prepares for separating local fluctuations.

[0270] (5) Extract detail components:

[0271]

[0272] in, For the first The latent intrinsic mode function (IMF) components extracted in the next iteration This is the error sequence for the current iteration. The mean of the envelope is used; this formula separates out the details of local high-frequency or periodic fluctuations by subtracting the trend component from the original sequence.

[0273] (6) Inspection Does it meet the IMF criteria?

[0274] Condition 1: The difference between the number of extreme points and the number of zero-crossing points does not exceed 1;

[0275] Condition 2: The absolute value of the mean of the upper and lower envelopes at any given time is less than the threshold. .

[0276] in, The components to be examined are detailed components; Condition 1 ensures the symmetry and uniform fluctuation of the components, which are fundamental characteristics of the IMF; Condition 2... The preset mean threshold is used to determine whether the component has been sufficiently stripped of its trend. If the absolute value of the mean is less than the threshold, it means that the component is symmetrical and tends to the zero mean.

[0277] (7) If the IMF conditions are met, then For the nth IMF component, calculate the residual: If not satisfied, assume , Return to step (2);

[0278] in, The first decomposition obtained One effective IMF component, Before removal The residual sequence after each IMF component This is the original error sequence. For the front The sum of all IMF components; if the condition is not met, then the current detailed component is... The new error sequence is used to continue iterative sieving, with the number of iterations... Add one.

[0279] (8) The remaining As a new Repeat steps (2)-(7) until the residual is a monotonic function or the amplitude is less than the preset threshold. .

[0280] in, The residual sequence after each decomposition serves as the original error sequence for the next iteration; The residual amplitude termination threshold is equal to 0.001 times the maximum absolute value of the original error signal; the decomposition process ends when the residual sequence is a monotonic function or its amplitude is less than this threshold.

[0281] After decomposition, we get:

[0282]

[0283] in, This is the original filling error signal. The total number of effective IMF components obtained from the decomposition. The summation of all IMF components. The final residual sequence (typically containing slow trends or system static biases) represents the fact that the original error can be completely decomposed into the sum of each IMF component and the final residual.

[0284] Classification based on the frequency characteristics of each IMF component:

[0285] High-frequency random error: the first 1-2 IMF components, average period Secondary filling;

[0286] Periodic error: the intermediate IMF component, average period Secondary filling;

[0287] Trend drift: the last few IMF components and residuals, Secondary filling or monotonous filling.

[0288] Among them, high-frequency random errors correspond to rapid fluctuations in the initial stage of decomposition, periodic errors correspond to regular fluctuations in the medium period, and trend drifts correspond to long-term slow changes or monotonic trends. The classification results are used for targeted compensation.

[0289] Frequency characteristics are determined by calculating the average period of each IMF:

[0290]

[0291] in, For the first The average period of each IMF component The total number of samples in the error sequence. For the first The number of zero-crossing points of each IMF component; this formula is based on the ratio of the number of zero-crossing points to the total number of samples to obtain the average period. The smaller the period, the higher the frequency, and the larger the period, the lower the frequency.

[0292] Step S2.2: Calculate the trend prediction value for the trend drift component using a prediction algorithm that combines linear extrapolation and curvature correction; perform frequency domain analysis on the periodic error component to extract the main frequency and amplitude parameters and calculate the period prediction value; integrate the trend prediction value and the period prediction value to determine the comprehensive error prediction value.

[0293] In this embodiment, the trend drift component Adopt the latest The filling data was used to fit a model using weighted least squares. The weights decay exponentially over time.

[0294]

[0295] in, For the first Weighting of sub-filling data It is a natural exponential function. The most recent fill count used for fitting. For the filling sequence number (the value range is...) to (corresponding to the last 30 fillings). This is the weight decay time constant, used to control the rate of weight decay. The smaller the value, the faster the weight decays, and the more emphasis is placed on recent filling data.

[0296] Fitting the objective function:

[0297]

[0298] in, Indicates the parameter , , Find the minimum value. , , The parameters of the trend model to be fitted ( The linear trend coefficient is... For constant terms, (nonlinear curvature coefficient). In order to address the recent Data from the first filling (from the first) Next to Summation (times), For the first Weighting of sub-filling data For the first Actual value of the trend drift component of the second filling. These are the predicted values ​​from the trend model. As a reference point, The objective function represents minimizing the sum of squared weighted errors, thus achieving the optimal fit, by representing the square of the prediction error.

[0299] Solving for the parameters , , Then, the trend forecast value:

[0300]

[0301] in, For the first Predicted values ​​of trend drift components for each filling cycle. To predict the step size (i.e., to predict the future step size) Trend error of the next filling This is the current filling sequence number. , , The model parameters obtained from the fitting are... As a reference point (usually taken) (that is, based on the current filling sequence number). The difference between the predicted index and the reference point is used to calculate the effect of the nonlinear curvature term.

[0302] For periodic error components Perform a Fast Fourier Transform:

[0303]

[0304] in, For periodic error components The frequency domain representation (Fourier transform result). This indicates the Fast Fourier Transform operation. For the first The actual value of the periodic error component of each filling. It is a frequency variable used to characterize the frequency characteristics of periodic errors.

[0305] Identify the largest amplitude in the spectrum Frequency components and its magnitude and phase :

[0306]

[0307]

[0308] in, For the first The amplitude of each major frequency component ( ), This is an absolute value operation used to extract the amplitude of a frequency domain signal; For the first The phase of the main frequency components, To perform phase operations, it is used to extract phase information from frequency domain signals; For the first The frequency values ​​of the three main frequency components are the three frequencies with the largest amplitude in the spectrum, representing the main fluctuation frequencies of the periodic error.

[0309] Cycle forecast value:

[0310]

[0311] in, For the first Predicted values ​​of periodic error components for each filling cycle. The number of main frequency components selected. To sum the contributions of the three main frequency components, For the first The amplitude of each frequency component, For the first The frequency of each frequency component For the first The phase of each frequency component, For the predicted filling sequence number, For the first The phase of each frequency component at the prediction index. It is a sine function used to reconstruct the time-domain signal of periodic errors.

[0312] Overall error prediction value:

[0313]

[0314] in, For the first The predicted value of the overall error of the first filling. For the first Predicted trend drift error value for each filling cycle. For the first The predicted periodic error values ​​for each filling cycle are added together to obtain the comprehensive error prediction result, which is used for subsequent control compensation.

[0315] For high-frequency random errors, due to their unpredictability, they are not included in the prediction model, but their statistical properties (standard deviation) are used to calculate the prediction confidence interval:

[0316]

[0317] in, This is the comprehensive error prediction value. To predict the upper and lower boundaries of the confidence interval (usually corresponding to a 95% confidence level). The standard deviation of the high-frequency random error component is used to characterize the fluctuation range of random errors and provide a basis for the fault-tolerant design of control strategies.

[0318] Step S2.3: Based on the comprehensive error prediction value, calculate the filling control adjustment amount by combining the feedback compensation component and the feedforward compensation component, and convert the filling control adjustment amount into the time compensation value.

[0319] In this embodiment, the feedback compensation component adopts an incremental PID algorithm:

[0320]

[0321] in, For the first The feedback compensation increment for each filling (i.e., the additional feedback control quantity that needs to be added this time). The proportional coefficient of the incremental PID controller. The integral coefficient is... These are the differential coefficients; For the first The actual error of a single filling (the difference between the actual filling volume and the target filling volume). For the first The actual error of each filling. For the first The actual error of each filling; the three terms in the formula correspond to the proportional increment term, integral increment term, and derivative increment term, respectively, which together constitute the incremental adjustment amount of this feedback control.

[0322] Cumulative control quantity:

[0323]

[0324] in, For the first Cumulative control amount of feedback compensation for each filling. For the first Cumulative control amount of feedback compensation for each filling. This is to compensate for the incremental change in the feedback; by accumulating incremental control quantities, continuous correction of errors is achieved, avoiding sudden changes in control quantities.

[0325] To prevent integral saturation, the cumulative control quantity is limited:

[0326]

[0327] in, To provide an upper limit threshold for the cumulative control quantity, This is the lower limit threshold for the cumulative control quantity. This indicates that the cumulative control amount will be limited to no more than the upper limit. This indicates that the control quantity is further limited to no less than the lower limit, ultimately ensuring that the cumulative control quantity is within a reasonable range and avoiding the control quantity exceeding the system's carrying capacity (i.e., integral saturation) due to the continuous accumulation of integral terms.

[0328] The feedforward compensation component is based on the prediction error and uses adaptive feedforward coefficients:

[0329]

[0330] in, For the first The adaptive feedforward coefficient for each filling (dynamically adjusted according to the prediction confidence level). The basic feedforward coefficients (preset initial feedforward weights). It is a natural exponential function. The adjustment factor (controlling the sensitivity of the feedforward coefficient as confidence level changes). For the first The prediction confidence index for the first filling (characterizing the reliability of the prediction error).

[0331] in, Based on the feedforward coefficients, As a regulating factor, To predict the confidence index:

[0332]

[0333] in, The confidence index is set to a value in the range of [0, 1], with the value being closer to 1 indicating a more reliable prediction. For the first The standard deviation of the prediction error for each filling (characterizing the range of fluctuation of the predicted value). The standard deviation of the total error (characterizing the overall fluctuation level of error throughout the filling process) is denoted by the formula. This formula quantifies the reliability of the prediction by the ratio of the standard deviation of the prediction error to the standard deviation of the total error. The smaller the ratio, the higher the confidence level.

[0334] Feedforward compensation amount:

[0335]

[0336] in, For the first The feedforward compensation amount for each filling cycle, with a negative sign indicating that the compensation direction is opposite to the prediction error direction (used to offset the prediction error). For adaptive feedforward coefficients, For the first The comprehensive error prediction value for each filling (trend error prediction value + periodic error prediction value).

[0337] Total compensation control amount:

[0338]

[0339] in, For the first Total compensation control amount for each filling. The cumulative control quantity is used for feedback compensation (to correct errors that have already occurred). The feedforward compensation amount (used to offset prediction errors in advance) and the two work together to achieve accurate error compensation.

[0340] Convert to time compensation value:

[0341]

[0342] in, For the first The time compensation value for each filling (i.e., the amount of filling time that needs to be adjusted). This is the total compensation control amount. This is the system's baseline filling time. The conversion factor from control quantity to time is calibrated based on the characteristics of the filling system, such as flow rate and pressure (usually with a value of 0.01-0.05), and is used to convert dimensionless control quantities into actual time adjustment quantities.

[0343] Final filling control time:

[0344]

[0345] in, For the first The final control time for each filling (the actual filling time). The nominal filling time (which has been initially corrected by temperature compensation in step S1 above). This is the time compensation value calculated in this experiment. By superimposing the compensation value, the prediction error and historical error can be accurately corrected to ensure filling accuracy.

[0346] Example 4

[0347] This embodiment refines step S3 based on embodiment 1. The step of obtaining the flow coefficient, time constant, and valve response delay parameters of the filling equipment, performing multi-objective optimization calculations on the flow coefficient, time constant, and valve response delay parameters, and outputting the optimal controller parameters specifically includes the following sub-steps:

[0348] Step S3.1: Determine the current working state of the filling system based on the corrected filling control parameters and the time compensation value, apply a preset excitation signal to the filling equipment, and collect the input and output response data of the filling equipment to the excitation signal.

[0349] In this embodiment, the current operating status includes current operating conditions such as temperature, pressure, and flow rate, as well as the current parameters of the controller. After confirming that the system is in a stable state (the standard deviation of the filling error is less than 0.5% for 10 consecutive filling cycles), the system identification process is initiated.

[0350] A PRBS (Pseudo-Random Binary Sequence) excitation signal is applied to the filling equipment. The PRBS sequence is generated using a 10-bit shift register, with a feedback polynomial of 1+x. 7 +x 10 A 1023-bit pseudo-random sequence is generated. The PRBS signal is superimposed on the current control quantity. Specifically, the nominal control quantity that is currently operating stably is used as the basis, and a PRBS signal with a fixed amplitude is superimposed to form the input control quantity used for system identification.

[0351] The nominal control quantity refers to the control quantity used when the system is currently in a stable operating state. The amplitude of the PRBS signal is 5% of the nominal control quantity. The value of each sampling point in the PRBS sequence is only -1 or +1. This binary pseudo-random signal can effectively stimulate the dynamic response of the system, making it easier to identify the system characteristics later.

[0352] The PRBS clock cycle is set to 4 fill cycles, meaning the PRBS signal value switches once every 4 fill cycles. Theoretically, a complete PRBS sequence requires 1023 × 4 = 4092 fill cycles to complete. However, in practical engineering applications, only 200-300 fill cycles are needed to obtain sufficient system recognition accuracy. Therefore, the complete sequence is not required; only the first 75 bits of the PRBS sequence are used for excitation.

[0353] Collect input and output response data. The specific data to be collected is as follows:

[0354] Input data: This refers to the control quantity after applying PRBS excitation. The pressure control quantity or time control quantity can be selected according to the characteristics of the filling system.

[0355] Output data: The actual filling volume output by the filling equipment under the corresponding input control quantity;

[0356] Auxiliary data includes directly observable system state variables such as current temperature and liquid level, which are used to correct the identification results and improve the identification accuracy.

[0357] All collected data are timestamped, and the sampling period is consistent with the filling period. Typically, a filling is completed every 2-5 seconds, meaning data is collected every 2-5 seconds.

[0358] Step S3.2: Use the recursive least squares algorithm to estimate the parameters of the input and output response data, and identify the flow coefficient, time constant and valve response delay parameters of the filling equipment.

[0359] This embodiment employs an autoregressive exogenous input model to describe the dynamic characteristics of the filling system. This model represents the current filling volume as a linear combination of past filling volumes, past control inputs, and the current modeling error. The model includes four key parameters that reflect the system's dynamic response characteristics.

[0360] Recursive least squares is an online parameter estimation method that continuously updates parameter estimates as new data arrives. Starting with initial guesses, the algorithm gradually adjusts the parameters to minimize the prediction error by comparing the model's predictions with the actual measurements. A covariance matrix is ​​introduced to track the uncertainty of the parameter estimates; a smaller covariance matrix indicates a more reliable parameter estimate.

[0361] To enable the algorithm to adapt to slow changes in system characteristics, this embodiment employs a variable forgetting factor strategy. The forgetting factor uses a variable forgetting factor strategy:

[0362]

[0363] in, , , , This represents the prediction error.

[0364] The forgetting factor determines how much importance the algorithm places on historical and new data. When the system is stable, The parameter remains stable when it is close to 0.99; however, when environmental changes cause the prediction error to increase, Lowering it to 0.95 speeds up parameter tracking.

[0365] After convergence, the parameters in the difference equation form need to be converted into system parameters with clear physical meaning. The flow coefficient reflects the steady-state effect of the control input on the filling volume; the time constant reflects the speed of the system response and is obtained by solving the characteristic equation of the model; the valve response delay is determined by analyzing the cross-correlation function of the input and output signals, reflecting the time lag between the control action and the actual effect.

[0366] Step S3.3: Calculate the first controller parameters based on the identified flow coefficient, time constant and valve response delay parameters, and perform change rate limiting and absolute boundary constraint processing on the first controller parameters to generate the second controller parameters.

[0367] Based on the identified system parameters, this embodiment employs an improved Ziegler-Nichols method to calculate preliminary controller parameters. This method first converts the identified discrete model into a continuous transfer function form, and then calculates the proportional gain, integral time, and derivative time according to specific empirical formulas based on the system's flow coefficient, time constant, and delay time.

[0368] To ensure the smoothness and safety of parameter adjustment, two layers of constraints are applied to the parameters obtained from the initial calculations:

[0369] The first layer: rate of change limitation, preventing parameters from changing too rapidly relative to their current values. If a newly calculated parameter increases by more than 20% or decreases by more than 20% relative to the current parameter, the change is limited to this range. This avoids system oscillations or instability caused by sudden parameter changes.

[0370] The second layer: absolute boundary constraints, used to ensure that controller parameters are always within physically reasonable and safe operating ranges. For each controller parameter, based on the system's own physical characteristics and safe operating requirements, corresponding upper and lower thresholds are preset to strictly limit the parameter values.

[0371] Specifically, the preset boundary ranges for each key controller parameter are as follows:

[0372] The proportional gain is set to a range of 0.1 to 10;

[0373] The integration time is set to range from 0.5s to 50s;

[0374] The range of the differential time is set from 0.01s to 5s.

[0375] After the above two layers of constraint processing, the generated second controller parameters will serve as the initial solution for the subsequent multi-objective optimization process, providing a foundation for the efficient implementation of the optimization process.

[0376] Step S3.4: Using filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives, the second controller parameters are optimized using a multi-objective genetic algorithm, and the optimal controller parameters are output.

[0377] This embodiment defines four mutually constraining optimization objectives:

[0378] (1) Filling accuracy target: The deviation between the actual filling volume and the target value should be as small as possible, which is measured by calculating the root mean square value of the relative error during the simulation process.

[0379] (2) Response speed target: The system is required to reach a steady state quickly, which is measured by measuring the settling time of the step response.

[0380] (3) Control smoothness target: It is required that the control action changes smoothly and avoids drastic fluctuations. It is measured by calculating the root mean square value of the control increment.

[0381] (4) Energy consumption target: The energy consumption of the control process should be as low as possible, which is measured by calculating the average value of the absolute value of the control quantity.

[0382] This embodiment employs the Non-Dominated Sorting Genetic Algorithm (NSGA-II) for multi-objective optimization. This algorithm maintains a population containing 100 candidate solutions and searches for the optimal solution set by simulating a natural evolutionary process. NSGA-II Algorithm Flow:

[0383] (1) Initialize the population: Randomly generate an initial population to ensure that the parameters of the second controller are included in the initial population, providing a good starting point for optimization.

[0384] (2) Evaluation of fitness: For each candidate solution, the algorithm uses the identified system model to perform closed-loop simulation and calculates the four objective function values.

[0385] (3) Non-dominated sorting: The population is sorted non-dominated according to the Pareto dominance relationship, and the candidate solutions are divided into different levels.

[0386] (4) Crowding distance calculation: Within the same level, calculate the crowding distance of each solution to maintain the diversity of solutions. Solutions with large crowding distances are more isolated in the target space. Retaining these solutions helps to maintain the wide distribution of the solution set.

[0387] (5) Selection: The algorithm selects parent individuals from the current population through a tournament selection mechanism, giving priority to individuals with lower ranks, and selecting individuals with larger crowding distances when ranks are the same.

[0388] (6) Parent generation: The selected parent generation generates offspring through a simulated binary crossover operator with a crossover probability of 0.9. The crossover process simulates the process of biological chromosome exchange.

[0389] (7) Offspring: The offspring also need to be processed by the polynomial mutation operator. The mutation probability is set to 0.1. The mutation introduces random perturbation to explore new solution space regions.

[0390] (8) Merge the parent and offspring generations, recalculate the non-dominated sorting and crowding distance, and select the top 100 individuals as the new generation population;

[0391] (9) Repeat steps (2)-(8) for 150 iterations or until the Pareto front converges (the front change is less than 1% for 20 consecutive iterations).

[0392] After optimization, a Pareto optimal solution set is obtained, where each solution achieves a different balance among the four objectives. The final solution is selected based on current production requirements.

[0393] High-precision mode: Selects the solution with the lowest precision target.

[0394] High-speed mode: The solution that minimizes the speed target.

[0395] Balanced mode: The TOPSIS method is used for comprehensive scoring, and the solution closest to the ideal point is selected.

[0396] This method defines an ideal point and an anti-ideal point, calculates the distance of each solution to these two points, and selects the solution that is closer to the ideal point and farther from the anti-ideal point.

[0397] Example 5

[0398] This embodiment refines step S4 based on embodiment 1. The step of real-time monitoring and analysis of the operating status data of the filling equipment through multi-level fault detection to obtain fault diagnosis results, and then activating the corresponding fault-tolerant control strategy based on the fault diagnosis results, specifically includes the following sub-steps:

[0399] Step S4.1: The operating status data of the filling equipment is sequentially subjected to signal-level detection, feature-level detection, and decision-level detection. Signal-level detection is used to perform threshold judgment and trend analysis on the original sensor signals, feature-level detection is used to extract signal feature parameters and identify anomalies, and decision-level detection is used to perform fusion judgment on multi-sensor data and output the fault type and fault severity level.

[0400] This embodiment establishes a three-tiered fault detection system at the signal level, feature level, and decision level to achieve comprehensive monitoring of the operating status of filling equipment.

[0401] Signal-level detection: The most basic detection layer, directly monitoring the raw sensor signals in real time. For critical sensor signals such as pressure, flow, and temperature, the system sets upper and lower thresholds; an alarm is immediately triggered if the signal exceeds the normal range. In addition to absolute threshold judgment, signal-level detection also performs trend analysis, using a sliding window linear regression method to calculate the rate of change of the signal. When the rate of change exceeds the maximum allowable value, it indicates that the system state is changing rapidly and an anomaly may exist, triggering a trend alarm. This method can provide early warning before the fault fully manifests.

[0402] Feature-level detection involves a deeper analysis of the signal, extracting statistical and frequency domain features that reflect its intrinsic characteristics. Time-domain statistical features include mean, standard deviation, kurtosis, and skewness. The mean reflects the average level of the signal, the standard deviation reflects the degree of fluctuation, kurtosis reflects the sharpness of the signal distribution, and skewness reflects the symmetry of the signal distribution. Frequency domain features are extracted through Fast Fourier Transform, including the dominant frequency and the energy ratio of specific frequency bands. These features can reveal hidden abnormal patterns in the signal.

[0403] Feature-level detection compares extracted features with normal benchmark values ​​and uses statistical tests to identify anomalies. For example, if the current mean deviates from the benchmark mean by more than three times the benchmark standard deviation, or the current standard deviation exceeds twice the benchmark standard deviation, or the kurtosis increases abnormally, these are all considered feature anomalies. This detection method based on statistical distribution has strong robustness and can adapt to reasonable fluctuations under normal operating conditions.

[0404] Decision-level detection: This approach integrates the judgment results from multiple sensors and employs evidence theory for comprehensive decision-making. First, a fault hypothesis space is defined, encompassing various possible specific fault types, such as pressure sensor failure, valve leakage, and pipeline blockage, as well as an uncertainty category. Each sensor, based on its signal and feature analysis results, provides a Basic Probability Assignment (BPA), representing the confidence level for various fault hypotheses.

[0405] The evidence theory's combination rule fuses evidence from different sensors. The fusion process considers the consistency and conflict between pieces of evidence; consistent evidence reinforces each other, while conflicting evidence weakens each other. By progressively fusing evidence from pressure, flow, and temperature sensors, a comprehensive fault confidence distribution is obtained. The fault type with the highest confidence is identified as the most probable fault.

[0406] Based on the highest confidence level, the system classifies the severity of the fault into four levels.

[0407] A confidence level greater than 0.5 indicates a possible anomaly, but the cause is uncertain; attention is advised.

[0408] Confidence level between 0.5 and 0.7: Level of concern, indicating a high probability of a fault and requiring close monitoring;

[0409] Confidence level between 0.7 and 0.9: Warning level, indicating a high probability of a fault, and countermeasures should be prepared;

[0410] Confidence level greater than or equal to 0.9: Alarm level, indicating that a fault is almost certain and immediate action must be taken.

[0411] Step S4.2: Match the fault type and fault severity level to a preset fault mode library to determine the fault diagnosis result.

[0412] In this embodiment, the fault mode library contains approximately 100 typical fault modes, and each fault mode record includes:

[0413] Fault codes: unique identifiers (e.g., F001, F002, ...);

[0414] Fault Description: Describe the fault phenomenon and cause in words;

[0415] Fault characteristics: typical sensor signal characteristics and statistical parameters;

[0416] Diagnostic methods: Recommended further diagnostic steps;

[0417] Fault tolerance measures: corresponding fault tolerance control strategies;

[0418] Maintenance recommendations: Recommended maintenance procedures and spare parts.

[0419] Fault codes are unique identifiers that facilitate quick retrieval and referencing. Fault descriptions use text to describe the phenomena, possible causes, and impacts of the fault. Typical characteristics record the characteristic manifestations of the fault on various sensor signals and statistical parameters, providing a basis for matching.

[0420] Fault mode matching employs a feature similarity calculation method. For a detected fault, the system extracts its feature vector, including the signal features and statistical parameters of each sensor. Then, it calculates the similarity between this feature vector and the feature vector of each fault mode in the fault mode library. The similarity calculation uses a weighted combination approach, assigning different weights to different features according to their importance, and employing a Gaussian similarity function to measure the similarity of individual features.

[0421] The system selects the fault pattern with the highest similarity as the matching result. If the highest similarity is lower than the threshold (e.g., 0.6), it indicates that the detected fault is not very similar to the known fault patterns and may be an unknown fault type. The system will record the characteristics of the fault and prompt manual diagnosis, while adding it to the database as a potential new fault pattern.

[0422] After determining the fault diagnosis result, the system outputs the fault type, confidence level, severity level, and recommended actions. The confidence level combines the fusion confidence of decision-level detection and the similarity of pattern matching, reflecting the reliability of the diagnosis result. Recommended actions are extracted from the fault mode library, including fault-tolerant control strategies and maintenance suggestions, providing guidance for subsequent processing.

[0423] Step S4.3: Based on the fault diagnosis results, switch to backup sensor or state estimation mode for sensor faults, activate redundant actuators for actuator faults, switch to degraded control mode for control algorithm faults, and generate the fault-tolerant control strategy.

[0424] This embodiment implements differentiated fault-tolerance strategies for different types of faults to ensure that the system can still operate safely and reliably in the event of a fault:

[0425] (1) Sensor fault tolerance:

[0426] If a pressure sensor failure is detected (e.g., constant output, excessive noise, or significant inconsistency with other sensors), the preferred solution is to switch to the backup pressure sensor. The filling system is equipped with dual redundant pressure sensors, automatically switching to the backup sensor to continue operation when the primary sensor fails. If there is no backup sensor or the backup sensor also fails, a state estimation mode is activated. This mode estimates the measured value of the failed sensor based on data from other available sensors and a system model. For example, pressure values ​​can be estimated using data from flow and temperature sensors, combined with regression models trained using fluid dynamics equations and historical data.

[0427] In state estimation mode, due to the uncertainty of the estimated value, the system adds estimation uncertainty compensation, automatically adjusting the controller parameters to more conservative settings (such as increasing the integral time and decreasing the proportional gain), and reducing the filling speed to improve reliability. Although these measures sacrifice some performance, they ensure that the system can still operate stably in the event of sensor failure.

[0428] (2) Actuator fault tolerance:

[0429] If the main filling valve is detected to be stuck or responding abnormally, the preferred solution is to activate the backup filling valve. The system is equipped with a dual-valve redundancy design, allowing switching to the backup valve to continue filling in the event of a main valve failure. If there is no backup valve or the backup valve also fails, the system will adjust its control strategy to adapt to degraded operation. Specific measures include reducing the filling speed to 70% of the normal speed, increasing the filling time compensation to 1.3 times the normal time, and adjusting the controller parameters to conservative settings. These measures reduce the performance requirements of the actuators, enabling the system to continue operating even if some actuators fail.

[0430] If unstable pressure is detected in the pneumatic system, the system will activate the pressure stabilization compensation mode, increase the integral gain of the pressure feedback control to speed up the pressure regulation response, and adjust the filling strategy to adopt a segmented filling method, which is fast at first and then slow, reducing the requirements for air pressure stability.

[0431] (3) Fault tolerance of control algorithm:

[0432] If an MPC solution timeout, parameter identification divergence, or abnormal controller output is detected, the system will automatically switch to a degraded control mode (from MPC to classic PID control). Specifically, the PID controller uses the most recently verified parameter configuration, while the system performance requirements are reduced, allowing a larger filling error range (from ±0.3% to ±0.8%) to ensure stable system operation.

[0433] If an abnormal parameter identification result is detected (such as the identified parameter exceeding the physical reasonable range), the system will refuse to use the new identification result, keep using the parameters that passed the previous verification, trigger the re-identification process, adjust the excitation signal parameters and increase the identification data volume, record the abnormal event and prompt manual inspection.

[0434] The generated fault-tolerant control strategies include:

[0435] Control Mode: The currently active control mode (normal mode, sensor fault-tolerant mode, actuator fault-tolerant mode, degraded control mode);

[0436] Parameter configuration: Currently used controller parameters and system parameters;

[0437] Performance limitations: Current performance requirements (accuracy range, speed limits, etc.);

[0438] Recovery conditions: The conditions for recovering from fault-tolerant mode to normal mode (such as fault elimination, continuous stable operation time, etc.);

[0439] Alarm information: Alarm and suggestion information sent to operators and maintenance personnel;

[0440] The fault-tolerant control strategy is executed automatically, and all fault-tolerant actions and system state changes are recorded for subsequent analysis and system optimization.

[0441] Example 6

[0442] This embodiment refines step S5 based on embodiment 1. The step of collaboratively optimizing the control parameters of the filling equipment and controlling the filling process through a hierarchical optimization strategy specifically includes the following sub-steps:

[0443] Step S5.1: Determine the current set of optimization variables and constraints based on the fault-tolerant control strategy and the operating status data of the filling equipment.

[0444] This embodiment dynamically determines the set of optimized variables and constraints based on the current fault-tolerant control strategy and system operating status. In normal mode, the set of optimized variables includes all adjustable parameters: target pressure, target filling time, controller proportional gain integral time and derivative time, temperature compensation coefficient, feedforward compensation gain, batch adaptation parameters, etc. These variables collectively determine the control performance of the filling system.

[0445] In fault-tolerant mode, some variables are fixed or subject to stricter constraints. In sensor fault-tolerant mode, due to the use of state estimation, the parameters of the estimation model are fixed, and only the controller parameters are optimized to accommodate the uncertainty of the estimated values. In actuator fault-tolerant mode, considering the limited performance of the actuator, the pressure and speed ranges are limited to more conservative ranges, and the compensation parameters are mainly optimized to compensate for the degradation of actuator performance. In degraded control mode, the controller parameters are fixed to a conservative configuration, and only the target setpoint is optimized.

[0446] Constraints encompass multiple levels. Physical constraints ensure parameters remain within the equipment's physical capabilities; for example, pressure cannot exceed the capacity of pipes and valves, and time cannot be less than the minimum time for mechanical action. Stability constraints, through analysis of the closed-loop system's characteristic equations, ensure that the real parts of all eigenvalues ​​are less than zero, guaranteeing system stability. Performance constraints limit the standard deviation of errors and settling time to acceptable ranges. Safety constraints limit the rate of pressure change and the magnitude of control action changes, preventing excessively rapid changes from causing equipment damage or safety accidents.

[0447] Depending on the current fault-tolerance strategy, constraints may need to be tightened. In sensor fault-tolerance mode, performance constraints are relaxed to avoid system instability caused by excessive pursuit of accuracy due to increased measurement uncertainty. In actuator fault-tolerance mode, physical constraints are tightened to protect damaged actuators. In degraded control mode, safety constraints are tightened to ensure system safety under simplified control.

[0448] Step S5.2: Use the Bayesian optimization algorithm to perform offline global parameter optimization on the historical filling data to determine the first optimization parameter.

[0449] In this embodiment, offline optimization is based on historical filling data from the past 7 days (approximately 100,000 filling records). The goal of offline optimization is to find the parameter configuration that performs optimally under historical operating conditions.

[0450] The optimization objective function is a weighted combination of four performance metrics: filling accuracy, response speed, control smoothness, and energy consumption. Each performance metric is calculated based on historical data, and its weight is set according to current production needs. For example, when pursuing high accuracy, the accuracy metric is given a larger weight; when pursuing high efficiency, the speed and energy consumption metrics are given larger weights.

[0451] Bayesian optimization is an efficient global optimization method, particularly suitable for situations where objective function evaluation is costly. The core idea is to construct a probabilistic surrogate model of the objective function, using this surrogate model to guide sampling, thus achieving a balance between exploring unknown regions and utilizing known optimal regions.

[0452] The Bayesian optimization algorithm first employs Latin Hypercube Sampling (LHS) to generate 10 uniformly distributed initial sampling points in the parameter space. For each sampling point, the objective function value is evaluated through simulation using historical data. Then, a Gaussian process surrogate model is constructed based on these initial data. The Gaussian process assumes that the objective function follows a Gaussian process distribution, completely determined by the mean function and covariance function (kernel function). This embodiment uses a constant mean function and a Matérn 5 / 2 kernel function, which can flexibly model functions with different smoothness levels.

[0453] The objective function is defined as:

[0454]

[0455] in, Configure parameters The corresponding offline comprehensive performance objective function value (scalar) is the core indicator to be optimized; for Parameter configuration vector in dimensional parameter space ( (for parameter dimensions) , , , These are four types of performance indicators (scalars) obtained from historical data statistical calculations: production efficiency, stability, maintenance cost, and energy consumption. , , , The weighting coefficients (non-negative scalars) for the four types of performance indicators satisfy the following conditions: It is dynamically set according to current production needs.

[0456] The mean function uses a constant mean form:

[0457]

[0458] in, The global constant mean (scalar) represents the average level of the objective function in the parameter space.

[0459] The covariance function uses the Matérn 5 / 2 kernel function, and its expression is:

[0460]

[0461] Among them, the distance term ; The signal variance (positive scalar) controls the overall amplitude of the kernel function, reflecting the global fluctuation of the objective function; The normalized Euclidean distance (scalar) represents the parameter configuration. and Similarity between them; , They are parameter vectors , The Dimensional components (scalars); For the first The length scale parameter (positive scalar) of the dimension parameter controls the range of influence of this dimension parameter on the relevance of the objective function.

[0462] The Gaussian process model contains several hyperparameters, such as the length scale of the kernel function and the signal variance. These hyperparameters are optimized by maximizing the marginal likelihood to best fit the observed data. Hyperparameter optimization is performed using the L-BFGS-B (Limited-memory Broyden–Fletcher–Goldfarb–Shanno-Bound) algorithm with 20 iterations.

[0463] With the Gaussian process surrogate model, the algorithm uses the expected improvement sampling function to determine the next evaluation point. The expected improvement quantifies the anticipated improvement in the objective function that sampling at a given point might bring, comprehensively considering both the predicted mean (utilization) and the predicted uncertainty (exploration) at that point. The sampling function takes on larger values ​​in regions with low predicted mean and high uncertainty, guiding the algorithm to sample in these promising regions.

[0464] The Differential Evolution (DE) algorithm is used to maximize the expected improvement of the acquisition function and find the next evaluation point:

[0465]

[0466] in, This represents the parameter configuration vector that maximizes the EI acquisition function. .

[0467] At this point, evaluate the objective function, add new data to the training set, re-optimize the Gaussian process hyperparameters, and update the model. Repeat this process, iterating 50 times or until the expected improvement is less than a threshold, indicating that the optimal solution has been found.

[0468] After convergence, the point with the minimum objective function value is selected as the first optimization parameter. Compared to traditional grid search or random search, Bayesian optimization can find a better solution with fewer evaluations, making it particularly suitable for situations with high parameter dimensionality.

[0469] Step S5.3: Based on the first optimized parameters, the parameters are fine-tuned online in real time according to the preset filling cycle using the gradient descent method to generate the second optimized parameters.

[0470] Based on the initial optimized parameters, this embodiment employs the gradient descent method for online real-time fine-tuning, updating the parameters every 10 filling cycles. The purpose of online optimization is to adapt the parameters to the current actual operating conditions and compensate for potential deviations in offline optimization based on historical data.

[0471] The gradient calculation employs a finite difference approximation method. For each parameter, a small perturbation is applied near its current value, and the change in the objective function is observed. The gradient in the direction of that parameter is approximated by dividing the change by the perturbation. To reduce computational cost, a stochastic coordinate descent strategy is used, calculating the gradient in only one or two coordinate directions for each update, iterating through all coordinates. This significantly reduces computational cost while ensuring convergence.

[0472] Parameter updates employ a momentum-based gradient descent method. The momentum term accumulates historical gradient information, giving parameter updates inertia, which accelerates convergence and reduces oscillations. The momentum coefficient is set to 0.9, and the learning rate is set to 0.02. The learning rate determines the step size for each update; a proper setting can achieve a balance between convergence speed and stability.

[0473] To prevent online optimization from deviating too far from the offline optimal solution, a regularization term is introduced. The regularization term penalizes the distance between the current parameter and the first optimized parameter, ensuring that online optimization adapts to the current operating conditions while maintaining a connection to the global optimal solution. The regularization coefficient is set to 0.1, balancing adaptability and stability.

[0474] Online optimization also implements adaptive learning rate adjustment, gradually decreasing the learning rate as the number of updates increases. Initially, the learning rate is larger, resulting in faster parameter adjustments; later, the learning rate is smaller, allowing for more refined parameter adjustments, which is beneficial for convergence to a local optimum. After online optimization, a second set of optimized parameters is generated. This parameter inherits the global optimum of offline optimization while also adapting to the specific characteristics of the current operating conditions.

[0475] Step S5.4: Based on the measurable disturbance data, the second optimization parameter is corrected through feedforward compensation calculation, a globally optimal filling control scheme is output, and filling control is executed according to the globally optimal filling control scheme.

[0476] This embodiment uses feedforward compensation based on measurable disturbance data to further refine the second optimization parameters. Measurable disturbances include changes in ambient temperature, liquid level in the storage tank, and gas supply pressure. These disturbances can affect the filling process, but since they can be measured in real time, their impact can be actively offset through feedforward compensation.

[0477] The system establishes a mapping relationship between disturbances and compensation amounts, obtained through offline experiments or simulations. For temperature disturbances, compensation includes pressure compensation and time compensation, with the compensation amount exhibiting a quadratic function relationship with the temperature deviation. For liquid level disturbances, pressure is primarily affected, and the compensation amount shows a linear relationship with the liquid level deviation. For gas supply pressure disturbances, response speed is primarily affected, requiring time compensation, with the compensation amount showing a linear relationship with the pressure deviation.

[0478] The coefficients of the mapping relationship are obtained by fitting historical data using the least squares method. In actual operation, the system measures each disturbance in real time, calculates the deviation relative to the reference value, calculates each compensation amount according to the mapping relationship, and then superimposes all compensation amounts to obtain the total feedforward compensation.

[0479] The corrected final control parameters equal the second optimized parameters plus feedforward compensation. The system outputs the globally optimal filling control scheme, including complete information such as target pressure, target time, controller parameters, feedforward compensation, and temperature compensation coefficients.

[0480] According to the scheme, the specific steps of the control system to perform filling are as follows: First, adjust the reference filling time according to the temperature compensation coefficient, then apply feedforward compensation to obtain the set value, the controller calculates the control quantity according to the set value and actual feedback, the control quantity is converted into valve opening or pressure command to drive the filling valve to perform filling, the filling volume is monitored in real time and the error is calculated, and the feedback is given to the control system for optimization of the next filling.

[0481] Through a three-tiered approach of offline optimization, online fine-tuning, and feedforward compensation, the system achieves global optimization of filling control parameters, ensuring high-precision filling under various operating conditions.

[0482] Example 7

[0483] This embodiment enhances steps S3.2 and S3.4 based on embodiment 4. Specifically, it uses a recursive least squares algorithm to estimate parameters from the input / output response data, identifying the flow coefficient, time constant, and valve response delay parameters of the filling equipment. Furthermore, it optimizes the second controller parameters using a multi-objective genetic algorithm, with filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives, and outputs the optimal controller parameters. The specific sub-steps include:

[0484] Step S3.2B: Construct a Gaussian process regression model, using the control input and observable state variables in the input-output response data as the input of the Gaussian process regression model, and the filling volume measurement value as the output of the Gaussian process regression model. Calculate the posterior mean and posterior variance of the flow coefficient, time constant, and valve response delay parameter through Bayesian inference.

[0485] This embodiment constructs a Gaussian process regression model to identify system parameters. Compared with the traditional recursive least squares method, Gaussian process regression can not only provide point estimates of parameters, but also quantify the uncertainty of the estimates and provide the probability distribution of the parameters.

[0486] The Gaussian process regression model models the relationship between filling volume and control inputs and observable state variables as a Gaussian process. The input vector includes variables such as current and historical control inputs, ambient temperature, and gas supply pressure. The output is the measured filling volume. The model assumes that the output equals an unknown function of the inputs plus measurement noise, and that the unknown function follows a Gaussian process distribution.

[0487] A Gaussian process is completely determined by its mean function and covariance function (kernel function). This embodiment uses a composite kernel function, combining the radial basis function (RBF) kernel and the Matérn kernel, to capture variations at different scales. The RBF kernel is suitable for modeling smooth variations, while the Matérn kernel is suitable for modeling variations with a certain degree of roughness. A white noise kernel is used to model measurement noise.

[0488] Composite kernel function:

[0489]

[0490] Radial basis function (RBF) kernel:

[0491]

[0492] in, The signal variance of the RBF kernel controls the overall amplitude of the control function. The length scale of the RBF kernel controls the smoothness of the function and the typical distance of the feature changes; For input samples and The Euclidean distance between them.

[0493] Matérn nucleus (ν=3 / 2):

[0494]

[0495] in , The signal variance of the Matérn kernel controls the amplitude of the function corresponding to the kernel function. The length scale of the Matérn kernel controls the correlation decay rate and local smoothing properties of the control function; The normalized sample distance eliminates the influence of length scale dimensions.

[0496] White noise kernel:

[0497]

[0498] in For the Kronecker delta function, The noise variance measures the intensity of independent Gaussian noise in the observed data. For the Kronecker delta function, when It takes a value of 1 when the condition is met and 0 otherwise, and is used to simulate irrelevant noise terms.

[0499] The kernel function contains several hyperparameters, such as length scale and signal variance. These hyperparameters are optimized by maximizing the marginal likelihood. Marginal likelihood is the probability of the data after integrating the function values, reflecting how well the model fits the data. Hyperparameter optimization uses a quasi-Newton method, iteratively updating the hyperparameters by calculating the gradient of the marginal likelihood.

[0500] With a trained Gaussian process model, the posterior distribution of the output can be predicted for a new input point. The posterior distribution is a Gaussian distribution, determined by the posterior mean and posterior variance. The posterior mean is the best estimate of the output, and the posterior variance quantifies the uncertainty of the estimate.

[0501] System parameters such as flow coefficient, time constant, and valve response delay have a definite mapping relationship with the input-output relationship. The probability distribution of these parameters can be estimated using Gaussian process regression. This embodiment employs variational inference for Bayesian inference. Assuming the posterior distribution of the parameters is Gaussian, the mean and covariance of the posterior distribution are optimized by minimizing the distance between the posterior distribution and the true distribution (KL divergence).

[0502] The optimization employs a stochastic gradient variational Bayesian algorithm, iteratively updating the posterior mean and covariance. After convergence, the posterior mean serves as the most probable estimate of the parameters, and the diagonal elements of the posterior variance represent the uncertainty of each parameter estimate. The smaller the uncertainty, the more reliable the parameter estimate.

[0503] Step S3.2C: Calculate the confidence index of each parameter based on the posterior variance, adjust the excitation signal spectrum corresponding to the parameter whose confidence is lower than the preset threshold, increase the signal energy of the sensitive frequency band of the parameter, and update the Gaussian process regression model after collecting supplementary response data.

[0504] In this embodiment, the excitation signal is adaptively adjusted based on the confidence level of the parameter estimation. The confidence index is defined as a value minus the ratio of the posterior standard deviation to the prior standard deviation. A confidence level close to one indicates that the parameter estimation is reliable, while a confidence level close to zero indicates that it is unreliable.

[0505] For parameters with confidence levels below a preset threshold, the excitation signal needs to be adjusted to improve identification accuracy. Different parameters are sensitive to excitation signals in different frequency bands. The flow coefficient mainly affects steady-state gain and is sensitive to low-frequency signals; the time constant mainly affects dynamic response and is sensitive to mid-frequency signals; the valve response delay mainly affects high-frequency response and phase lag and is sensitive to high-frequency signals.

[0506] The excitation signal adjustment strategy involves increasing the signal energy in the frequency band sensitive to the corresponding parameters. For cases with low confidence in the flow coefficient, the clock period of the pseudo-random sequence is extended to increase low-frequency components. For cases with low confidence in the time constant, the clock period is maintained while the amplitude is increased to enhance mid-frequency components. For cases with low confidence in the valve response delay, the clock period is shortened to increase high-frequency components.

[0507] The specific implementation involves generating a new pseudo-random sequence and weighting its spectrum. The weighting function takes a larger value in the frequency bands that need enhancement and a unit value in other frequency bands. Supplementary response data is collected using the adjusted excitation signal, with a data volume of fifty to one hundred fills.

[0508] New data is added to the training set, and the Gaussian process model and parameter posterior distribution are re-optimized to obtain updated posterior means and covariances. This process is iterated until the confidence index of all parameters exceeds the threshold or the maximum number of iterations is reached. This adaptive strategy can specifically improve the identification accuracy of low-confidence parameters and improve the overall identification quality.

[0509] Step S3.4B: Perform Monte Carlo sampling on the posterior mean and the posterior variance to generate multiple sets of parameter samples. Calculate the stability margin and dynamic performance index of the closed-loop system for each set of parameter samples, and determine the expected value and variance of each optimization objective function.

[0510] In this embodiment, Monte Carlo sampling is performed on the posterior distribution of the parameters to generate multiple sets of parameter samples, and the robust performance of the controller under parameter uncertainty is evaluated.

[0511] Two hundred parameter samples are randomly selected from the posterior distribution of the parameters. Each sample represents a possible system parameter configuration, and the probability of its occurrence is determined by the posterior distribution. For each set of parameter samples, a corresponding system model is constructed, and a closed-loop system is built by combining it with the given candidate solutions for the controller parameters.

[0512] Calculate the stability margin for a closed-loop system, including gain margin and phase margin. Gain margin indicates how much the system gain can be increased while remaining stable, and phase margin indicates how much the system phase can lag while remaining stable. A larger margin indicates stronger system robustness.

[0513] Apply a unit step input to the closed-loop system, calculate the step response, and extract dynamic performance indicators, including overshoot, settling time, and rise time. Overshoot represents the maximum excess of the response, settling time represents the time it takes for the response to reach and remain near its steady-state value, and rise time represents how quickly the response rises from low to high.

[0514] For each set of parameter samples, four objective function values ​​are also calculated: filling accuracy, response speed, control smoothness, and energy consumption. These objectives are calculated through simulation of the filling process.

[0515] Statistically analyze the objective function values ​​for all samples, and calculate the expected value and variance. The expected value represents the average performance of the objective function under parameter uncertainty, while the variance represents the degree of fluctuation of the objective function. Large variance indicates that the performance is sensitive to parameter changes and has poor robustness; small variance indicates stable performance and good robustness.

[0516] A robust optimization objective is constructed by weighting the expected value and variance. The weighting coefficient is called the risk aversion coefficient, ranging from 0.5 to 1.0. A larger coefficient indicates a greater emphasis on robustness, willing to sacrifice average performance for performance stability. The robust optimization objective aims to achieve good average performance while controlling performance fluctuations, ensuring acceptable performance even under parameter uncertainty.

[0517] Step S3.4C: Construct a weight matrix based on the variance values ​​of each objective function. Assign a first update step size to the objective function with a large variance value and a second update step size to the objective function with a small variance value. Use weighted Jacobi iteration to locally refine the candidate solutions generated by the multi-objective genetic algorithm. The second update step size is greater than the first update step size.

[0518] First, calculate the normalized variance of each objective function, which is the variance divided by the square of the expected value, to eliminate the influence of dimensions. Then, construct a diagonal weight matrix, where the diagonal elements are the reciprocals of the normalized variance. Thus, objectives with larger variances have smaller weights, and objectives with smaller variances have larger weights.

[0519] The physical meaning of the weight matrix is ​​as follows: for objectives with high uncertainty (large variance), a smaller step size is used during optimization updates, with careful adjustments to avoid significant performance fluctuations; for objectives with high determinism (small variance), a larger step size can be used to quickly approach the optimum. This adaptive step size allocation strategy can accelerate convergence while ensuring robustness.

[0520] Weighted Jacobian iteration is an iterative method for solving nonlinear equation systems. In multi-objective optimization, the objective function vector is linearized at the current solution to obtain the Jacobian matrix. Each element of the Jacobian matrix is ​​a partial derivative of a certain objective function with respect to a certain decision variable, reflecting the influence of the decision variable on the objective function.

[0521] The Jacobian matrix is ​​approximated using the finite difference method. For each decision variable, a small perturbation is applied near its current value, the change in all objective functions is calculated, and the partial derivative of that column is obtained by dividing the change by the perturbation.

[0522] In standard Jacobian iteration, the inverse of the diagonal portion of the Jacobian matrix is ​​multiplied by the residual vector each time the decision variable is updated. Weighted Jacobian iteration adds a weight matrix to the diagonal portion, adjusting the update step size according to the weights. The expanded update formula is the diagonal portion plus the inverse of the weight matrix, multiplied by the product of the weight matrix and the current solution, minus the product of the transpose of the Jacobian matrix and the objective function value.

[0523] This update method comprehensively considers the gradient information and uncertainty information of the objective function. For high-variance objectives, the weights are small, and the update step size is automatically reduced to avoid large adjustments; for low-variance objectives, the weights are large, and the update step size is relatively large to accelerate convergence.

[0524] The local refinement process involves performing weighted Jacobian iterations on each Pareto front candidate solution generated by the multi-objective genetic algorithm NSGA-II. Each solution is iterated 10 times or until convergence, with the convergence criterion being that the change in solution between two consecutive iterations is less than a threshold. The refined solution is of higher quality and closer to the true local optimum.

[0525] The hybrid optimization strategy combines the global search capability of the NSGA-II genetic algorithm with the local refinement capability of weighted Jacobi iteration. The NSGA-II genetic algorithm first runs for 50 generations to perform a global search, generating an initial Pareto front. Then, local refinement is applied to the higher-quality solutions in the front (the top 20% sorted by crowding distance). The refined solutions are reintroduced into the population, and the genetic algorithm continues to run for 50 generations. This process of local refinement and global search is repeated for 3-5 rounds until convergence.

[0526] Experiments show that, compared with NSGA-II alone, the hybrid method can improve the quality of the solution (hypervolume index) by 35%, increase the convergence speed by 50%, and significantly enhance robustness.

[0527] Step S3.4D: Calculate the performance confidence interval of each solution in the locally refined candidate solution set, filter the candidate solutions by interval dominance relationship, retain robust solutions whose lower bound of performance confidence interval dominates the upper bound of other solution confidence intervals, and select the solution that meets the current production requirements from the robust solutions as the optimal controller parameters.

[0528] In this embodiment, the performance confidence interval of the locally refined candidate solution is calculated, and robust solutions are selected by the interval dominance relationship.

[0529] For each candidate solution, a 95% confidence interval is calculated based on the performance distribution obtained from Monte Carlo sampling. The lower bound of the confidence interval is equal to the expected value minus 1.96 times the standard deviation, and the upper bound is equal to the expected value plus 1.96 times the standard deviation. The confidence interval reflects that, under parameter uncertainty, the performance value has a 95% probability of falling within this interval.

[0530] Define interval dominance: Solution A robustly dominates solution B if and only if the upper bound of the confidence interval of solution A on all objectives is less than the lower bound of the confidence interval of solution B. This means that even considering uncertainty, solution A outperforms solution B in the worst case, and solution A is clearly superior to solution B.

[0531] The robust Pareto front selection process involves checking whether each candidate solution is robustly dominated by other solutions. If it is not robustly dominated by any solution, it is added to the robust solution set. Solutions in the robust solution set are all robustly optimal, maintaining good performance even under parameter uncertainty.

[0532] The final solution is selected from the robust solution set based on production requirements. A demand weight vector is defined, with different weight allocations for different modes. The high-precision mode prioritizes accuracy, the high-speed mode prioritizes speed, the balanced mode assigns similar weights to each objective, and the energy-saving mode prioritizes energy consumption.

[0533] For each solution in the robust solution set, a weighted composite score is calculated. The score is the sum of the products of the normalized values ​​of each objective and their corresponding weights. Normalization maps objectives of different dimensions to the interval between zero and one, making them comparable. The solution with the highest score is selected as the optimal controller parameter.

[0534] Along with outputting the optimal controller parameters, the system also outputs the expected performance and uncertainty information of the solution, including the expected values ​​and confidence intervals of each objective. This information helps operators understand the expected performance and possible fluctuation range of the controller, providing a reference for practical applications.

[0535] By using parameter identification based on Gaussian process regression and multi-objective optimization based on weighted Jacobian iteration, the system can quantify the uncertainty of parameter estimation, actively adjust the excitation signal to improve identification accuracy, consider the impact of parameter uncertainty in multi-objective optimization, accelerate convergence through intelligent weight allocation, and screen out robust controller parameters.

[0536] Practice shows that, compared with traditional methods, the method in this embodiment can improve the parameter identification accuracy to within ±3%, improve the quality of multi-objective optimization solutions by 30-50%, and improve the robustness of the controller under parameter uncertainty by more than 40%.

[0537] Example 8

[0538] This embodiment enhances step S5 based on embodiment 6. The step involves collaboratively optimizing the control parameters of the filling equipment and controlling the filling process using a hierarchical optimization strategy. It introduces a fourth-order tensor synchronization technique based on Tucker decomposition, specifically including the following sub-steps:

[0539] Step S5.1B: Organize the control parameters of the filling equipment into a fourth-order control parameter tensor according to the control type dimension, parameter type dimension, operating condition type dimension, and optimization target dimension.

[0540] In this embodiment, the control parameters of the filling equipment are organized into fourth-order tensors according to four dimensions to achieve a structured representation of the parameters.

[0541] The first dimension is the control type dimension, which includes six control types: pressure control, time control, temperature control, flow control, valve control, and compensation control. Each control type is responsible for one aspect of the system, collectively forming a complete control system.

[0542] The second dimension is the parameter type dimension, which includes eight parameter types: target setpoint, proportional gain, integral time, derivative time, feedforward coefficient, compensation coefficient, filter parameter, and limiting parameter. These parameter types cover all adjustable parameters of the controller.

[0543] The third dimension is the operating condition type dimension, which includes ten operating condition types: high viscosity liquid, medium viscosity liquid, low viscosity liquid, liquid containing particles, high temperature environment, low temperature environment, large volume filling, small volume filling, high-speed production, and precision filling. These operating condition types cover various operating conditions that filling systems may encounter.

[0544] The fourth dimension is the optimization objective dimension, which includes four optimization objectives: accuracy priority, speed priority, smoothness priority, and energy consumption priority. Parameter configurations will differ under different optimization objectives to achieve the corresponding performance requirements.

[0545] Each element of the fourth-order tensor represents the optimal value of a specific parameter under a specific optimization objective, specific operating condition, and specific control type. The total dimension of the tensor is 6 × 8 × 10 × 4, with a total of 1920 elements.

[0546] The tensor filling process extracts known parameter configurations from historical databases and parameter libraries. For operating conditions and target combinations that have been actually run, the corresponding tensor elements are filled with historically optimal parameter values. For unobserved operating condition combinations, the tensor elements are initially empty, and subsequent inference is made using tensor completion techniques.

[0547] To facilitate subsequent decomposition and optimization, the tensors are normalized. For each parameter type, its mean and standard deviation are calculated across all control types, operating conditions, and objectives. Then, the mean is subtracted from all elements of that parameter type, and the result is divided by the standard deviation. Normalization eliminates differences in dimensions and numerical ranges between different parameter types, making them comparable on the same scale.

[0548] Step S5.2B: Perform Tucker decomposition on the fourth-order control parameter tensor, decomposing the fourth-order control parameter tensor into the product of a core tensor and four factor matrices, wherein the dimension of the core tensor is lower than the dimension of the fourth-order control parameter tensor, and each factor matrix represents the main direction of change of the corresponding dimension.

[0549] This embodiment performs Tucker decomposition on the fourth-order control parameter tensor to achieve dimensionality reduction and structured representation of the parameter space.

[0550] Tucker decomposition approximates the original fourth-order tensor as a multilinear product of a low-dimensional core tensor and four factor matrices. The core tensor has a much smaller dimension than the original tensor, and the column vectors of the factor matrices represent the main directions of change for each dimension. This decomposition reveals the inherent low-dimensional structure of high-dimensional data.

[0551] This embodiment employs a high-order singular value decomposition (SVD) algorithm for Tucker decomposition. The algorithm first expands the fourth-order tensor into a matrix along a first mode, where rows correspond to the first dimension and columns correspond to all combinations of the other three dimensions. Singular value decomposition is then performed on the expanded matrix, and the first few left singular vectors are used to form the factor matrix of the first mode. Similarly, singular value decomposition is performed on the expanded matrices of the second, third, and fourth modes to obtain the other three factor matrices.

[0552] With the four factor matrices in place, the core tensor is calculated through a multilinear product of the original tensor and the transposes of the four factor matrices. The dimension of the core tensor is determined by the number of singular vectors preserved in each mode.

[0553] The selection of core tensor dimensions requires a trade-off between compression ratio and reconstruction accuracy. This embodiment employs a multi-criteria selection strategy. The reconstruction error criterion requires that the relative error between the reconstructed tensor and the original tensor does not exceed 10%. The energy proportion criterion requires that the sum of squared singular values ​​retained by each mode accounts for no less than 95% of the total sum. The cross-validation criterion divides known elements into training and validation sets, selecting the dimension combination with the smallest prediction error on the validation set.

[0554] Based on the three criteria, this embodiment selects the core tensor dimension as 3×4×5×2, achieving approximately 60% dimension reduction, from the original 1920 dimensions (6×8×10×4) to 120 dimensions (3×4×5×2), significantly reducing the optimization complexity.

[0555] The column vectors of the factor matrix have explicit physical interpretations. The column vectors of the first mode factor matrix represent the primary control type modes, such as temperature-dominated, pressure-dominated, and hybrid control modes. The column vectors of the second mode factor matrix represent the primary parameter type modes, such as target value-dominated, controller parameter-dominated, compensation coefficient-dominated, and constraint parameter-dominated modes. The column vectors of the third mode factor matrix represent the primary operating condition modes, such as high viscosity characteristic modes, low viscosity characteristic modes, temperature-sensitive modes, capacity-related modes, and speed-related modes. The column vectors of the fourth mode factor matrix represent the primary optimization objective modes, such as quality-first modes and efficiency-first modes.

[0556] The elements of the core tensor represent the interaction strength between different patterns. Larger absolute values ​​of the elements indicate strong interactions between corresponding pattern combinations; these patterns often occur simultaneously and significantly influence each other. Smaller absolute values ​​indicate weaker interactions between corresponding pattern combinations, with less mutual influence. By analyzing the core tensor, key higher-order interaction relationships can be identified, guiding the design of control strategies.

[0557] Step S5.2C: Determine the dimension of the core tensor based on the reconstruction error threshold and the cumulative energy ratio, analyze the column vectors of each factor matrix, and extract the control mode features of each dimension.

[0558] This embodiment analyzes the column vectors of the factor matrix to extract the control mode features of each dimension.

[0559] For each factor matrix column vector, analyze the weight distribution of its elements. The element with the largest weight corresponds to the dominant feature of the pattern. For example, in the first pattern factor matrix, the third element in the first column (temperature control) has the largest weight of 0.86, hence it is named the temperature-dominant pattern. The characteristic description of this pattern is: temperature control is dominant, pressure control and compensation control are auxiliary, suitable for filling temperature-sensitive liquids.

[0560] The correlation between column vectors of different mode factor matrices is calculated to identify frequently occurring mode combinations. The correlation is calculated by dividing the vector inner product by the product of the vector magnitudes, and its value ranges from -1 to +1. A correlation close to one indicates that the two modes are highly correlated and frequently occur simultaneously. For example, the correlation between the temperature-dominant control mode and the temperature-sensitive operating condition mode is 0.82, indicating that a temperature-dominant control strategy is typically used under temperature-sensitive operating conditions.

[0561] A control mode feature library is constructed, recording the feature vector, physical meaning, applicable scenarios, related modes, and recommended parameter configurations for each mode. The feature library provides guidance for practical applications; when the system identifies a certain operating condition, it can query the feature library to find the corresponding control mode and recommended parameters.

[0562] Analyze the core tensor to identify pattern combinations with high interaction intensity. Set the interaction intensity threshold to 0.5; pattern combinations where the absolute value of a core tensor element is greater than the threshold are considered strong interaction combinations. For example, the pattern combination corresponding to a core tensor element of 0.85 is temperature-dominated control, controller parameter-dominated control, high viscosity conditions, and quality-priority objectives—this is a strong interaction combination. The extracted rule is: when prioritizing quality under high viscosity conditions, a temperature-dominated control strategy should be adopted, focusing on optimizing controller parameters.

[0563] Identify anomalous or redundant patterns. If a pattern has very low interaction strength across all core tensor elements, it indicates that the pattern has a weak impact and may be redundant. It can be considered for removal or merging with other patterns to simplify the model.

[0564] Step S5.2D: In the low-dimensional space corresponding to the core tensor, the global control parameters are optimized offline using the Bayesian optimization algorithm to determine the low-dimensional optimal solution. The low-dimensional optimal solution is then reconstructed into the optimized parameters of the original parameter space through the factor matrix.

[0565] This embodiment performs global parameter optimization in the low-dimensional space corresponding to the core tensor, which greatly reduces the optimization complexity.

[0566] The original high-dimensional optimization problem's decision variables are the entire fourth-order tensor, with a dimension of 1920. Through Tucker decomposition, the optimization problem is transformed into optimization within the core tensor space, where the decision variables are the core tensor, with a dimension of only 120. This 94% reduction in the dimensionality of the decision variables significantly reduces the optimization difficulty.

[0567] The objective function is a weighted combination of four performance metrics. For any configuration of the core tensor, it is reconstructed into a tensor of the original parameter space using a factor matrix. Then, historical data simulations or actual tests are performed based on this tensor configuration to evaluate the performance metrics.

[0568] A Bayesian optimization algorithm is employed in the low-dimensional core tensor space. The core tensor is vectorized into a 120-dimensional vector, and a Gaussian process surrogate model is constructed to model the objective function. The Matrn kernel function is used, and 15 initial points are generated using Latin hypercube sampling. The acquisition function employs an upper confidence bound, and the optimization is performed iteratively 80 times.

[0569] Bayesian optimization efficiently searches in a low-dimensional space to find the optimal configuration of the core tensor. Then, it reconstructs the tensor of the original parameter space using a factor matrix, yielding the complete parameter configuration.

[0570] Experiments show that, compared with optimization directly in the original space, dimensionality reduction optimization based on Tucker decomposition reduces computation time by about 85%, increases optimization convergence speed by about 3 times, and improves the quality of the optimal solution by about 15%, while avoiding the local optimum trap in high-dimensional space.

[0571] Step S5.3B: Perform online parameter fine-tuning and feedforward compensation correction on the filling equipment according to the optimized parameters, output the globally optimal filling control scheme, and execute the filling control according to the globally optimal filling control scheme.

[0572] In this embodiment, based on the optimized parameter tensor, combined with the current operating conditions and the target, the corresponding parameter configuration is extracted, and online fine-tuning and feedforward compensation are performed to generate the final control scheme.

[0573] Current operating condition identification is based on real-time sensor data to determine liquid viscosity, ambient temperature, and filling capacity, categorizing the current operating condition into a specific position in the third dimension of the tensor. Optimization objective determination is based on current production needs, determining whether accuracy, speed, smoothness, or energy consumption is prioritized, corresponding to a specific position in the fourth dimension of the tensor.

[0574] Extract parameter values ​​for all control types and parameter types under the corresponding operating conditions and objectives from the optimized tensor to form a preliminary parameter configuration.

[0575] Online fine-tuning calculates the parameter adjustments based on error data from several recent filling cycles. A gradient descent method is employed, with the gradient approximated using finite difference. The updated parameters are then fine-tuned based on the initial configuration to adapt to the current situation.

[0576] Feedforward compensation calculates the compensation amount based on measurable disturbances. The deviations of the current ambient temperature and supply pressure from reference values ​​are measured, and pressure and time compensations are calculated based on a pre-established mapping relationship between disturbances and compensation amounts. The coefficients of the mapping relationship are obtained by fitting historical data.

[0577] The final control parameters equal the extracted parameters plus the online fine-tuning amount plus the feedforward compensation. The system outputs the globally optimal filling control scheme, which includes all necessary control parameters and compensation information.

[0578] According to the control scheme, the filling control system performs specific control actions: setting target values, applying compensation, adjusting the controller based on feedback, driving the actuator, monitoring the results and recording them for the next optimization.

[0579] By using a tensor optimization method based on Tucker decomposition, the system achieves structured representation and dimensionality reduction optimization of high-dimensional control parameter space, automatic extraction and physical interpretation of control modes, collaborative optimization of different dimensions, parameter inference of unobserved conditions, improves computational efficiency by 5-10 times, and improves optimization quality by more than 15%.

[0580] Example 9

[0581] Based on Example 8, this embodiment further enhances multi-timescale coordinated control by performing Tucker decomposition on the fourth-order control parameter tensor, specifically including the following sub-steps:

[0582] Step S5.3C: Divide the control decision of the filling equipment into four time scales according to the rapid control layer, single filling layer, batch adjustment layer and long-term optimization layer, and construct the control decision tensor corresponding to each time scale.

[0583] In this embodiment, the control decisions of the filling equipment are divided into four levels according to the time scale, and control decision tensors are constructed for each level.

[0584] The fast control layer operates at the millisecond level, with an update cycle of 10-100 milliseconds. The decision tensor of this layer contains three dimensions: control variables such as valve opening and pressure regulation, state variables such as current pressure, flow rate, and temperature, and a prediction time domain encompassing the next 10-50 sampling points. The function of this layer is to adjust valves and pressure in real time, responding quickly to disturbances and ensuring a stable filling process.

[0585] The single-fill layer operates on a second-by-second basis, updating once after each fill. This layer's decision tensor contains three dimensions: control parameters such as target pressure, target time, and controller parameters; the filling stage, including acceleration, constant speed deceleration, and stopping; and container types, including different bottle sizes. The purpose of this layer is to set optimal parameters for each fill, adapting to the characteristics of different containers.

[0586] The batch adjustment layer operates on a minute-by-minute basis, updating once per batch or every 100 fills. The decision tensor of this layer contains three dimensions: adaptive parameters such as compensation coefficients and filtering parameters; product batches reflecting the characteristic differences between different batches of liquid; and environmental conditions including temperature and humidity ranges. The role of this layer is to adapt to batch-to-batch differences and environmental changes, maintaining long-term stability.

[0587] The long-term optimization layer operates on an hourly basis, updating every 8 hours or per shift. The decision tensor of this layer comprises three dimensions: control strategies such as algorithm selection and mode switching; optimization objectives such as the weighted allocation of accuracy, speed, and energy consumption; and equipment status such as new equipment, aging equipment, and equipment after maintenance. The role of this layer is global strategy optimization and full lifecycle management of equipment.

[0588] The four time scales constitute a complete hierarchical control system. Each level optimizes different decision variables at different time scales to jointly achieve high-performance control of the filling system.

[0589] Step S5.3D: Define the scale mapping function between control decision tensors of adjacent time scales, and establish the synchronization constraints between control decision tensors of each time scale through the scale mapping function.

[0590] This embodiment defines a scale mapping function between adjacent time scales and establishes synchronization constraints to ensure the coordination and consistency of decisions at all levels.

[0591] The mapping from the rapid control layer to the single-fill layer aggregates the control trajectory of the rapid control layer into the parameter settings of the single-fill layer. Specifically, the target pressure equals the average of all pressure control values ​​during a single filling process in the rapid control layer, and the target time equals the number of rapid control sampling points multiplied by the sampling period. This mapping ensures that the settings of the single-fill layer are consistent with the actual execution of the rapid control layer.

[0592] The mapping from the single-fill layer to the batch adjustment layer uses the parameters of the single-fill layer as a compensation coefficient for the batch adjustment layer. Specifically, the compensation coefficient equals one plus the average deviation rate of the pressure setpoints of all single-fill layers within a batch relative to the nominal values. This mapping reflects the overall offset trend of parameters within a batch and is used for batch-level adaptive adjustment.

[0593] The mapping from the batch adjustment layer to the long-term optimization layer extracts the adaptation parameter trends of the batch adjustment layer as the strategy adjustments for the long-term optimization layer. Specifically, it involves performing trend analysis on the adaptation parameters of multiple batches, such as linear regression, and extracting the slope of the parameters over time, which is then used as the strategy weights for the long-term optimization layer. This mapping captures long-term evolutionary trends and is used for equipment aging compensation and strategy adjustment.

[0594] Synchronization constraints ensure consistency in decision-making across all levels. The consistency constraint requires that the average behavior of the rapid control layer be consistent with the settings of the single-fill layer, with deviations not exceeding the tolerance. The transitivity constraint requires that the composite of multi-level mappings maintain consistency; that is, the rapid control layer, through multi-level mapping to the long-term optimization layer and then back mapping, should produce results close to the original values, with deviations not exceeding the total tolerance.

[0595] These constraints ensure that each timescale level maintains overall coordination while optimizing itself, avoiding conflicts between levels that could lead to a decline in system performance.

[0596] Step S5.3E: The optimization problem with synchronization constraints is decomposed into independent subproblems at each time scale using the alternating direction multiplier method. The controller at each time scale independently optimizes its own objective function, and the coordination between time scales is maintained through the synchronization update step and the Lagrange multiplier update step.

[0597] This embodiment uses the alternating direction multiplier method to decompose the optimization problem with synchronization constraints into independent subproblems at each time scale, thereby achieving distributed collaborative optimization.

[0598] The goal of the optimization problem is to minimize the sum of the objective functions across all time scales, with constraints consisting of synchronization constraints between the levels. This is a large-scale constrained optimization problem; solving it directly is computationally intensive and difficult to implement in a distributed manner.

[0599] The alternating direction multiplier method transforms the constrained optimization problem into the optimization of an augmented Lagrangian function by introducing Lagrange multipliers and penalty terms. The augmented Lagrangian function comprises the original objective function, linear penalty terms for constraint violations, and quadratic penalty terms. The linear penalty terms are weighted by Lagrange multipliers, while the strength of the quadratic penalty terms is controlled by penalty parameters.

[0600] The algorithm's iterative process comprises three steps. The first step involves each timescale level independently optimizing its own decision variables, while fixing the decision variables and Lagrange multipliers at other levels. Each level solves a subproblem with a penalty term; this subproblem only involves the decision variables at its own level, is relatively small in scale, and can be solved efficiently. The optimization at each level can be performed in parallel, significantly improving computational efficiency.

[0601] The second step is to update the Lagrange multipliers. Based on the constraint violations in the current iteration, the Lagrange multipliers are adjusted. If a constraint is violated (left side is not equal to right side), the corresponding Lagrange multiplier is increased by a factor equal to the penalty parameter multiplied by the constraint violation amount. This update rule gradually brings the Lagrange multipliers closer to the optimal dual variable.

[0602] The third step is convergence checking and penalty parameter adjustment. The original residual and dual residual are calculated. The original residual measures the degree of constraint violation, and the dual residual measures the degree of change in the decision variable. If both residuals are less than a preset threshold, the algorithm converges and stops iterating. If the original residual is much larger than the dual residual, it indicates a severe constraint violation, and the penalty parameter is increased to strengthen the constraint; if the dual residual is much larger than the original residual, it indicates a drastic change in the decision variable, and the penalty parameter is decreased to accelerate convergence.

[0603] The distributed implementation of the algorithm deploys controllers at different time scales on different computing nodes. The fast control layer is deployed on real-time controllers such as programmable logic controllers or embedded systems, the single-fill layer is deployed on edge computing nodes, the batch adjustment layer is deployed on local servers, and the long-term optimization layer is deployed on cloud servers. Each node exchanges synchronization information via a lightweight messaging protocol. At the end of each fill cycle, the fast control layer sends the average control value to the single-fill layer; the single-fill layer sends statistical information to the batch adjustment layer every 100 fill cycles; and the batch adjustment layer sends trend data to the long-term optimization layer every eight hours.

[0604] This distributed architecture makes full use of computing resources at all levels, achieves collaborative optimization across multiple time scales, and maintains the system's scalability and fault tolerance.

[0605] Step S5.3F: Identify the first interaction mode and the second interaction mode based on the numerical values ​​of each element in the core tensor, relax the synchronization constraints corresponding to the second interaction mode, and maintain the constraint strength of the synchronization constraints corresponding to the first interaction mode.

[0606] This embodiment identifies strong and weak interaction modes based on the numerical values ​​of each element in the core tensor, and relaxes the synchronization constraints corresponding to the weak interaction mode.

[0607] Analyze the core tensor, calculate the absolute value of each element, and set the interaction strength threshold to 0.6. Elements with absolute values ​​greater than or equal to the threshold correspond to strong interaction modes, while elements with absolute values ​​less than the threshold correspond to weak interaction modes. Strong interaction modes indicate significant mutual influence between the corresponding control mode, parameter mode, operating condition mode, and target mode, requiring strict coordination. Weak interaction modes indicate less mutual influence and can be relatively independent.

[0608] The constraint relaxation strategy is as follows: For synchronization constraints corresponding to strong interaction patterns, the original constraint strength is maintained, and the penalty parameter is set to the baseline value. For synchronization constraints corresponding to weak interaction patterns, relaxation is applied, and the penalty parameter is set to the baseline value multiplied by a relaxation coefficient, with the relaxation coefficient set to 0.3. This reduces the penalty for weak interaction constraints, allowing for greater constraint violations and providing more optimization freedom at each level.

[0609] The adaptive constraint strength is dynamically adjusted based on actual performance. If a violation of a weak interaction constraint causes a performance drop of more than 5%, it indicates that the constraint is actually more important than expected, and the penalty parameter of that constraint is increased, up to the level of a strong interaction constraint. If a strong interaction constraint is too strict, causing optimization difficulties (e.g., the subproblem has no feasible solution or converges extremely slowly), the penalty parameter of that constraint is appropriately relaxed, down to the level of a weak interaction constraint.

[0610] This adaptive strategy provides the system with more flexibility and improves optimization efficiency and performance while ensuring that key constraints are met.

[0611] Step S5.3G: For the missing elements corresponding to unobserved operating conditions in the fourth-order control parameter tensor, perform tensor completion calculation based on the low-rank structure and physical feasibility constraints of the tensor to infer the control parameters under unobserved operating conditions, and use the completed control parameters for the operating condition switching and filling control of the filling equipment.

[0612] This embodiment performs tensor completion calculation on the missing elements corresponding to unobserved operating conditions in the fourth-order control parameter tensor to infer the control parameters under unobserved operating conditions.

[0613] Missing pattern recognition determines which elements in a tensor are observed and which are missing. Observed elements come from historical operational data, while missing elements correspond to combinations of operating conditions that have never been actually run.

[0614] The goal of the low-rank tensor completion optimization problem is to find a low-rank tensor that is as close as possible to the original tensor at observed positions, while satisfying both physical feasibility and smoothness constraints. The low-rank constraint leverages the tensor's inherent low-dimensional structure, assuming that missing elements should maintain a consistent pattern with observed elements. The physical feasibility constraint ensures that the completed parameter values ​​are within a reasonable range, not exceeding physical upper or lower limits. The smoothness constraint requires that the parameter values ​​of adjacent elements (such as adjacent operating conditions or adjacent targets) change smoothly, avoiding unreasonable jumps.

[0615] An alternating least squares algorithm iteratively optimizes the core tensor and factor matrices. Each iteration fixes the factor matrices to optimize the core tensor, then fixes the core tensor and other factor matrices to optimize a specific factor matrix. The optimization process only calculates error at observed locations; missing locations are not included in the error calculation. After each iteration, the reconstructed tensor is projected onto the constraint set, ensuring all elements satisfy the physical feasibility constraints. Iteration continues until convergence, the convergence criterion being that the tensor change between two consecutive iterations is less than a threshold.

[0616] Verification of the completed parameters involves performing simulation validation. Virtual filling is conducted in a simulation environment using the completed parameters to evaluate performance metrics. If the performance metrics are within acceptable limits, the completed parameters are considered reasonable and can be adopted; if the performance metrics are abnormal, they are marked as high-risk parameters, and small-scale testing is recommended before actual use.

[0617] The operating condition switching application identifies the location of a new operating condition in the tensor when it is detected, and checks if there is any observation data at that location. If so, the historically optimal parameters are used directly; otherwise, the inferred parameters are completed using the tensor. The inferred parameters are applied to initiate the modeling process, and performance is monitored in real time. If the error exceeds a threshold, online optimization is triggered to quickly adjust the parameters to adapt to the actual situation. After optimization convergence, the new parameters are added to the tensor as observation data to update and complete the model, allowing it to continuously learn and improve.

[0618] By employing a fourth-order tensor synchronization technique based on Tucker decomposition, the system achieves structured representation and efficient optimization of multi-dimensional control parameters, coordinated consistency of control decisions across multiple time scales, adaptive constraint adjustment of strong and weak interaction modes, parameter inference and rapid adaptation for unobserved operating conditions, and generation and execution of the globally optimal filling control scheme.

[0619] Experimental verification shows that, compared with traditional methods, the parameter optimization calculation time is reduced by 85%, the working condition switching adaptation time is shortened from 50-100 times to 5-10 times, the coordination of multi-timescale control is improved by 60%, and the filling accuracy is stably maintained within ±0.3% under various working conditions.

[0620] Example 10

[0621] This embodiment provides an adaptive control device for filling equipment, such as... Figure 2 As shown, it includes:

[0622] The temperature compensation module is used to collect real-time temperature data of liquid temperature, gas temperature, sensor body temperature and ambient temperature of the filling equipment. It calculates the influence coefficient of each temperature change on the filling volume through a hybrid compensation model that combines physical model and neural network, and outputs the corrected filling control parameters.

[0623] The dynamic time compensation module is used to perform filling and collect filling error data according to the corrected filling control parameters. It performs empirical mode decomposition on the filling error data, separating the filling error into random error components, periodic error components and trend drift components. Based on each error component, it predicts the error trend and calculates the time compensation value.

[0624] The parameter self-tuning module is used to obtain the flow coefficient, time constant and valve response delay parameters of the filling equipment based on the corrected filling control parameters and time compensation value, perform multi-objective optimization calculations on the flow coefficient, time constant and valve response delay parameters, and output the optimal controller parameters.

[0625] The fault diagnosis and fault tolerance module is used to control the operation of the filling equipment according to the optimal controller parameters. It monitors and analyzes the operating status data of the filling equipment in real time through multi-level fault detection, obtains fault diagnosis results, and starts the corresponding fault tolerance control strategy based on the fault diagnosis results.

[0626] The collaborative optimization module is used to collaboratively optimize the control parameters of the filling equipment and control the filling process based on the fault-tolerant control strategy and the operating status data of the filling equipment through a hierarchical optimization strategy.

[0627] The apparatus in this embodiment can be used to execute the technical solutions of the methods described in embodiments 1-9. Its implementation principle and technical effect are similar, and will not be repeated here.

[0628] The above description is only a preferred embodiment of the present invention and is not intended to limit the present invention. Any modifications, equivalent substitutions, and improvements made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

Claims

1. A method of adaptive control of a filling apparatus, characterized in that, include: Real-time temperature data of liquid temperature, gas temperature, sensor body temperature and ambient temperature of filling equipment are collected. The influence coefficient of each temperature change on the filling volume is calculated by a hybrid compensation model that combines physical model and neural network, and the corrected filling control parameters are output. Filling is performed according to the modified filling control parameters and filling error data is collected. Empirical mode decomposition is performed on the filling error data to separate the filling error into random error components, periodic error components and trend drift components. Error trend prediction is performed based on each error component and time compensation value is calculated. Based on the corrected filling control parameters and the time compensation value, the flow coefficient, time constant and valve response delay parameters of the filling equipment are obtained, and multi-objective optimization calculations are performed on the flow coefficient, time constant and valve response delay parameters to output the optimal controller parameters; The filling equipment is controlled to operate according to the optimal controller parameters. The operating status data of the filling equipment is monitored and analyzed in real time through multi-level fault detection to obtain fault diagnosis results. Based on the fault diagnosis results, the corresponding fault-tolerant control strategy is activated. Based on the fault-tolerant control strategy and the operating status data of the filling equipment, the control parameters of the filling equipment are collaboratively optimized and the filling is controlled through a hierarchical optimization strategy.

2. The method of claim 1, wherein, The hybrid compensation model, which combines a physical model with a neural network, calculates the influence coefficient of temperature changes on the filling volume and outputs corrected filling control parameters, including: The real-time temperature data is filtered, and the filtered temperature data is input into the corresponding component of the physical model to calculate the individual influence coefficient of each temperature point on the filling volume. The real-time temperature data is input into the neural network for nonlinear compensation calculation. The neural network takes the liquid temperature, gas temperature, sensor body temperature and ambient temperature as inputs and outputs nonlinear coupling compensation coefficients. Based on the individual influence coefficients of each temperature point and the nonlinear coupling compensation coefficient, a comprehensive compensation coefficient is generated by fusion calculation according to the adaptive weights of each temperature point. The filling control time is then corrected based on the comprehensive compensation coefficient, and the corrected filling control parameters are output.

3. The method of claim 1, wherein, The step of performing empirical mode decomposition on the filling error data, separating the filling error into random error components, periodic error components, and trend drift components, and predicting the error trend and calculating the time compensation value based on each error component includes: The filling error data is subjected to iterative screening. The intrinsic mode function is extracted by identifying extreme points and the mean of the envelope. The filling error data is decomposed into multiple intrinsic mode function components and residual components. According to the frequency characteristics of each component, the intrinsic mode function components and residual components are classified as the random error component, the periodic error component, and the trend drift component. The trend drift component is calculated using a prediction algorithm that combines linear extrapolation and curvature correction. The periodic error component is analyzed in the frequency domain to extract the main frequency and amplitude parameters and calculate the periodic prediction value. The trend prediction value and the periodic prediction value are integrated to determine the comprehensive error prediction value. Based on the comprehensive error prediction value, the filling control adjustment amount is calculated by combining the feedback compensation component and the feedforward compensation component, and the filling control adjustment amount is converted into the time compensation value.

4. The method of claim 1, wherein, The process involves acquiring the flow coefficient, time constant, and valve response delay parameters of the filling equipment, performing multi-objective optimization calculations on these parameters, and outputting optimal controller parameters, including: The current working state of the filling system is determined based on the corrected filling control parameters and the time compensation value. A preset excitation signal is applied to the filling equipment, and the input and output response data of the filling equipment to the excitation signal are collected. The input and output response data are evaluated using a recursive least squares algorithm to estimate parameters and identify the flow coefficient, time constant, and valve response delay parameters of the filling equipment. The first controller parameters are calculated based on the identified flow coefficient, time constant and valve response delay parameters. The first controller parameters are then subjected to rate of change limitation and absolute boundary constraint processing to generate the second controller parameters. Using filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives, a multi-objective genetic algorithm is used to optimize the parameters of the second controller, and the optimal controller parameters are output.

5. The method according to claim 1, characterized in that, The process involves real-time monitoring and analysis of the operating status data of the filling equipment through multi-level fault detection to obtain fault diagnosis results. Based on these results, a corresponding fault-tolerant control strategy is initiated, including: The operating status data of the filling equipment are sequentially subjected to signal-level detection, feature-level detection, and decision-level detection. Signal-level detection is used to perform threshold judgment and trend analysis on the raw sensor signals, feature-level detection is used to extract signal feature parameters and identify anomalies, and decision-level detection is used to fuse and judge multi-sensor data, and output the fault type and fault severity level. The fault diagnosis result is determined by matching the fault type and fault severity level with a preset fault mode library. Based on the fault diagnosis results, the fault-tolerant control strategy is generated by switching to a backup sensor or state estimation mode for sensor faults, activating redundant actuators for actuator faults, and switching to a degraded control mode for control algorithm faults.

6. The method according to claim 1, characterized in that, The method of collaboratively optimizing the control parameters and controlling the filling process of the filling equipment through a hierarchical optimization strategy includes: The current set of optimization variables and constraints are determined based on the fault-tolerant control strategy and the operating status data of the filling equipment. The Bayesian optimization algorithm was used to perform offline global parameter optimization on historical filling data to determine the first optimization parameter; Based on the first optimized parameters, the gradient descent method is used to fine-tune the parameters online in real time according to the preset filling cycle to generate the second optimized parameters; Based on measurable disturbance data, the second optimization parameter is corrected through feedforward compensation calculation, a globally optimal filling control scheme is output, and filling control is executed according to the globally optimal filling control scheme.

7. The method according to claim 4, characterized in that, The input and output response data are evaluated using a recursive least squares algorithm for parameter estimation to identify the flow coefficient, time constant, and valve response delay parameters of the filling equipment. Furthermore, the second controller parameters are optimized using a multi-objective genetic algorithm with filling accuracy, response speed, control smoothness, and energy consumption as optimization objectives, outputting the optimal controller parameters, including: A Gaussian process regression model is constructed, with the control input and observable state variables in the input-output response data as the input of the Gaussian process regression model, and the filling volume measurement value as the output of the Gaussian process regression model. The posterior mean and posterior variance of the flow coefficient, time constant and valve response delay parameter are calculated by Bayesian inference. The confidence index of each parameter is calculated based on the posterior variance. The spectrum of the excitation signal corresponding to the parameter whose confidence is lower than the preset threshold is adjusted to increase the signal energy of the sensitive frequency band of the parameter. The Gaussian process regression model is updated after collecting supplementary response data. Monte Carlo sampling is performed on the posterior mean and the posterior variance to generate multiple sets of parameter samples. The stability margin and dynamic performance index of the closed-loop system are calculated for each set of parameter samples to determine the expected value and variance of each optimization objective function. A weight matrix is ​​constructed based on the variance values ​​of each objective function. A first update step size is assigned to the objective function with a large variance value, and a second update step size is assigned to the objective function with a small variance value. The candidate solutions generated by the multi-objective genetic algorithm are locally refined using weighted Jacobi iteration. The second update step size is larger than the first update step size. For the locally refined candidate solution set, calculate the performance confidence interval of each solution. Filter the candidate solutions by interval dominance relationship, retain robust solutions whose lower bound of the performance confidence interval dominates the upper bound of the confidence intervals of other solutions, and select the solution that meets the current production requirements from the robust solutions as the optimal controller parameters.

8. The method according to claim 1, characterized in that, The method of collaboratively optimizing the control parameters and controlling the filling process of the filling equipment through a hierarchical optimization strategy includes: The control parameters of the filling equipment are organized into a fourth-order control parameter tensor according to the dimensions of control type, parameter type, operating condition type, and optimization objective. The fourth-order control parameter tensor is decomposed into a core tensor and four factor matrices. The core tensor has a lower dimension than the fourth-order control parameter tensor, and each factor matrix represents the main direction of change of the corresponding dimension. The dimensions of the core tensor are determined based on the reconstruction error threshold and the cumulative energy ratio. The column vectors of each factor matrix are analyzed to extract the control mode features of each dimension. In the low-dimensional space corresponding to the core tensor, the global control parameters are optimized offline using the Bayesian optimization algorithm to determine the low-dimensional optimal solution. The low-dimensional optimal solution is then reconstructed into the optimized parameters of the original parameter space through the factor matrix. Based on the optimized parameters, the filling equipment is fine-tuned online and the feedforward compensation is corrected to output a globally optimal filling control scheme, and the filling control is executed according to the globally optimal filling control scheme.

9. The method according to claim 8, characterized in that, After performing Tucker decomposition on the fourth-order control parameter tensor, the method further includes: The control decision of the filling equipment is divided into four time scales: rapid control layer, single filling layer, batch adjustment layer and long-term optimization layer. Control decision tensors corresponding to each time scale are constructed respectively. Define a scaling function between control decision tensors of adjacent time scales, and establish synchronization constraints between control decision tensors of each time scale through the scaling function. The optimization problem with synchronization constraints is decomposed into independent subproblems at each time scale by using the alternating direction multiplier method. The controller at each time scale independently optimizes its own objective function, and the coordination between time scales is maintained through synchronization update steps and Lagrange multiplier update steps. The first interaction mode and the second interaction mode are identified based on the numerical values ​​of each element in the core tensor. The synchronization constraints corresponding to the second interaction mode are relaxed, while the synchronization constraints corresponding to the first interaction mode are kept in strength. For the missing elements corresponding to unobserved operating conditions in the fourth-order control parameter tensor, tensor completion calculation is performed based on the low-rank structure and physical feasibility constraints of the tensor to infer the control parameters under unobserved operating conditions. The completed control parameters are then used for the operating condition switching and filling control of the filling equipment.

10. An adaptive control device for a filling machine, characterized in that, include: The temperature compensation module is used to collect real-time temperature data of the liquid temperature, gas temperature, sensor body temperature and ambient temperature of the filling equipment, calculate the influence coefficient of each temperature change on the filling volume through a hybrid compensation model that combines physical model and neural network, and output the corrected filling control parameters. The dynamic time compensation module is used to perform filling and collect filling error data according to the corrected filling control parameters, perform empirical mode decomposition on the filling error data, separate the filling error into random error components, periodic error components and trend drift components, predict the error trend based on each error component and calculate the time compensation value. The parameter self-tuning module is used to obtain the flow coefficient, time constant and valve response delay parameters of the filling equipment based on the corrected filling control parameters and the time compensation value, perform multi-objective optimization calculations on the flow coefficient, time constant and valve response delay parameters, and output the optimal controller parameters. The fault diagnosis and fault tolerance module is used to control the operation of the filling equipment according to the optimal controller parameters, monitor and analyze the operating status data of the filling equipment in real time through multi-level fault detection, obtain fault diagnosis results, and start the corresponding fault tolerance control strategy based on the fault diagnosis results. The collaborative optimization module is used to perform collaborative optimization and filling control on the control parameters of the filling equipment based on the fault-tolerant control strategy and the operating status data of the filling equipment through a hierarchical optimization strategy.

Images