Indoor temperature real-time regulation and control method of heat distribution pipeline and control system thereof

Abstract

The invention discloses an indoor temperature real-time regulation and control method of a heat distribution pipeline and a control system of the indoor temperature real-time regulation and control method, relates to the technical field of dynamic control of a heat distribution pipe network, and solves the problems of hydraulic oscillation and temperature control hysteresis caused by the fact that local valve regulation neglects whole-network coupling and a first-order linear model is difficult to describe multi-order thermal inertia and large heat capacity in the prior art. According to the scheme, on the basis of pipe network distributed PDE / lumped parameter hybrid modeling and in combination with extended Kalman filtering and unscented Kalman filtering on-line identification, a feedforward decoupling compensation item is generated through spectral decomposition, a self-adaptive multi-model predictive control and iterative learning compensation closed-loop structure is constructed, a control instruction is issued according to a pump-first and valve-second serialization strategy, and a self-adaptive multi-model predictive control and iterative learning compensation closed-loop structure is constructed. Meanwhile, the model weight and the prediction time domain are dynamically adjusted; according to the method, the global balance capability and the temperature tracking precision of heat distribution pipeline regulation and control are remarkably improved.

Description

Technical Field

[0001] The present invention relates to the technical field of dynamic control of thermal pipe networks, and more particularly to a method for real-time control of indoor temperature of a thermal pipe and a control system thereof. Background Art

[0002] Real-time indoor temperature control in centralized heating networks primarily relies on terminal temperature sensors, return water temperature, and network flow measurements. Feedback signals are fed into controllers such as PID, self-correcting fuzzy, or dynamic matrix prediction, driving valve openings and pump speeds to achieve temperature tracking. In recent years, edge-cloud collaboration and artificial intelligence algorithms have also been introduced to achieve a balance between energy conservation and comfort. However, the long lags and strongly coupled multi-node nature of pipeline networks pose greater challenges to control accuracy and real-time performance.

[0003] Existing real-time room temperature control systems often rely on the collaboration of modular hardware and online computing platforms. Patent CN222783576U, for example, installs a temperature-controlled valve and actuator on the household pipes. The valve opening is dynamically adjusted based on feedback from indoor temperature sensors. This solution, combined with heat meters, enables household heat usage monitoring and real-time control. Its advantages lie in its ease of deployment and integration of household metering. Another typical solution, described in CN119642254A, utilizes an edge-cloud collaborative architecture. A heat exchanger station gateway uploads historical meteorological and indoor temperature data in real time. A cloud platform trains a heating control model using algorithms such as LSTM neural networks, then distributes the optimized parameters to the local heat exchanger station for execution. This approach, powered by big data, captures load characteristics and enables rolling iterative optimization. Furthermore, there are a dynamic heat dissipation rate adjustment method based on building thermal parameters (CN119267994A) and a room temperature feedback linear regression control scheme (CN113390126B). These approaches, respectively, utilize simple mathematical models to adjust the water supply temperature based on individual building models and the relationship between water supply and room temperature at the heating station. These technologies are effective in small-scale or simulation scenarios, but they often focus on the refined design of local valve characteristics or single load models. In complex systems, they are exposed to the following shortcomings: For example, patent CN108916984B proposes to manage branch valves in sequence according to heating distance to eliminate horizontal imbalances. However, this local flow distribution scheme ignores the mutual restraint between the upstream pressure of the pipeline network and the partitions. Distributed regulation often causes pump-valve linkage conflicts and hydraulic oscillations, making it difficult for the pipeline network to maintain overall balance under load fluctuations, thereby causing temperature interference and reduced control quality. Secondly, CN113390126B constructs a linear regression relationship between the water supply temperature and the mean room temperature. However, the first-order linear model is difficult to describe the multi-order hysteresis and huge heat capacity characteristics exhibited by long-distance pipelines and floor heating systems. The controller output is seriously mismatched with the actual temperature response, and obvious hysteresis, overshoot, and oscillation often occur, which not only weakens user comfort, but also causes energy waste and aggravates actuator wear. Summary of the Invention

[0004] In view of the deficiencies in the prior art, the present invention discloses a method for real-time control of indoor temperature of a thermal pipeline and a control system thereof, aiming to solve the problems in the background technology.

[0005] In order to achieve the above technical effects, the present invention adopts the following technical solutions: A method for real-time control of indoor temperature of a thermal pipeline, comprising: Step 1: Based on the node temperature, pressure and flow data collected in real time from the pipe network, the PDE / lumped parameter hybrid method is used to construct an initial full-network model including second-order thermal inertia and hydraulic coupling terms. , and output the initial state equations; Step 2: The extended Kalman filter and the unscented Kalman filter are used alternately to identify the water impedance and heat capacity parameters online and obtain the updated model. ; Step 3: Based on the updated model The parameter vector and state equation group of the network are used to calculate the hydraulic coupling coefficient matrix and the thermal coupling coefficient matrix, and generate the corresponding feedforward decoupling compensation terms; Step 4: Based on the feedforward decoupling compensation term and The parameter vector is used to construct a multi-model prediction framework including the long-distance pipeline model and the terminal partition model, and the weighted prediction model group is output through the model pool weight scheduling method. ; Step 5: Use online quadratic programming method according to , the upper and lower bounds of the actuator and the future load forecast curve are used to solve the nonlinear optimization problem including temperature error, energy consumption and actuator constraints to generate the initial control sequence; Step 6: Generate actual pump speed and valve opening control instructions according to the initial control sequence and decoupling compensation terms in accordance with the pump-first-valve-later sequencing strategy and in combination with dead zone and rate limit; Step 7: After the control is executed, the temperature residual and flow oscillation data at the previous moment are used as input to dynamically adjust the multi-model weights and the prediction time domain length, and output the updated control configuration parameters.

[0006] As a further technical solution of the present invention, a real-time indoor temperature control system for a thermal pipeline includes: The data acquisition module is used to collect the original signals of temperature, flow, pressure and valve position of each heat exchange station and the terminal room and align the timestamps according to the unified clock to obtain a standardized time series data set. ; Timing preprocessing module is used to Perform wavelet packet multi-scale filtering and Bayesian interpolation processing to separate noise and fill missing values ​​to obtain a smoothed and denoised time series. ; Parametric modeling module, used to Mapping to the pipe network topology and building a full network model with second-order thermal inertia and hydraulic coupling based on a PDE / lumped parameter hybrid strategy ; Online identification module, used to identify the and As input, the water impedance and node heat capacity parameters are corrected alternately by extended Kalman filtering and unscented Kalman filtering to obtain the updated model ; Coupling and decoupling modules for Perform spectral decomposition on the hydraulic and thermal coupling matrix of the network, extract the subnetwork coupling subspace and generate the feedforward decoupling operator ; Predictive control module for , operator The parallel multi-model is constructed in the rolling time domain with the future load forecast sequence, and the initial control vector is solved by linear quadratic programming. ; Execution monitoring module, used to issue orders according to the sequence strategy of pump first and valve later , and collect room temperature residual and flow deviation series in real time ; Compensation iteration module, used to And the historical residual library is used as input to run the iterative learning algorithm to generate deviation compensation increments And output the corrected control vector ; Closed-loop coordination module, used to 、 The subnetwork interface error is input into the augmented Lagrangian iterative process together to update the model parameters and MPC weights and output the configuration for the next cycle.

[0007] Based on the above technical solutions, the positive and beneficial effects of the present invention are: Through full-network PDE / lumped parameter hybrid modeling and online recursive identification, the hydraulic and thermodynamic coupling characteristics of the pipeline network are captured in real time and dynamically corrected, thereby maintaining a global balance between upstream pressure and flow in each partition when the load fluctuates violently, eliminating pump-valve linkage conflicts and hydraulic oscillations caused by local valve adjustment.

[0008] With the help of the synergistic effect of the feedforward decoupling compensation term generated by spectral decomposition and the incremental multi-model predictive control, multi-order thermal inertia and huge heat capacity effects are synchronously compensated in the prediction time domain, so that the hysteresis and overshoot of the temperature response can be fundamentally suppressed.

[0009] Under the premise of considering the upper and lower bounds of the actuator and the rate limit, the online quadratic programming integrates the temperature error, energy consumption and actuator constraints into the optimization framework, so that the control sequence is highly matched with the actual dynamics of the pipeline network, thereby reducing the mismatch between the controller output and the field response.

[0010] Adaptive weight scheduling and prediction time domain adjustment based on residual monitoring form a closed-loop feedback. The pipeline network can continuously update the model and control parameters under different working conditions, thereby maintaining stable temperature control in complex systems and effectively delaying actuator wear. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] In order to more clearly illustrate the embodiments of the present invention or the technical solutions in the prior art, the following briefly introduces the drawings required for use in the embodiments or the description of the prior art. Obviously, the drawings described below are only some embodiments of the present invention. Those skilled in the art can also derive other drawings based on these drawings without inventive work, among which: Figure 1 This is a structural diagram of a method for real-time control of indoor temperature of a thermal pipeline according to the present invention; Figure 2 This is a flowchart of the working method of step 1 of the present invention; Figure 3 This is a step diagram of the calculation method of step 3 of the present invention; Figure 4 Schematic diagram of indoor temperature error and control signal change of the present invention; Figure 5 This is a schematic diagram of the real-time indoor temperature control system of the thermal pipeline of the present invention. DETAILED DESCRIPTION

[0012] To make the purpose, technical solutions, and advantages of the embodiments of this application more clear, the technical solutions of this application will be clearly and completely described below in conjunction with the embodiments. Obviously, the embodiments described are part of the embodiments of this application, not all of them. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of this application.

[0013] To facilitate understanding of this embodiment, firstly, a method for real-time control of indoor temperature of a thermal pipe disclosed in the embodiment of this application is described in detail. Figure 1 The figure shows a schematic diagram of the steps of a method for real-time control of indoor temperature of a thermal pipeline, which includes: Step 1: Based on the node temperature, pressure and flow data collected in real time from the pipe network, the PDE / lumped parameter hybrid method is used to construct an initial full-network model including second-order thermal inertia and hydraulic coupling terms. , and output the initial state equations; specifically, Figure 2 As shown, the system relies on a network of three high-precision sensors to continuously monitor the operational status of the pipeline network. First, pressure sensors are deployed along key nodes such as the main trunk line. Typical models include piezoresistive or capacitive sensors with a range of 0–2 MPa and an accuracy of 0.1% FS. The sampling frequency is generally between 10 Hz and 0.1 Hz, which can be flexibly adjusted to meet system response requirements. The specific sampling frequency is determined based on actual conditions and is not specified. The pressure gradient distribution within the region is obtained by comparing the pressure difference between adjacent measuring points, which is used to determine the fluid drive and damping characteristics. Second, flow meters are installed at branch outlets or bifurcations, typically using vortex or electromagnetic solutions. When the flow pulsation amplitude per unit time falls below the turbulence threshold (generally taken as the lower velocity limit corresponding to Re = 2300, approximately 0.2 m / s; the flow threshold value can be determined based on pipe diameter and temperature conditions and is not specified), the branch is dynamically classified into the lumped parameter domain. Third, distributed temperature sensing uses fiber Bragg grating (FBG) or distributed optical scattering (DAS) technology, with measurement points set every 0.5 to 5 meters along the pipeline to capture the axial temperature gradient field. The resolution and spacing of the measurement points can be selected based on the insulation thickness and the expected temperature gradient sensitivity. The specific resolution and spacing can be determined based on actual conditions and are not limited. Next, the raw signal needs to be denoised through multi-stage filtering. Wavelet threshold denoising is preferably used for pressure and flow signals to eliminate high-frequency interference. Temperature data can be smoothed using sliding average or Kalman filtering. All data must be clock-aligned. PTP (Precision Time Protocol) or IEEE1588 precise synchronization protocol is recommended to achieve microsecond timestamp consistency to ensure accurate correspondence of spatiotemporal information in subsequent hybrid models.

[0014] After data acquisition is complete, the central processing unit reads the continuous time-series pressure gradient distribution sequence and introduces the Morris sensitivity analysis algorithm to dynamically partition the network's hydraulic domain. Morris sensitivity analysis, a factor-by-factor screening method based on parameter perturbations, calculates the sensitivity index of pressure gradient changes by sequentially perturbing hydraulic parameters such as pipe diameter, length, and flow resistance coefficient for each pipe segment (for example, with a relative perturbation amplitude set to ±10%). If the sensitivity index of a pipe segment or node group exceeds a preset threshold, the segment's hydraulic dynamics are considered to have a significant impact on the overall network coupling and should be classified into the PDE domain; otherwise, it is classified into the lumped parameter domain. The threshold value S_{thresh} can be adjusted based on the network size, pipe diameter distribution, and load fluctuation amplitude. The specific value can be determined based on actual conditions and is not limited.

[0015] At the hydraulic solution level, the core of the implementation plan using the PDE / lumped parameter hybrid domain method lies in dynamic domain division. First, a Morris global sensitivity analysis is performed on the geometric parameters (inner diameter D, length L, roughness ε) and current flow velocity u and other input variables of each pipe section in the pipeline network to evaluate the contribution of the pipe section to the local pressure gradient change. When the sensitivity index is higher than the preset threshold (such as the normalized sensitivity index exceeds 0.1, it can be flexibly set in combination with the total number of nodes in the pipeline network and the expected amount of calculation, and can be determined according to the actual situation without limitation), the section is dynamically included in the PDE solution domain; otherwise, it is considered to have limited impact on the overall pressure distribution and is included in the Lumped domain. The second-order thermal inertia term is reflected in the subsequent coupled state equation as an additional term related to heat capacity and flow velocity acceleration, but it does not participate in the judgment criteria in the domain division stage. The main numerical solution in the PDE domain relies on the characteristic line method to solve the one-dimensional transient flow conservation equation, which is simplified as follows: in represents the cross-sectional area of ​​the pipe, u is the flow velocity, h is the head height, g is the acceleration due to gravity, f is the friction factor, and D is the pipe diameter. The characteristic line method maintains the physical integrity of information transmission by integrating along the curve in both the positive and negative characteristic directions. Simultaneously, only a few point values ​​along the characteristic line need to be stored, enabling high-precision reconstruction of the main flow within a time scale of seconds. To ensure solution stability, the time step Δt and the spatial step Δx must satisfy the CFL condition: In addition, when high-frequency oscillations or numerical fluctuations occur in the PDE domain, the TVB (Total Variation Bounded) limiter can be enabled to control the fluctuations. However, this strategy and the CFL condition can also be adjusted according to the real-time response requirements of the system without any restrictions.

[0016] The lumped domain uses the node method to construct the equivalent resistor-capacitor network model, and its basic equation is: p represents the equivalent node pressure, Q is the flow input, R is the hydraulic impedance, which corresponds to the pressure drop ratio caused by the coal, water quality and wall state of the pipeline; C is the capacity coefficient, which reflects the momentum storage characteristics of the water body and the pipe wall. R and C can be obtained through historical calibration conditions or on-site The model is obtained through experimental calibration and can also be used for online identification using field sampling data in the early stages of system operation. This lumped model only needs to process the flow and pressure relationship between nodes, without the need for spatial discretization, and the computational complexity is much lower than that of PDE solution.

[0017] In order to achieve close coupling of the hybrid domain, the system monitors the pressure p_b and flow Q_b at the boundary nodes of the PDE domain and the Lumped domain at the same time, and sets the synchronization trigger condition: when When the boundary variables are immediately injected into the contralateral model through the data interface for correction. and The value of can be determined based on the maximum operating pressure level of the pipeline network and the branch flow scale, for example The pressure can be set between 0.005 MPa and 0.02 MPa, and ΔQ_sync can be set between 2 L / min and 10 L / min. The specific range is determined based on actual conditions and is not limited. If the boundary deviation persists for more than N_sync time steps (generally 3 to 5 steps), the system triggers a re-division mechanism, reassesses the sensitivity of the pipe segment or region, and may adjust its domain affiliation to ensure that the model is consistent with the actual state.

[0018] The thermal modeling module uses the axial temperature field T(x,t) obtained by distributed temperature sensing as input and uses Karhunen–Loève decomposition to extract its dominant spatial mode. Specifically, the temperature field is regarded as a random process and its spatial covariance function is solved. The Fredholm integral equation of , we get the characteristic function and the corresponding eigenvalues , and then reconstruct the high-dimensional temperature field as: in = , r is the mode cutoff number. Truncation ceases when the cumulative energy of the r-order preceding mode reaches 95% or above. This can also be adjusted based on the pipe's insulation performance and expected temperature response time. The specific value is determined based on actual conditions and is not limited. After modeling, thermal dynamics are compressed into an r-dimensional state space, significantly reducing computational pressure.

[0019] The online parameter correction relies on the instantaneous heat flux q_w(t) measured by the wall heat flow meter, which is consistent with the heat flux predicted by the model. For comparison. When , the parameter update subroutine is triggered, and the wall specific heat capacity C_w and heat transfer coefficient h are iteratively optimized through extended Kalman filtering or recursive least squares method. It can be set between 5W / m² and 20W / m². The specific value varies with the pipe diameter and insulation thickness. It can be determined according to the actual situation and is not limited. If the deviation exceeds the limit for M_sync consecutive sampling periods (usually M_sync = 3 to 5), the filter convergence step size can be increased or the suboptimal model update mode can be switched to speed up the adjustment.

[0020] In terms of hardware and software architecture, the PDE domain dynamic partitioning and solver are deployed on edge computing units, leveraging multi-core CPUs or FPGAs to achieve millisecond-level responses. KL decomposition and online parameter identification run on central servers or the cloud, leveraging parallel computing resources for efficient processing. Modules exchange data via real-time industrial bus protocols such as OPC UA, MQTT, or DDS, connecting to on-site DCS or SCADA systems. Under real-world operating conditions, whether experiencing transient flow changes in the trunk pipeline or sudden disturbances in heating and cooling loads, the system can complete model updates and state estimation within multiple control cycles (typically no more than ten seconds), providing precise support for subsequent model predictive control (MPC) or optimized scheduling.

[0021] Step 2: The extended Kalman filter and the unscented Kalman filter are used alternately to identify the water impedance and heat capacity parameters online and obtain the updated model. ; In the process of online identification of water impedance and node heat capacity parameters, the system first uses the coupled water-heat state model generated in step one and the real-time collected pressure, flow and temperature data to build a complete state-observation mapping relationship. After the sensor network is low-pass filtered and Kalman smoothed, the pressure difference, branch flow pulsation and distributed temperature modal coefficients are interpolated and synchronized in the form of equally spaced time windows (for example, approximately five to ten control cycles) to form a smooth observation sequence. The model prediction module uses the last updated parameter value as the initial input to drive the water-heat numerical simulation and generate a prediction output that matches the current observation time step. After subtracting the two, a high-dimensional residual vector can be obtained. The covariance matrix of this vector in the time window and its Frobenius norm are used as the main basis for judging the degree of nonlinearity of the system.

[0022] When the residual norm falls within the pre-set linear approximation interval (typically between 0.1 and 0.2 in normalized dimension), indicating that the deviation between the current parameter estimate and the field state is within the acceptable error band of first-order Taylor linearization, the system automatically selects the extended Kalman filter (EKF) for rapid iteration. The nonlinear water-thermal model then calculates local sensitivities in the state variable and parameter directions using perturbation or numerical difference methods. This method predicts the rate of change in the output for a 1% increase in each water impedance or heat capacity parameter at the current estimated point. The resulting finite-difference sensitivity matrix is ​​equivalent to a numerical approximation of the Jacobian matrix, avoiding the high cost of manual derivation of complex models while meeting the requirements of high-precision linearization. Based on this sensitivity matrix, the filter can construct the Kalman gain and perform parameter correction with minimal computational effort (typically equivalent to several hundred multiplication-addition operations) after performing a single state prediction. After the parameter update is completed, the covariance matching reconstruction of the noise matrix is ​​performed according to the set forgetting factor (such as 5%) using the deviation between the residual and the predicted covariance to ensure that the EKF adaptively tracks the statistical characteristics of the process and observation noise.

[0023] If the residual norm exceeds this linear interval, it indicates significant model nonlinearity, and simple first-order linearization will result in excessive error. In this case, the system seamlessly switches to the unscented Kalman filter (UKF). The UKF does not rely on explicit derivatives. Instead, it constructs a set of weighted sigma points based on the current state estimate covariance, using a sampling point generation strategy where the state dimension n is multiplied by 2 plus one. This sampling point generation process uses the Cholesky decomposition as its numerical basis, and the distribution and concentration of the sigma points in the state space are controlled by adjusting the parameters α and κ (e.g., α is 0.1 and κ is 3−n). Each sigma point is then input into the complete nonlinear simulation model to obtain the corresponding predicted and observed sigma points. The prior mean and covariance are then synthesized based on preset weights. Finally, the real-time residuals are combined to calculate the UKF gain and calibrate the state parameters. The UKF is inherently robust to strong nonlinear scenarios, maintaining stable estimation accuracy despite large perturbations or sudden changes in model structure.

[0024] After each EKF or UKF parameter update, the system immediately employs a posterior smoothing strategy based on Rauch–Tung–Striebel to reversely fuse the latest estimate with the previous one to eliminate potential jumps in estimates caused by filter switching or nonlinearity. Posterior smoothing uses the sensitivity matrix obtained from the previous EKF step or the Sigma propagation information from the UKF to calculate smoothing gains and make weighted corrections to the forward estimate, ultimately outputting a smoothed and continuous time series of water impedance and heat capacity parameters. This process completes in just a few milliseconds on the edge computing node.

[0025] At the same time, online noise reconstruction runs in parallel after each filtering iteration. By blending the outer product of the prediction residuals for the current cycle with the current noise matrix at a set ratio (e.g., 5% vs. 95%), the system continuously updates the covariance estimate of process and observation noise, enabling the filter to adaptively respond to uncertainties such as heat source load fluctuations, pipeline wear, and changes in the external environment. If the statistical characteristics of the field noise suddenly change, such as a sudden increase in flowmeter signal noise, the reconstruction ratio can be dynamically increased to accelerate estimation convergence. The specific reconstruction rate can be flexibly set within a range of 1% to 10%, with no hard limit.

[0026] After the EKF and UKF run alternately and complete smoothing, the latest water impedance parameters and heat capacity parameters are brought back to the water-heat coupling simulation model in step one to predict key indicators such as the terminal return water temperature and the maximum main pressure difference. If the maximum relative error between these predicted results and the real-time observation exceeds the pre-set tolerance (such as the pressure tolerance is 1% of the full scale and the temperature tolerance is 0.5K), the system determines that the current online parameter identification has lost convergence or the model has drifted, and then triggers the model reconstruction process. Model reconstruction first reinitializes the noise covariance to the design calibration value and returns the estimated parameter to the recent smooth average level. If convergence is still not possible, global optimization methods such as particle swarm or genetic algorithm are run in parallel on the central server. Through multiple rounds of iterations, the optimization is sought in the parameter space, and the new parameters are finally sent to the edge node to complete the so-called model reconstruction action in the manual.

[0027] In the above process, it should be noted that the residual norm criterion is used for filter switching; the perturbation ratio provides a balance between accuracy and numerical stability for the numerical approximation of sensitivity; the forgetting factor determines the adaptive rate of noise reconstruction; the smoothing gain ensures that the posterior smoothing achieves a compromise between debounce and retaining dynamic response; and the tolerance of the consistency detection corresponds to the actual control requirements. These values ​​can be adjusted according to the scale of the pipeline network, the amplitude of flow and temperature difference fluctuations, and the performance of the edge computing hardware. The specific values ​​can be determined by the site environment and project requirements without rigid restrictions. Through the alternating execution of EKF and UKF, the organic combination of noise adaptive reconstruction, posterior smoothing and model consistency detection, the system realizes real-time, accurate and robust online identification of water impedance and node heat capacity parameters, providing a reliable and continuously self-correcting model support for subsequent model predictive control and optimized scheduling.

[0028] Step 3: Based on the updated model The parameter vector and state equation group of the network are used to calculate the network hydraulic coupling coefficient matrix and thermal coupling coefficient matrix, and generate the corresponding feedforward decoupling compensation terms; Figure 3 As shown, specifically, the calculation steps of step 3 include: s100, based on the water impedance parameters and node heat capacity parameters in the updated model output in step 2, construct the hydraulic coupling coefficient matrix through the node admittance matrix analysis method ,in For nodes The pressure difference, For branch The node admittance matrix analysis method solves the sensitivity of the pressure gradient of the pipeline network node to the flow change through the Jacobi iteration method, and outputs a diagonal-dominated sparse matrix ; s200, according to the state equations of the heat capacity chain state space model, the state matrix Perform singular value decomposition and retain the first r significant singular values ​​to reconstruct the reduced-order matrix , extract the partial derivative of heat flux to stratification temperature and generate the thermal coupling coefficient matrix ,in is the temperature state of the kth layer, is the heat flux; s300, calculate hydraulic coupling matrix The off-diagonal norm of If it is greater than the preset threshold , calculate the pressure disturbance suppression term, and then generate the valve opening correction through the feedforward compensator. The formula is: (1) In formula (1), is the gain matrix, which is used to map the pressure fluctuation into the valve opening correction; Real-time pressure gradient residual, in Pa / m, used to reflect the uncompensated oscillation energy; s400, synchronous calculation of thermal coupling matrix condition number ,exist Greater than the preset threshold When , the temperature lag compensation term is generated based on the pseudo-inverse matrix of the heat flux temperature transfer function through the thermal inertia decoupling algorithm. The formula is: (2) In formula (2), is the pseudo-inverse matrix of the heat flux-temperature transfer function, which is used to offset the multi-order lag effect of the building heat capacity; is the target heat flux deviation; s500, based on the real-time pipeline network pressure oscillation spectrum energy , and thermal inertia residual Ratio , weighted fusion of compensation terms generates the final feedforward decoupling instruction In the scenario of real-time temperature control of indoor thermal network, step 3 requires to build the hydraulic coupling coefficient matrix of the network based on the model M1 output in step 2 and the updated water impedance parameters and node heat capacity parameters. and the thermal coupling coefficient matrix , and based on these two, a feedforward decoupling compensation term is generated to address the uncertainty and cross-interference caused by the water-heat synergistic coupling. The entire process must be completed within each control cycle (typically 1 to 5 seconds) to meet the real-time requirements of indoor temperature control. The numerical thresholds and matrix truncation orders must balance computing resources and comprehensive control accuracy. The specific determination can be made based on the scale of the on-site pipeline network, the amplitude of load fluctuations, and the edge computing capabilities, and is not limited.

[0029] Hydraulic coupling coefficient matrix The calculation of is based on the inverse operation of the node admittance matrix. , essentially reflects the hydraulic admittance between any two nodes in the pipeline network in the frequency domain or static state, and is constructed through online updated water impedance parameters; its inverse matrix corresponds to the sensitivity mapping of the pressure difference of the pipeline network nodes to the flow rate. System call After obtaining the preliminary coupling matrix, in order to characterize the more realistic sensitivity under transient conditions, the Jacobi iteration method is used to fine-tune the derivative of the pressure gradient to the flow change. The specific steps are: for each branch , based on the current traffic Add a small disturbance (For example, take 5L / h or 1%FS), use the fast hydraulic simulation module or the table lookup method to calculate the pressure difference increment of each node ; and feed the disturbance-response pair back to the initial point in the iteration until the change converges to the preset accuracy (for example, the error is less than 0.01Pa under a 0.1Pa change), forming a numerically approximate first-order partial derivative This fills the matrix . It should be noted that the "iterative convergence accuracy" in this application is different from the 0.1% error standard commonly used in offline simulation. It can be adjusted within the range of 0.01Pa~1Pa according to the pressure level of the pipeline network and the control accuracy requirements. The specific accuracy can be determined by the on-site measurement accuracy and is not limited. After completing the above steps, the obtained AHA_HAH is usually a sparse matrix dominated by row diagonals, and its diagonal elements dominate the "self-flow-self-pressure difference" coupling of each node, and the remaining non-diagonal elements reflect the auxiliary coupling between nodes.

[0030] Thermal coupling coefficient matrix The extraction of heat capacity chain state space model relies on this model. This model can be divided into layered or partitioned temperature states. (e.g. multi-point temperature measurement data from the heat exchange station outlet, terminal heat dissipation surface and return water outlet) and heat flux input (The product of the flow rate of hot water or steam and the temperature difference) is described by a set of linear ordinary differential equations: Among them, the matrix Contains the thermal inertia and thermal diffusion characteristics of each section of the pipeline and the indoor heat exchange equipment, the matrix Describe the direct driving effect of heat flux input on state change. In order to reduce the dimension and extract the main coupling mode, the system Perform singular value decomposition: Keep the first r singular values (The energy contribution rate corresponding to these values ​​reaches about 95% cumulatively), and the reduced-order matrix is ​​reconstructed as follows: Get the corresponding column.

[0031] The thermomechanical coupling coefficient matrix is ​​then defined as the instantaneous partial derivatives of each input flux with respect to each temperature state: It should be noted that the "singular value energy contribution rate" in this application is different from the variance explanation rate of traditional principal component analysis. It takes into account the impact of thermal hysteresis dynamics on spatial coupling; singular value decomposition can be based on a direct QR algorithm or a randomized SVD acceleration scheme, which can be determined by the number of nodes m and computing power without limitation.

[0032] After the coupling matrix is ​​constructed, the system enters the stage of generating the feedforward decoupling compensation term. First, for the hydraulic disturbance suppression, it is necessary to measure the off-diagonal energy of the coupling matrix - that is, the remaining matrix after excluding the diagonal elements. The Frobenius norm of : , when this norm exceeds the preset threshold (For example, a dimensionless value of 0.1 to 0.5 is used to evaluate the degree of cross-node coupling. The specific value can be set according to the scale of the pipeline network and the coupling strength, and there is no limit.) This indicates that the flow-pressure difference coupling between nodes is significant and needs to be suppressed through feedforward control. The control quantity is based on the formula Generate, where the vector Represents the transient pressure gradient residual of each node (unit: Pa / m, which can be obtained by dividing the difference between the readings of two adjacent pressure sensors by the installation spacing), the matrix is the feedforward gain matrix, which is usually calculated by engineers based on The inverse or pseudo-inverse of is combined with the proportional adjustment coefficient (the value ranges from 0.5 to 2.0 and can be flexibly set according to the regulator response speed and system stability requirements, without limitation) to determine: Calculation completed It represents the correction amount required for each valve opening or pump frequency. During specific operation, the controller will write these instructions to the target register of each valve position controller or frequency converter through MODBUS or OPC UA protocol. The actuator completes the opening or pump speed adjustment at a fixed ramp rate (such as 0.1% / s) to avoid mechanical vibration.

[0033] The triggering conditions of thermal decoupling compensation are based on the matrix condition number when >ϵ2 (usually 10³~10 4 , to measure the numerical pathological degree of the thermal coupling matrix, which can be set according to the simulation test results and is not limited), it indicates that there is a high-order lag or amplification effect in the coupling of heat flux to temperature, and it is necessary to start the thermal inertia decoupling algorithm. This algorithm relies on the pseudo-inverse matrix of the temperature transfer function To offset the multi-order lag, the formula is ,(2) where the vector Indicates the target deviation of the heat load of the heat exchange station or terminal (unit can be kW or m³ / h), the matrix ∈ Each row of is usually obtained by system identification or empirical model, reflecting the multi-order delay relationship from heat flux input to temperature response. Converted into the water supply temperature increment or heat pump load adjustment instruction of the corresponding heat exchange station, and sent through the building automation BACnet or Modbus interface.

[0034] In order to balance the hydraulic and thermal compensation quantities, the system adopts dynamic weighted fusion: .in To target the main pressure signal of the pipeline network in the oscillation frequency band The integral of the power spectral density over the frequency band (e.g. 0.1Hz2Hz) represents the hydraulic suppression demand; it is defined as the norm of the difference between the latest temperature prediction and the actual measurement (unit K); when ≫ Thermal inertia residual, If the value is close to 1, the system prioritizes hydraulic suppression; if it is lower, it prioritizes thermal compensation. It should be noted that "spectral energy integration" in this application is distinct from instantaneous energy measurement and must be implemented using a real-time FFT or bandpass filter and energy detection module. The calculation window length can be selected between 5s and 30s, and can be determined based on the dynamic characteristics of the pipeline network and the control cycle synchronization strategy. There is no limitation on this.

[0035] The final fused feedforward decoupling instruction Contains valve opening or pump speed offset and heat exchange water temperature difference or power offset , sent to each actuator through the industrial bus in the form of a continuous signal, to achieve water-heat collaborative decoupling control of the indoor thermal pipe network. The implementation of the entire step 3 not only covers the complete link from model data acquisition, sensitivity matrix calculation to feedforward compensation generation and dynamic weight fusion, but also provides specific instructions for various trigger conditions, numerical thresholds, matrix operation methods and execution interfaces, and provides flexible value selection space with "specific values ​​can be determined according to actual conditions without limitation" where necessary, to ensure that the solution has good applicability and engineering feasibility in indoor thermal pipe systems of different scales and topologies. It should be noted that the "node admittance matrix analysis method" in this application differs from the existing static hydraulic admittance method in that it can obtain the pressure-flow coupling sensitivity between pipeline nodes through comparative iterative calculations under real-time operating conditions, and then output the mapping coefficients required for hydraulic decoupling. The "heat capacity chain state space model" differs from the traditional single-pipe thermal balance model in that it can uniformly describe the thermal inertia of multiple pipeline sections and heat exchangers in the form of a state matrix, facilitating the extraction of major coupling modes through singular value decomposition. The "thermal inertia decoupling algorithm" differs from conventional PID or first-order lag compensation mechanisms in that it can actively generate multi-order thermal lag compensation based on the inverted pseudo-inverse matrix of the temperature transfer function. The "pressure oscillation spectrum energy" differs from the instantaneous pressure measurement value in that it is the integral of the power spectral density within a specific oscillation frequency band after frequency domain analysis of the pressure signal at the main or branch outlet of the pipeline network, which can quantify the intensity of system oscillations. The "dynamic weighted relaxation" differs from the static weighting coefficient in control theory in that it is adjusted online based on the ratio of the real-time measured hydraulic oscillation energy to the thermal inertia residual, ensuring that feedforward compensation can suppress oscillations without losing temperature response.

[0036] Step 4: Based on the feedforward decoupling compensation term and The parameter vector is used to construct a multi-model prediction framework including the long-distance pipeline model and the terminal partition model, and the weighted prediction model group is output through the model pool weight scheduling method. In the scenario of real-time temperature control of indoor thermal network, the implementation of step 4 first depends on the key operating condition indicators collected online by the system: the sudden change rate of flow in any trunk or branch of the network; (unit can be m³ / s², when the value exceeds 2m³ / s², it is determined to be a severe flow disturbance) and the spatial variance of the terminal temperature gradient (Measured in °C² / m², a value exceeding 0.5°C² / m² indicates significant temperature unevenness). The operating condition feature classifier monitors the above two indicators. When the flow rate suddenly changes, , the system automatically determines that it is a "hydraulic-dominated working condition" and calls the long-distance pipeline model based on transmission line theory. This model uses the one-dimensional wave equation: As the core, P(x, t) and Q(x, t) are the coordinates along the pipeline respectively. Pressure and flow on is the fluid density, is the wave velocity (can be calibrated online according to the pipe material, filling fluid and pressure level). In actual operation, the system will generate the feedforward decoupling compensation value in step 3. Feedback is fed back to the model's boundary condition module to adjust the model's initial pressure and flow boundary values, ensuring that the model predictions reflect the compensation action immediately. At the moment of the operating condition switch, the model adds the corresponding compensation amount to the original boundary pressure or flow value and outputs the predicted pressure waveform in the next cycle through parallel calculations such as finite differences or characteristic line methods.

[0037] If the terminal temperature gradient variance ) exceeds the 0.5°C² / m² threshold, the system determines it as a "thermal dominant condition" and switches to the terminal partition model. This model is based on the distributed heat capacity chain state equation where index and Traverse the same area or adjacent partitions, For the Node equivalent heat capacity, is the equivalent thermal resistance between nodes, both of which can be obtained from the online identification results in step 3 or the initial calibration data. The system also uses the feedforward compensation The model is superimposed on the heat exchanger output temperature or heat flow boundary, driving the model state to converge to the new operating value. The model uses implicit Euler or Gear methods to ensure numerical stability.

[0038] Regardless of which dominant model is used, the system will retain the prediction results of the other model in parallel for residual comparison and dynamic weighting. Generate residuals When , the normalized residual amplitude is calculated. In order to achieve online update of model weights, the system uses recursive least squares method with forgetting factor to adjust the weights , where the forgetting factor λ∈(0,1) can be set to typically 0.01~0.1 to balance the influence of new and old residuals. It can be set according to the dynamic frequency of the pipeline network and the measurement noise level without limitation. In this way, the model weight The model that best suits the current working conditions will be automatically selected based on the size of the residual, ensuring smooth prediction and controlled errors when the working conditions are switched or coexisting.

[0039] Ultimately, the prediction outputs from multiple models are fused in a weighted sum. This fusion result is then used as input for model predictive control (MPC) to optimize the heat source's outlet water temperature, pump speed, and valve opening. During action execution, the system simultaneously sends the weighted, fused control increments to the corresponding execution units. The system then verifies command arrival and execution feedback via the SCADA or BACnet network. If a discrepancy is detected between the command delivery and execution feedback, the system automatically triggers a compensation fault-tolerance mechanism, re-adjusting weights or amplifying the compensation amount.

[0040] The "condition-adaptive model pool" in this application is different from the traditional single model or static switching mechanism. It can trigger the optimal model online according to the real-time flow mutation and temperature gradient characteristics, and maintain the continuity of prediction and control through weighted fusion; the "covariance matching method" is different from the standard least squares fitting. It introduces a forgetting factor in the recursive process to balance the influence of historical and current residuals on the model weights; the "long-distance pipeline model" and the "end partition model" are different from the uncoupled model. They are based on transmission line theory and distributed heat capacity chain theory, respectively, and accurately characterize the dominant dynamics of hydraulics and thermals under different spatial scales and working conditions; the "normalized residual" is different from the direct residual value. It achieves comparable error measurement between different models by dividing the residual by a preset accuracy or range.

[0041] Step 5: Use online quadratic programming method according to , actuator upper and lower bounds and future load forecast curves to solve the nonlinear optimization problem involving temperature error, energy consumption and actuator constraints, and generate an initial control sequence; It should be noted that the "online quadratic programming method" in this application is different from traditional real-time optimization or rolling optimization. It not only targets the multivariable control problem of the water-heat coupling system of the pipe network, but also integrates the temperature tracking error, energy consumption index and actuator physical constraints into the same quadratic optimization framework, which can generate a smooth and implementable control sequence in a short time; For details, please refer to Figure 4 The comparison shows the indoor temperature error and control signal changes within the same time window.

[0042] During implementation, the weighted prediction model group generated in the previous step and the heat load and cooling load prediction curves within a certain time domain in the future (for example, the next ten-minute control window, divided into several steps) are first used as input, and the latest actuator status is retrieved, including the opening of each valve, pump frequency and speed, and heat source water supply temperature. After parallel simulation based on the prediction model group, the system obtains a multi-model joint prediction of equipment behavior and indoor terminal temperature at each prediction step. In order to meet the needs of multi-objective control, the system forms a temperature tracking error sequence with the difference between the predicted temperature and the set temperature, and forms an energy consumption penalty term with the sum of the control quantity fluctuation and the accumulated energy consumption, and weaves the two into a quadratic objective function according to the pre-set comfort weight and energy efficiency weight. In this process, users can set the relative weight coefficients of comfort and energy efficiency according to the building type and energy consumption cost. This application does not impose a rigid limit on their values.

[0043] Once the objective function and constraints are defined, the system must integrate the upper and lower bounds of the actuators: each valve has a physical limit between fully closed and fully open, each pump group has a safety limit between the minimum and maximum frequency, and the heat source has equipment and comfort boundaries for the minimum and maximum outlet water temperatures. To prevent control instructions from exceeding the actuator's capabilities, the algorithm adds these constraints to the optimization problem in the form of inequalities. It also considers the maximum variation in the control quantity at adjacent moments to avoid water hammer or thermal shock caused by sudden changes. For specific configuration, refer to the actuator's manual for the technical specifications regarding response time and maximum ramp rate. The specific values ​​can be determined based on the performance and service life requirements of the on-site equipment and are not limited.

[0044] This optimization problem is essentially a quadratic programming (QP) problem coupled with multi-stage dynamic prediction. To achieve real-time solution with limited computing resources, this application proposes a solution strategy that combines Krylov subspace projection with preconditioned conjugate gradients. First, the system constructs a set of linear equations based on the KKT conditions, which include the Hessian matrix, gradient vector, and constrained Jacobian matrix. Due to the large number of pipeline network nodes and the number of future prediction steps, directly solving this linear system is extremely expensive. To this end, the system generates a Krylov subspace only for the current gradient and Hessian during each optimization cycle: the initial residual vector and the Hessian-residual product are repeatedly multiplied to generate a basis set with an equal dimension much smaller than the number of original variables. This process eliminates the need for full matrix storage and only requires multiple matrix-vector multiplications, fully leveraging the sparse matrix properties of the pipeline network's physically coupled structure. The subspace dimension can be flexibly selected between 10 and 50 dimensions, depending on the performance of the edge computing node and the control cycle requirements. The specific dimension can be determined based on the on-site hardware and system scale, and is not limited.

[0045] After obtaining a Krylov subspace basis, the system projects the KKT system into this subspace, forming a small-scale equivalent linear system. A preconditioned conjugate gradient method is then applied for iterative solution. The preconditioner typically uses an approximate inverse matrix based on a diagonal or incomplete LU decomposition to further reduce the condition number and accelerate convergence. In each iteration, the step size is estimated using the Barzilai–Borwein formula. This adaptive step size strategy eliminates the need for multiple objective function evaluations required by traditional line searches, significantly improving iteration speed while ensuring the validity of the search direction.

[0046] During the iteration process, the system continuously monitors constraint consistency. When a conflict is detected in some constraints—for example, when the temperature tracking objective and the energy minimization objective cannot be simultaneously met within given actuator limits—the system automatically activates a relaxation factor mechanism to dynamically relax the inequality bounds. This relaxation is achieved by adding or subtracting a certain percentage of the slack from the bounds. The size of the buffer factor can be adjusted online based on the severity of the conflict and the tolerance for comfort. The specific amount can be determined based on safety regulations and comfort standards and is not limited. The principle of minimizing the relaxation range is to recover a feasible solution in real time while avoiding excessive relaxation that could lead to a significant decrease in comfort.

[0047] During conjugate gradient iteration, if the residual rate of decrease remains below a preset threshold after several consecutive iterations (for example, if the residual decreases by 10% or less per step, which can be set based on a pre-defined iteration efficiency target and is not a limit), the algorithm triggers a heuristic restart mechanism: the current residual is used as the new search direction, the conjugate direction is reset, and the Krylov subspace basis is updated to avoid the potential long slot dilemma and ensure continued effective optimization. The restart criteria and restart frequency can also be flexibly adjusted based on the dynamic characteristics of the pipeline network and the control frequency, and are not limited.

[0048] When the preconditioned conjugate gradient iteration reaches the termination condition—including the residual norm falling below the absolute threshold or the number of iteration steps reaching the maximum limit—the algorithm maps the optimal solution in the output dimensionality reduction subspace back to the original decision variable space to form a complete control sequence To prevent mechanical or thermal shock caused by sudden execution changes, the system applies a smoothing filter to each control variable before issuing the control sequence. The system subtracts the value from the previous cycle's execution value and multiplies it by a smoothing factor (typically 0.2-0.5, but can be determined based on comfortable jitter tolerance and is not limited). This yields the final execution instruction. These instructions are issued to specific actuators via industrial protocols—such as the water supply temperature setpoint for a boiler or heat pump, the target frequency for a pump inverter, or the opening percentage of an electronic proportional valve—and continuously obtains execution feedback to ensure closed-loop reliability.

[0049] As shown in Table 1, the sensitivity analysis table of the parameters of the online quadratic programming method is used to evaluate the impact of each parameter on the convergence performance and comfort: Table 1 Sensitivity analysis table Parameter name Optional range Impact on convergence performance Impact on comfort Traffic threshold 0.1-1.0m³ / h The smaller the threshold, the slower the convergence may be. Too small will lead to too frequent control and affect comfort Forgetting factor λ 0.9-0.999 The larger λ is, the greater the influence of historical data is, and the convergence may be more stable. Too large a λ value may cause system response to be slow and affect temperature tracking. BB step size α 0.5-2.0 α affects the search step size, and a suitable value accelerates convergence Inappropriate α may lead to increased fluctuation of the control signal Krylov subspace dimensions 10-100 The higher the dimension, the worse the dimensionality reduction effect and the more computational effort Too low a dimension may affect the convergence accuracy and increase the temperature error It should also be noted that the "weighted prediction model group" in this application is different from the prediction of a single model. It is composed of the weighted outputs of multiple sub-models and can provide continuous and robust state prediction under working condition switching or coexistence scenarios; "Krylov subspace projection" is different from simple dimensionality reduction or feature decomposition. It uses the quadratic Hessian information of the objective function to construct a low-dimensional search subspace, which greatly reduces the computational complexity of solving the KKT system; "preconditioned conjugate gradient method" is different from the unpreconditioned conjugate gradient. It accelerates the convergence of the Krylov method in the pipeline coupling optimization problem by targeted selection of preconditioned matrices; "BB step size" is different from the fixed step size or line search step size. It adaptively estimates the approximate optimal step size in each iteration through the Barzilai–Borwein formula, significantly improving the iteration efficiency; "heuristic restart" is different from the fixed period restart. It dynamically determines when to reset the conjugate direction based on the residual decrease rate to avoid search stagnation or entering a pathological direction.

[0050] Step 6: Generate actual pump speed and valve opening control instructions according to the initial control sequence and decoupling compensation terms in accordance with the pump-first-valve-later sequencing strategy and in combination with dead zone and rate limit; The core of step 6 is to convert these intentions into specific instruction sets that can be directly issued to the pump frequency inverter and electronic proportional valve. This not only meets the physical limits and response characteristics of the actuators, but also ensures that the water-heat coupling compensation logic is implemented in the control action. This process mainly includes instruction extraction and classification, timestamp sorting and trigger criteria, deadband determination and action delay, rate planning and feedforward correction, and labeled instruction output.

[0051] First, the system parses two types of instructions from the initial control sequence: one is the pump speed instruction , representing the target speed or power of the inverter of each pump group; the other type is the valve opening instruction , representing the target opening percentage of each control valve. During the parsing process, the controller will specifically superimpose the feedforward decoupling compensation term on the valve instruction set to ensure that after the pump speed changes, the valve can compensate for the pressure or flow coupling error introduced accordingly without adjusting the pump frequency again. It should be noted that the "decoupling compensation term" in this application is different from the traditional feedback correction. It reduces the burden on the closed-loop controller by pre-offsetting the coupling interference at the instruction level; the specific value of the compensation amount can be dynamically determined in combination with the matrix product result generated in step 3, and is not limited.

[0052] After parsing is completed, all instructions will be timestamped, that is, each pump speed or valve opening instruction carries a time tag. , indicating that the instruction should be After the timestamp is marked (usually an integer multiple of the control cycle), it is sent for execution. The controller can strictly sort the pump speed command and the valve command: first, the pump speed command is sent according to the timestamp. Arrange in ascending order and send to the pump frequency converter interface in sequence; after each pump speed instruction is sent, the system will call the pressure gradient detector module in real time to monitor the pressure gradient of adjacent measuring points in the main pipeline and its rate of change When the pressure change rate is detected to enter the steady-state dead zone range (For example =0.001 MPa / s), the system triggers the valve opening instruction of the corresponding timestamp; otherwise, the valve instruction will be delayed until the steady-state condition is met or the maximum delay time is reached. (The typical value is 5 seconds, which can be determined based on the pipe network length and response requirements and is not limited.) This strategy of pumping first and then activating the valve, combined with the dead zone criterion, can prevent pressure waves caused by pump speed fluctuations from being reflected back to the main trunk and causing water hammer when the valve is not closed or opened in time.

[0053] For the specific execution of valve opening instruction, the system uses segmented S-type rate planner to generate a smooth opening curve. To target opening The change is divided into three stages: acceleration stage (the opening increment is accelerated from 0 to the maximum rate ), uniform speed stage (maintain ), and the deceleration phase (from decelerate uniformly to 0). Maximum allowable rate of change Hydraulic diameter of pipe section Determined by the current flow rate Q, usually based on the empirical formula in is the empirical coefficient (e.g. =0.1~0.5, which can be adjusted according to the valve actuator characteristics and the system water hammer tolerance, without limitation), to ensure that the opening changes faster when the pipe diameter is large and the flow rate is high, while the change rate needs to be reduced to avoid shock when the pipe diameter is small and the flow rate is low. The sliding window constraint calculates the actual rate in each small time step (such as 100ms) and ensures When the system detects a sudden change in traffic ( ,For example =2 ), indicating that the coupling compensation may be insufficient, the controller will enable the feedforward compensation term to correct the rate threshold: Multiply by a magnification factor >1 (range: 1.1~1.5) to respond to sudden disturbances more quickly. The specific factor can be determined according to the system safety constraints and comfort requirements and is not limited.

[0054] After the rate curve is generated, valve opening commands are queued according to their timestamps. Each opening curve is discretized into several steps at the underlying actuator interface and sent in batches, balancing the controller's communication bandwidth with the actuator's minimum motion resolution. For actuators that support batch writes to registers, the system can send the node values ​​for the entire S-curve all at once and initiate hardware self-execution. For devices that don't support batch writes, the next opening value is sent at each small time step via a timed interrupt to ensure complete tracking of the opening curve.

[0055] After all pump speed and valve commands have been sent and feedback has been received, the controller will recollect the actual actuator feedback values ​​in the next control cycle, including the pump speed achievement rate and actual valve opening, and compare them with the command values ​​to calculate the execution error. If the execution error exceeds the set threshold (pump speed error > 2%FS or valve opening error > 1%FS, which can be set based on actuator accuracy and is not limited), the system will automatically increase the corresponding compensation amount in step 3 or step 5 to ensure closed-loop robustness.

[0056] From the perspective of system behavior, step 6 ensures that the pump speed change is stabilized first through the sequential strategy of pump first and then valve, and then the valve action compensates for the coupled pressure error; the steady-state dead zone judgment is used to avoid the valve from malfunctioning before the pressure wave subsides; smooth and controllable valve dynamics are achieved through S-shaped rate planning and sliding window constraints; sudden flow disturbances are dealt with through the feedforward rate correction mechanism; sequential execution and real-time synchronization are guaranteed through timestamp marking and delay mechanism; the final output control instruction set can not only meet the physical characteristics and safety constraints of the equipment, but also accurately implement the decoupling compensation and optimization strategies generated in steps 3 and 5, providing a solid execution guarantee for the temperature control of the next cycle.

[0057] It should be noted that the "pump-first, valve-later sequencing strategy" in this application is different from the traditional simultaneous drive mode. It can eliminate the transient coupling of mutual interference between pump speed and valve through the time sequence relationship, thereby ensuring the predictability of pressure wave propagation and flow distribution; the "steady-state dead zone" is different from the dead zone in general control. It is defined as when the pressure change rate of the main pipe of the pipeline network is within a certain small range (such as ±0.001MPa / s), the system is considered to have entered a steady state and can safely execute the downstream valve action without causing secondary oscillations; the "S-type rate planner" is different from the linear ramp rate control. It generates a smooth opening change through a speed curve of three stages: acceleration-uniform speed-deceleration to reduce water hammer and mechanical shock; the "pressure gradient detector" is different from single-point pressure monitoring. It compares the pressure difference between adjacent measuring points in real time and calculates the gradient change rate to determine whether the system is in an allowed state.

[0058] Step 7: After the control is executed, the temperature residual and flow oscillation data of the previous moment are used as input to dynamically adjust the multi-model weights and prediction time domain length, and output the updated control configuration parameters. Specifically, the difference between the indoor temperature prediction and the measured temperature at the previous moment constitutes the temperature residual sequence , the flow sensor synchronously records the instantaneous flow of each branch or trunk The system first takes the residual sequence as input and calls the Hilbert-Huang transform module to Perform empirical mode decomposition (EMD): The algorithm decomposes the original residual signal into several IMF components through extreme point envelope construction and iterative screening (k),…, (k) and the residual trend term. Then, the time domain energy of each IMF is calculated And normalize to get the energy spectrum weight The energy ratio of all low-frequency IMFs (usually those with larger numbers and time domain periods exceeding twice the control period) in the system (For example, 0.6~0.8, which can be set according to the building insulation characteristics and thermal inertia time constant, and is not limited). If the proportion of low-frequency energy exceeds the threshold, it is determined that the temperature residual is mainly dominated by thermal inertia lag. At this time, the time domain length of the multi-model prediction needs to be extended to capture the slow-changing trend, so as to avoid the short-window prediction from missing significant thermal lag components. The prediction time domain length extension mechanism will extend the current prediction window length. Heat capacity time constant Incremental correction is made by 1.2 times: in The ratio of the node heat capacity parameter to the heat transfer coefficient, obtained from the online identification in step 2, represents the time it takes for the system to reach steady-state at the 1 / e level. It should be noted that the "heat capacity time constant" in this application is distinct from a general RC circuit constant; it specifically refers to the lumped response time of the thermal inertia of the entire thermal network system. The expansion range (1.2x) and the ultimate upper limit prediction window can be flexibly adjusted based on the system's computing power and control accuracy requirements (e.g., a maximum of 20 minutes) and are not limited.

[0059] At the same time, the system takes flow oscillation data as input and targets a typical measuring point or multi-point combined signal. Calculate the autocorrelation function And normalized to The main lobe width of the autocorrelation function Defined as the signal from Decay to critical level (such as 1 / e or 0.5) the required delay time. (The typical threshold is between 1s and 5s, which can be set according to the response characteristics and oscillation period of a single pump or valve, and is not limited). This indicates that the flow oscillation has high-frequency components and a short duration. It is determined that hydraulic interference is dominant in this area. At this time, the multi-model weights need to be redistributed to reduce the weight of the hydraulic model in this area and improve the confidence of the thermal model. The system uses the oscillation phase recognition algorithm to identify the Hilbert instantaneous phase. Extract the main frequency component of the traffic signal and divide it into phase-locked subnets, and regard the branches or partitions with synchronized phases and sudden amplitude changes as the same associated subnet. Weights of the (hydraulic dominated model) Perform a linear downscaling: And will reduce the ∈(0.1,0.3) is evenly distributed to all involved areas of the subnet; at the same time, the thermal model of the corresponding area Weight = +β. It should be noted that "oscillation phase identification" in this application is different from simple time-domain oscillation detection. It combines the Hilbert transform with the phase continuity criterion to distinguish between synchronous and asynchronous behaviors of multi-region flow oscillations. The weight adjustment ratio β can be adjusted within the range of 0.05 to 0.5 based on the oscillation intensity and control target, without limitation.

[0060] For areas where the oscillation criterion is not triggered, the system increases the confidence of the thermal model in the overall weight distribution: that is, all remaining model weights are renormalized proportionally and the thermal model weights are Multiply by the amplification factor (1+γ), γ∈(0.05,0.2) represents the trust increase in the thermodynamic model. In the same control cycle, the weights of all regions need to be renormalized after adjustment to ensure This dynamic update mechanism enables the model pool to maintain immediate preferences for different physical processes in scenarios where hydraulic and thermal dominance switch or coexist, thereby improving the robustness of prediction and control.

[0061] Once the update is complete, the system injects the new model pool weight vectors and prediction time-domain parameter sets into the multi-model prediction and online quadratic programming modules, providing the latest structured configuration for the next steps of predictive modeling and control optimization. This process is automatically completed within each control cycle (which can be executed periodically between 30 and 120 seconds). This ensures that the system can quickly adjust the prediction window and model confidence level under various disturbance conditions, such as sudden changes in indoor heat load, main valve adjustments, or switching of heat sources, achieving high-precision, dynamic, adaptive control of the thermal pipeline network.

[0062] It should be noted that the "Hilbert-Huang transform" in this application is different from the traditional Fourier or wavelet time-frequency analysis. It can adaptively decompose several intrinsic mode functions (IMFs) in the context of nonlinear and non-stationary signals and extract instantaneous frequency and amplitude information, which can reflect the multi-scale oscillation characteristics implicit in the indoor temperature residual signal; "intrinsic mode function energy spectrum" is different from the spectrum density. It refers to the proportion of the time domain energy of each IMF (i.e., the square integral of the IMF) to the total energy, which can intuitively quantify the low-frequency lag and high-frequency fluctuation components; "prediction time domain length" is different from the control period length. It refers to the model in online quadratic programming or multi-model The size of the future time window used in the prediction can affect the prediction accuracy and computational burden; the "autocorrelation function decay rate" is different from the correlation coefficient. It reflects the signal's memory length and coupling duration through the delay time required for the autocorrelation sequence peak to decay to 1 / e or 1 / 2; the "main lobe width" is different from the concept of frequency domain main lobe. It corresponds to the half-width of the main peak in the time domain of the autocorrelation function and can indicate the duration of flow oscillation; "oscillation phase identification" is different from ordinary phase identification. It uses Hilbert instantaneous phase analysis or phase envelope method to group the flow oscillations in each partition of the pipeline network, so as to identify the coupling area when the multi-model weights are redistributed.

[0063] See also Figure 5As a possible implementation of a multi-parameter coupled intelligent control system for the vanadium smelting process, the specific working principle is as follows: In actual operation, each module closely coordinates with the message bus through a clearly defined data interface to ensure that the output of each link can directly serve as the input of the downstream link. Specifically, the data acquisition module first completes the raw measurement of temperature, flow, pressure, and valve position through the edge node. After local filtering and clock synchronization, it pushes the standardized time series data set D1 to the central message bus. The time series preprocessing module uses this message bus as the subscription port, automatically pulls D1, performs wavelet packet multi-scale denoising and Bayesian interpolation padding on it, and then forwards the smoothed and denoised sequence S1 to the parameter modeling module via the same bus.

[0064] Upon receiving S1, the parametric modeling module immediately calls the network topology database service to load the static topological relationships and physical properties of nodes and pipe segments. It then constructs the full-network model M0 using a hybrid PDE / lumped parameter strategy. Once M0 is constructed, its state equations and parameter vector θ0 are encapsulated as a "model initialization" message and sent to the input queue of the online identification module.

[0065] The online identification module monitors the "model initialization" and "smoothing sequence" inputs. It first performs state prediction and residual calculation based on the interpolated S1 and M0. It then automatically switches between the EKF and UKF subprocesses based on the amplitude threshold of the residual covariance matrix. Through recursive calculation, it outputs the latest water impedance and heat capacity parameters θ1, as well as the updated model M1. This updated model and parameters are then published to the coupling and decoupling modules via a "model update" message.

[0066] After receiving M1, the coupling-decoupling module extracts the network hydraulic and thermal coupling matrices from it, performs spectral decomposition on them, screens the dominant coupling subspace, and generates a feedforward decoupling operator C. This operator, as a decoupling compensation term before quadratic programming, is pushed to the predictive control module along with M1 via a "decoupling operator" message.

[0067] The predictive control module monitors "decoupling operator" messages and simultaneously pulls the future load curve L(t+i) from an external load forecasting service. It takes M1, C, and L as inputs and, through a parallel multi-model adaptive ensemble and a linear quadratic programming solver, outputs the initial control vector U0 with decoupling compensation. U0 is then sent to the execution monitoring module via an "initial control sequence" message.

[0068] After receiving U0, the execution monitoring module first splits it into pump speed commands and valve opening commands, and then sends them to the on-site pump control unit and valve driver unit in sorted order by timestamp. The pump control unit responds to the valve driver unit's opening command only after executing the pump speed command and verifying, via the pressure gradient detector, that the main pipeline pressure has entered the dead zone. The execution monitoring module also continuously collects feedback from the equipment—including real-time pump speed values ​​and actual valve travel—and sends this feedback as a residual sequence R to the compensation iteration module.

[0069] The compensation iteration module subscribes to the “execution feedback” message, uses the residual sequence R and the historical residual library to run the iterative learning control algorithm, generates the compensation increment ΔU_ILC, and superimposes ΔU_ILC with U0 of the next cycle to form a corrected control vector U′, which is then input into the closed-loop coordination module again through the “corrected control” message to achieve pre-compensation of future control sequences.

[0070] The closed-loop coordination module monitors two messages: "corrected control" and "subnetwork interface error." The latter comes from the interface flow and pressure monitors of each local subnetwork. Taking ΔU_ILC, R, and the interface error E as input, the closed-loop coordination module coordinates the model parameters and MPC weights of the multiple subnetworks using an augmented Lagrangian iterative algorithm. Upon convergence, it outputs a new set of model parameters θ2 and MPC weight vector ω. These, along with the control configuration for the next cycle, are sent back to the predictive control module via a "new cycle configuration" message, forming an end-to-end closed-loop adaptive control process.

[0071] Each of these message push and subscriptions runs on a reliable industrial message bus (such as OPC UA or Kafka Industrial Edition) to ensure real-time and reliable data transmission. All modules are horizontally scalable to accommodate smelting systems of varying sizes, from single furnace islands to complex multi-furnace parallel operations, enabling end-to-end intelligent control of multi-parameter coupling in the vanadium ore smelting process.

[0072] During implementation, temperature, flow, pressure, and valve position sensors distributed across each heat exchange station and terminal room continuously upload raw signals with unified clock timestamps to the data acquisition module. The acquisition module aligns these signals in memory according to their timestamps and then pushes them as standardized time series data sets to the time series preprocessing module. Upon receiving the aligned raw sequences, the time series preprocessing module immediately activates a wavelet packet multiscale filter to separate high-frequency mechanical vibration noise from electrical interference. It also uses Bayesian interpolation to fill in missing values ​​caused by maintenance or communication delays, ultimately outputting a smoothed, noise-reduced, and complete time series signal for downstream use.

[0073] The parameter modeling module continuously receives this set of trusted time series from the preprocessing module and associates each measurement point to a node or edge segment of the model based on the physical topology of the pipe network. Using a hybrid PDE / lumped parameter strategy, it constructs an initial set of network-wide state equations, including second-order thermal inertia and hydraulic coupling terms. This model, along with the latest flow and temperature data, is then submitted to the online identification module. The online identification module uses this model as the basis for prediction, taking the smoothed time series data as observation input. By alternating between extended Kalman filtering and unscented Kalman filtering, it continuously corrects the water impedance and node heat capacity parameters. It then synchronously returns the updated coupling model to the parameter modeling module and subsequent coupling and decoupling modules to ensure that the models used by each module reflect the latest physical characteristics of the site.

[0074] After obtaining the latest model, the coupling-decoupling module first extracts the network hydraulic coupling coefficient matrix and the thermal coupling coefficient matrix from it, and performs spectral decomposition on them respectively - the inverse operation of the node admittance matrix combined with Jacobian iteration generates the hydraulic subspace, and the singular value decomposition reconstructs the thermal subspace; then, these two subspaces are jointly operated to form a feedforward decoupling operator, and this operator and the updated model are pushed to the predictive control module together to compensate for the coupling interference in the rolling time domain optimization.

[0075] The predictive control module takes the decoupling operator, the updated model, and the future load forecast sequence as input, constructing a multi-model parallel prediction framework. After generating weighted prediction results, it solves the initial control vector using a rolling-horizon quadratic programming approach and then delivers this initial sequence to the execution monitoring module. After receiving these control commands, the execution monitoring module, following a pump-first, valve-last sequencing strategy, timestamps the pump speed and valve opening commands, sorts them, and issues them. Simultaneously, it monitors the response state using a pressure gradient detector and flowmeter. Valve action is triggered only when steady-state criteria are met. The room temperature residual and flow deviation signals after execution are then fed back to the compensation iteration module in real time.

[0076] The compensation iteration module uses the latest temperature residual, flow deviation and historical residual library as the training set, runs an iterative learning algorithm (such as recursive least squares or online reinforcement learning) to generate deviation compensation increments, and pushes the corrected control vector back to the execution monitoring module to complete the closed-loop compensation; at the same time, it reports the compensation increment together with the residual after execution to the closed-loop coordination module.

[0077] The closed-loop coordination module inputs the compensation increment, the interface errors of each subnetwork, and the current model weights into the augmented Lagrangian iterative process. By dynamically adjusting the Lagrangian multipliers and relaxation factors, it not only updates the model parameters and MPC weights, but also sends the new round of control configuration parameters back to the predictive control module and the execution monitoring module, starting the closed-loop optimization cycle for the next period.

[0078] Although specific embodiments of the present invention have been described above, those skilled in the art will appreciate that these specific embodiments are merely illustrative, and that those skilled in the art may omit, substitute, and modify the details of the methods and systems described above without departing from the principles and spirit of the present invention. For example, combining the above method steps to perform substantially the same functions in substantially the same manner to achieve substantially the same results falls within the scope of the present invention. Accordingly, the scope of the present invention is limited solely by the appended claims.

Claims

1. A method for real-time control of indoor temperature of a thermal pipeline; characterized by: include: Step 1: Based on the node temperature, pressure and flow data collected in real time from the pipe network, the PDE / lumped parameter hybrid method is used to construct an initial full-network model including second-order thermal inertia and hydraulic coupling terms. , and output the initial state equations; Step 2: The extended Kalman filter and the unscented Kalman filter are used alternately to identify the water impedance and heat capacity parameters online and obtain the updated model. ; Step 3: Based on the updated model The parameter vector and state equation group of the network are used to calculate the hydraulic coupling coefficient matrix and the thermal coupling coefficient matrix, and generate the corresponding feedforward decoupling compensation terms; Step 4: Based on the feedforward decoupling compensation term and The parameter vector is constructed to build a multi-model prediction framework including the long-distance pipeline model and the terminal partition model, and the weighted prediction model group is output through the model pool weight scheduling method. ; Step 5: Use online quadratic programming method according to , the upper and lower bounds of the actuator and the future load forecast curve are used to solve the nonlinear optimization problem including temperature error, energy consumption and actuator constraints to generate the initial control sequence; Step 6: Generate actual pump speed and valve opening control instructions according to the initial control sequence and decoupling compensation terms in accordance with the pump-first-valve-later sequencing strategy and in combination with dead zone and rate limit; Step 7: After the control is executed, the temperature residual and flow oscillation data at the previous moment are used as input to dynamically adjust the multi-model weights and the prediction time domain length, and output the updated control configuration parameters.

2. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The working principle of the PDE / lumped parameter hybrid method is as follows: the pressure gradient distribution of the pipeline network nodes is collected in real time through the pressure sensor network, and based on the dynamic hydraulic domain partitioning mechanism of Morris sensitivity analysis, the main pipeline that meets the pressure gradient threshold is classified into the PDE domain, and the transient flow conservation equation is discretely solved using the characteristic line method; at the same time, the branch pulsation amplitude is identified through the flow meter data. If it is lower than the turbulence critical value, the branch subnet is classified into the Lumped domain, and the lumped parameter model is constructed using the node method; the PDE domain and the Lumped domain trigger bidirectional data synchronization based on the real-time exchange of pressure and flow interface conditions at the coupling boundary; for thermal modeling, the axial temperature distribution is obtained through distributed temperature sensors, and the dominant spatial mode is extracted using the KL decomposition method and reconstructed into a low-dimensional state-space equation; at the same time, the pipe wall heat capacity parameters are corrected online using the actual measured values ​​of the heat flow meter, and finally the hydraulic PDE solution and the thermal dimensionality reduction model are integrated to output a coupled state equation group.

3. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The steps of the alternating identification method of extended Kalman filtering and unscented Kalman filtering include: Through the initial full network model Generate state prediction value, compare it with the interpolated and smoothed sensor data and output the residual sequence; The extended Kalman filter is triggered when the Frobenius norm of the residual covariance matrix is ​​lower than the preset threshold. Perform first-order Taylor linearization, construct the Kalman gain using the state Jacobian matrix, and recursively derive the water impedance and node heat capacity parameters in the prediction-update loop to output the extended Kalman filter correction value; If the Frobenius norm of the residual covariance matrix exceeds the preset threshold, the weighted Sigma point is generated by the unscented Kalman filter and the prior covariance is calculated on the state prediction and observation mapping respectively using the Unscented transform, and then the parameter vector is updated in combination with the observation residual; After each filtering, the noise matrices of the extended Kalman filter and the unscented Kalman filter are dynamically adjusted based on the covariance reconstruction algorithm of the observation noise and process noise; The extended Kalman filter and unscented Kalman filter outputs of adjacent periods are fused through the posterior smoothing algorithm to filter out mutation errors and maintain parameter continuity, and the updated model is output after model consistency detection.

4. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The calculation steps of step 3 include: s100, based on the water impedance parameters and node heat capacity parameters in the updated model output in step 2, construct the hydraulic coupling coefficient matrix through the node admittance matrix analysis method ,in For nodes The pressure difference, For branch The node admittance matrix analysis method solves the sensitivity of the pressure gradient of the pipeline network node to the flow change through the Jacobi iteration method, and outputs a diagonal-dominated sparse matrix ; s200, according to the state equations of the heat capacity chain state space model, the state matrix Perform singular value decomposition and retain the first r significant singular values ​​to reconstruct the reduced-order matrix , extract the partial derivative of heat flux to stratification temperature and generate the thermal coupling coefficient matrix ,in is the temperature state of the kth layer, is the heat flux; s300, calculate hydraulic coupling matrix The off-diagonal norm of If it is greater than the preset threshold , calculate the pressure disturbance suppression term, and then generate the valve opening correction through the feedforward compensator. The formula is: （1） In formula (1), is the gain matrix, which is used to map the pressure fluctuation into the valve opening correction; Real-time pressure gradient residual, in Pa / m, used to reflect the uncompensated oscillation energy; s400, synchronous calculation of thermal-mechanical coupling matrix condition number ,exist Greater than the preset threshold When , the temperature lag compensation term is generated based on the pseudo-inverse matrix of the heat flux temperature transfer function through the thermal inertia decoupling algorithm. The formula is: （2） In formula (2), is the pseudo-inverse matrix of the heat flux-temperature transfer function, which is used to offset the multi-order lag effect of the building heat capacity; is the target heat flux deviation; s500, based on the real-time pipeline network pressure oscillation spectrum energy , and thermal inertia residual Ratio , the compensation terms are weighted and fused to generate the final feedforward decoupling instruction. The fusion formula is: .

5. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The multi-model prediction framework has a pipeline flow mutation rate greater than When the temperature gradient variance at the terminal is greater than When it is determined to be a thermally dominant operating condition, the terminal partition model is activated; the long-distance pipeline model is constructed based on the fluid mechanics wave equation, and the pressure wave propagation process is simulated by solving the one-dimensional unsteady Navier-Stokes equation; the terminal partition model adopts a variable-order heat capacity chain state space architecture, integrating the thermal coupling coefficient matrix of step 3 and compensation items Analyze the temperature response hysteresis effect caused by the thermal inertia of building envelope structures.

6. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The multi-model prediction framework corrects the initial parameters of each model based on the feedforward decoupling compensation term output in step 3, and then calculates the model residual sensitivity by the covariance matching method, where: Long-distance pipeline model residuals ; Terminal partition model residuals ; in and denote the norm of the deviation between measured and predicted pressure respectively; and Represent the measured and predicted temperature deviation norms respectively; then the recursive least squares algorithm with forgetting factor is used to update the model weights and output the weighted prediction model group; the update formula is: （3） In formula (3), Weight the model pool, and Respectively Hedi Normalized residual for each scene.

7. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The steps of the online quadratic programming method for solving a nonlinear optimization problem including temperature error, energy consumption, and actuator constraints include: By weighted prediction model group Generate the state prediction equation and construct a quadratic objective function with the goal of minimizing the control variable fluctuation and temperature tracking error; An inequality constraint matrix is ​​constructed based on the upper and lower bounds of the actuators, and the future load forecast curve is introduced as the right-hand side term of the equality constraint. The Krylov subspace projection method is used to compress the dimension of decision variables, and the preconditioned conjugate gradient method is triggered to solve the KKT system after dimensionality reduction. If a constraint conflict is detected, the inequality boundary tolerance is dynamically adjusted through the relaxation factor, and the BB step size is used to accelerate the iterative convergence; The residual norm is calculated after each iteration. If the residual decrease rate is lower than the preset threshold, the heuristic restart mechanism is triggered to reinitialize the conjugate direction vector.

8. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The working principle of step 6 is as follows: Extract the pump speed instruction set and the valve opening instruction set through the initial control sequence, so that the decoupling compensation term acts on the valve opening set to offset the coupling interference; A timestamp tagging mechanism is used to sequence instructions, giving priority to executing pump speed instructions, and a pressure gradient detector is used to monitor the main pipeline pressure propagation status; If the pressure change rate enters the steady-state dead zone, the valve command enable signal is triggered, otherwise the valve action is delayed; For valve opening instructions, a segmented S-type rate planner is used to calculate the maximum allowable change slope based on the hydraulic diameter of the pipe section, and the instantaneous adjustment rate is constrained by a sliding window. When a sudden change in traffic is detected, the feedforward compensation item is enabled to correct the rate limit threshold; Output actual pump speed and valve opening control instruction set with timing tags.

9. The method for real-time control of indoor temperature of a thermal pipeline according to claim 1, characterized in that: The dynamic update method of step 7 is as follows: the intrinsic mode function energy spectrum of the temperature residual sequence at the previous moment is extracted through the Hilbert-Huang transform; if the low-frequency IMF energy ratio exceeds a preset threshold, it is determined to be a thermal inertia-dominated error, and a prediction time domain length extension mechanism is used to correct it by an increment of 1.2 times the heat capacity time constant; at the same time, the autocorrelation function decay rate is calculated based on the flow oscillation data. If the main lobe width is lower than a preset critical value, the multi-model weight redistribution mechanism is triggered, and the associated subnet is identified through the oscillation phase, and the hydraulic model weight of the corresponding area is reduced; For non-oscillating areas, the confidence of the regional thermal model is improved; finally, the updated model pool weight vector and the predicted time domain parameter set are output.

10. A real-time control system for indoor temperature of a thermal pipeline, characterized by: A method for real-time control of indoor temperature of a thermal pipeline as claimed in any one of claims 1 to 9, comprising: The data acquisition module is used to collect the original signals of temperature, flow, pressure and valve position of each heat exchange station and the terminal room and align the timestamps according to the unified clock to obtain a standardized time series data set. ; Timing preprocessing module is used to Perform wavelet packet multi-scale filtering and Bayesian interpolation processing to separate noise and fill missing values ​​to obtain a smoothed and denoised time series. ; Parametric modeling module, used to Mapping to the pipe network topology and building a full network model with second-order thermal inertia and hydraulic coupling based on a PDE / lumped parameter hybrid strategy ; Online identification module, used to interpolate and smooth and As input, the water impedance and node heat capacity parameters are corrected alternately by extended Kalman filtering and unscented Kalman filtering to obtain the updated model ; Coupling and decoupling modules for Perform spectral decomposition on the hydraulic and thermal coupling matrix of the network, extract the subnetwork coupling subspace and generate the feedforward decoupling operator ; Predictive control module for , operator The parallel multi-model is constructed in the rolling time domain with the future load forecast sequence, and the initial control vector is solved by linear quadratic programming. ; Execution monitoring module, used to issue orders according to the sequence strategy of pump first and valve later , and collect room temperature residual and flow deviation series in real time ; Compensation iteration module, used to And the historical residual library is used as input to run the iterative learning algorithm to generate deviation compensation increments And output the corrected control vector ; Closed-loop coordination module, used to 、 The subnetwork interface error is input into the augmented Lagrangian iterative process together to update the model parameters and MPC weights and output the configuration for the next cycle.

Images