Method and system for energy efficiency distribution of multi-heat source heating and hot water

Abstract

This invention belongs to the field of water source heat pump technology, specifically relating to a method and system for energy efficiency allocation in multi-heat source heating and hot water systems. The method includes: constructing a multi-dimensional objective function; generating the load prediction mean and covariance matrix using Gaussian process regression; constructing a simplex time-varying solution space based on the extreme values ​​of the heat source health state, and constructing a comprehensive expected objective function in conjunction with load uncertainty; constructing a Riemann metric tensor through energy balance mapping, extracting the Riemann natural gradient descent to search for the optimal expected allocation point; and searching for the feasible allocation point with the minimum Riemann distance in the hard-constrained subspace and issuing the execution. This invention, by considering both equipment degradation and load uncertainty, eliminates space distortion and solves the problem of easy oscillation and boundary overflow during optimization under time-varying constraints, achieving stable and smooth optimal energy efficiency allocation under safety constraints.

Description

Technical Field

[0001] This invention relates to the field of water source heat pump technology. More specifically, this invention relates to a method and system for multi-heat source heating and hot water energy efficiency distribution. Background Technology

[0002] Multi-source heat pump systems typically integrate various types of heating equipment, such as solar energy, air-source heat pumps, gas boilers, and ground-source heat pumps, to meet the growing heating and domestic hot water comfort demands of various buildings. In actual operation, different heat source devices exhibit significant differences in energy efficiency characteristics, operating costs, response speed, and environmental emissions. Furthermore, the terminal heating and hot water loads are easily affected by changes in weather conditions, seasonal alternations, and users' random energy consumption behaviors, exhibiting obvious time-varying and uncertainties. Therefore, how to allocate the total heating load to various heat source devices in real-time and rationally under complex operating conditions to achieve comprehensive optimization of the system in terms of energy efficiency, economy, and environmental benefits has become a crucial problem that urgently needs to be solved in the fields of thermal system control and building energy management.

[0003] Existing multi-heat source allocation strategies mostly rely on static priority rules or simple empirical feedback, making it difficult to fully utilize front-end load forecasting results and effectively address operational fluctuations caused by forecast uncertainty. Furthermore, traditional methods often fail to adequately consider the health degradation of heat source equipment during long-term operation and the corresponding changes in the extreme values ​​of actual available output. This can easily lead to load allocation commands exceeding the safe operating boundaries of the equipment, resulting in risks such as decreased system energy efficiency, increased equipment wear, or heating interruptions. In addition, during multi-dimensional comprehensive objective optimization, existing technologies lack an effective characterization of the geometric structure of the multivariate allocation solution space and the changing patterns of its constraint boundaries. Conventional optimization algorithms are prone to getting trapped in local extrema within time-varying constraint spaces or experiencing computational oscillations near the feasible region boundaries, making it difficult to balance real-time performance, stability, and feasibility. Therefore, the industry urgently needs a multi-heat source energy efficiency allocation method that can integrate equipment health status, load forecast uncertainty, and multi-dimensional operational performance indicators, and can stably and in real-time solve problems within time-varying constraint spaces, to improve the comprehensive operational value of multi-heat source combined heating systems. Summary of the Invention

[0004] To address the technical problem that existing technologies fail to consider both equipment degradation and load uncertainty, leading to easy oscillations and out-of-bounds errors during optimization under time-varying constraints, this invention provides solutions in the following aspects.

[0005] In a first aspect, the present invention provides a method for multi-heat source heating and hot water energy efficiency allocation, comprising: S1, acquiring real-time operating parameters and environmental conditions of the multi-heat source system, constructing a multi-dimensional objective function representing energy efficiency, economy, and environment; processing historical operating data and weather forecast data using a Gaussian process regression model to generate a mean vector and a load forecast covariance matrix for terminal load prediction; S2, assessing the available output extreme values ​​based on the health status of the heat sources, and constructing a single simplex or a set of simplexes as a time-varying solution space, using the extreme allocation point where each heat source satisfies the system's total load under the constraint of the available output extreme values ​​as the vertex; and generating the corresponding time node uncertainty based on the load forecast covariance matrix. S3. Using the energy balance mapping matrix, the inverse matrix of the load prediction covariance matrix is ​​projected onto the time-varying solution space to construct a Riemann metric tensor. The Euclidean gradient obtained by adjusting the partial derivative of the comprehensive expected objective function with respect to the Riemann metric tensor is used to obtain the Riemann natural gradient. The comprehensive expected objective function is then subjected to a descent search to converge to the optimal expected allocation point. S4. Within the subspace that satisfies the hard constraints of the equipment, a feasible allocation point with the smallest Riemann distance to the optimal expected allocation point is searched. The feasible allocation point is then inversely mapped to the power setpoints of each heat source and sent for execution.

[0006] This invention determines the available output extremum based on the health status of the heat source to construct a simplex solution space, taking into account equipment degradation. It utilizes the uncertainty scalar of the load prediction covariance matrix to adjust the objective function and constructs the Riemann natural gradient for search, taking into account load uncertainty and eliminating space distortion. It solves the problem that existing technologies are prone to computational oscillations and allocation out-of-bounds when seeking optimization under time-varying constraints, and achieves stable and smooth optimal energy efficiency allocation under safety constraints.

[0007] Preferably, the construction of the multidimensional objective function representing energy efficiency, economy, and environment includes: using the ratio of the total heat supply of the system to the total calorific value of the fuel consumed as the energy efficiency sub-objective function; using the sum of the energy consumption of each heat source multiplied by the real-time time-of-use electricity price or natural gas price as the economic sub-objective function; using the sum of the energy consumption of each heat source multiplied by the corresponding energy carbon emission factor as the environmental sub-objective function; and combining the reciprocal of the energy efficiency sub-objective function, the economic sub-objective function, and the environmental sub-objective function after performing dimensionless normalization processing to form the multidimensional objective function.

[0008] This invention achieves the unification of multidimensional optimization objectives under different physical dimension attributes by performing dimensionless normalization processing on the reciprocal of energy efficiency, economic cost and environmental emissions, thus eliminating the interference of single objective differences on the comprehensive optimization direction.

[0009] Preferably, generating the mean vector and load prediction covariance matrix of the terminal load forecast includes: inputting real-time weather forecast data into a Gaussian process regression model that optimizes hyperparameters using the maximum likelihood estimation method, calculating the posterior probability distribution of the terminal load at each time node in the prediction time domain; extracting the mean of the posterior probability distribution as the mean vector, calculating the covariance between each prediction node based on the variance of the posterior probability distribution and the kernel function, and constructing a joint distribution matrix as the load prediction covariance matrix.

[0010] This invention constructs a covariance matrix by using the variance of the posterior probability distribution output by a Gaussian process regression model and the kernel function. This matrix characterizes the temporal correlation and fluctuation range between prediction time points, providing reliable probabilistic boundary data for the system to cope with future load shocks.

[0011] Preferably, the construction of a single simplex or a set of simplexes as the time-varying solution space includes: multiplying the rated power of each heat source device by the device capacity attenuation coefficient calculated in conjunction with environmental conditions to obtain the available output extreme value; constructing a convex feasible region by combining each heat source under the constraint of the available output extreme value to satisfy the extreme distribution points of the current total load of the system; generating a single simplex or a set of simplexes composed of multiple local simplex elements as the time-varying solution space through convex envelope calculation, wherein any point inside the time-varying solution space is represented by a non-negative linear combination of vertices and the sum of the combination coefficients is 1.

[0012] Preferably, the construction of the comprehensive expected objective function includes: extracting the main diagonal elements of the load forecast covariance matrix, dividing them by the variance benchmark value, and performing normalization to obtain the uncertainty scalar; inputting the uncertainty scalar into an exponential time decay function as a decay rate adjustment parameter to generate the decay kernel; using the decay kernel to perform a weighted summation of the expected values ​​of the multidimensional objective function at each discrete time point, and using multi-objective preference weights to perform a scalar combination of the multidimensional objective function to obtain the comprehensive expected objective function in single-valued scalar form.

[0013] This invention combines an uncertainty scalar with an exponential time decay function to generate an asymmetric decay kernel, which automatically reduces the decision weight of nodes with later prediction times and lower confidence levels in global optimization, thus significantly enhancing the system's stability against medium- and long-term prediction errors.

[0014] Preferably, obtaining the Riemann natural gradient includes: constructing an energy balance mapping equation from the centroid coordinate sequence to the terminal load trajectory space; obtaining the Jacobian matrix as the energy balance mapping matrix by taking the partial derivative of the energy balance mapping equation; calculating the product of the transpose of the energy balance mapping matrix, the inverse of the load prediction covariance matrix, and the energy balance mapping matrix to generate the Riemann metric tensor constructed on the time-varying solution space; and multiplying the inverse matrix of the Riemann metric tensor with the Euclidean gradient to obtain the Riemann natural gradient.

[0015] Preferably, the search for the feasible allocation point with the minimum Riemann distance to the optimal desired allocation point includes: obtaining the pipeline flow limit, equipment ramp rate limit, and minimum start-up / stop time limit of the multi-heat source system; performing state locking based on the current start-up / stop status and duration of the equipment; converting the allowed output boundary, pipeline flow limit, and equipment ramp rate limit after state locking into multidimensional inequality constraint boundaries in the time-varying solution space as the hard constraints of the equipment; using the multidimensional inequality constraint boundaries to cut the time-varying solution space to obtain the closed subspace; and solving the problem in the subspace with the objective function of minimizing the quadratic distance between two points based on the Riemann metric tensor.

[0016] This invention transforms complex nonlinear constraints such as underlying physical pipeline limitations, equipment ramp-up rates, and start-up / shutdown times into geometric cutting boundaries within the simplex time-varying solution space, significantly reducing the complexity of quadratic programming solutions under hard constraints and avoiding equipment overload risks from the underlying logic.

[0017] Preferably, the output of the energy balance mapping equation is selected as the load balance residual at each time node in the prediction time domain, and the load balance residual is the difference between the effective heat supply of the terminal after considering efficiency, loss and delay and the mean vector.

[0018] Preferably, after calculating the product of the transpose of the energy balance mapping matrix, the inverse of the load prediction covariance matrix, and the energy balance mapping matrix to generate the Riemann metric tensor, the method further includes: when the Riemann metric tensor has insufficient rank, adding a regularization term consisting of the product of a small constant and the identity matrix to the Riemann metric tensor to construct a positive definite modified Riemann metric tensor.

[0019] Secondly, the present invention provides a multi-heat source heating and hot water energy efficiency distribution system, including a processor and a memory, wherein the memory stores computer program instructions, and when the computer program instructions are executed by the processor, the above-mentioned multi-heat source heating and hot water energy efficiency distribution method is implemented.

[0020] By adopting the above technical solution, the above-mentioned multi-heat source heating and hot water energy efficiency distribution method is generated into a computer program and stored in a memory so that it can be loaded and executed by a processor. In this way, a terminal device can be made based on the memory and the processor for convenient use.

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

[0022] This invention constructs an objective function encompassing multiple dimensions such as energy efficiency, economy, and environment, and uses a Gaussian process regression model to model multi-source operating data, generating the mean vector and covariance matrix of terminal load forecasts. This achieves a joint characterization of load change trends and forecast uncertainties. Based on the health status of heat source equipment, the actual available output extrema are determined, and a centroid coordinate simplex is constructed as a time-varying solution space, effectively reducing the complexity of the optimization boundary under multiple constraints. Simultaneously, the time decay kernel is dynamically adjusted based on the uncertainty scalar generated from the load forecast covariance matrix, forming a comprehensive expected objective function, effectively improving the stability and adaptability of the expected optimization process.

[0023] Furthermore, this invention constructs a Riemann metric tensor through energy balance mapping projection and obtains the Riemann natural gradient accordingly. This successfully transforms the complex output allocation optimization problem into a descent search within the manifold space, significantly improving the computational efficiency and solution stability of the system's output allocation. Within the feasible subspace satisfying the underlying hard constraints of the equipment, the invention further searches for and determines the feasible allocation point with the minimum Riemann distance to the optimal desired allocation point, and then executes it. This improves the operational safety of the underlying equipment and the stability of system scheduling, thereby ensuring the comprehensive optimal operation of multi-heat source heating and hot water systems, achieving economic efficiency, energy saving, and low emissions. Attached Figure Description

[0024] Figure 1 This is a flowchart of the multi-heat source heating and hot water energy efficiency distribution method in this invention; Figure 2 This is a schematic diagram of Gaussian process load prediction in this invention; Figure 3 This is a schematic diagram of the simplex solution space and extreme value constraint boundary of the multi-heat source simplex in this invention. Detailed Implementation

[0025] The technical solutions of the present invention will be clearly and completely described below with reference to the accompanying drawings of the embodiments of the present invention. Obviously, the described embodiments are some embodiments of the present invention, but not all embodiments.

[0026] This invention discloses a method for multi-heat source heating and hot water energy efficiency distribution, referring to... Figure 1 This includes steps S1-S4: S1. Construct a multidimensional objective function and generate load forecasts.

[0027] In an optional embodiment, real-time operating parameters and environmental conditions of the multi-heat source system are obtained, and a multi-dimensional objective function representing energy efficiency, economy and environment is constructed. A Gaussian process regression model is used to process historical operating data and real-time weather forecast data to generate the mean vector of terminal load in the prediction time domain and the load prediction covariance matrix.

[0028] The Modbus-TCP protocol is used to call the PyModbus library to read the supply water temperature, return water temperature, flow rate, and compressor operating frequency of each heat source as real-time operating parameters from the programmable logic controller. The Requests library is used to call a third-party meteorological RESTful application programming interface to obtain the outdoor dry-bulb temperature, relative humidity, and solar radiation intensity within the forecast time domain as environmental conditions. Energy efficiency, economic, and environmental sub-objective functions are constructed separately. The energy efficiency sub-objective function is set as the total system heat supply divided by the total system power consumption (electricity or gas calorific value); the economic sub-objective function is set as the sum of the energy consumption of each heat source multiplied by the real-time time-of-use electricity price or natural gas price; the environmental sub-objective function is set as the sum of the energy consumption of each heat source multiplied by the corresponding energy carbon emission factor. For gas-fired boilers, the gas calorific value is calculated based on the lower heating value of the gas and the gas consumption; for air-source heat pumps, ground-source heat pumps, or electric heating sources, the electric power is directly used to calculate the coefficient of performance (COP), or converted according to a preset primary energy conversion factor to uniformly calculate the system's primary energy utilization efficiency. To ensure a consistent minimization optimization direction in subsequent descent searches, the energy efficiency sub-objective function is transformed into an energy efficiency cost term when constructing the multidimensional objective function for optimization. The energy efficiency cost term is either the reciprocal of the energy efficiency sub-objective function or a dimensionless value obtained through monotonically decreasing normalization. The economic and environmental sub-objective functions maintain the cost-oriented objective direction. Therefore, higher energy efficiency corresponds to a smaller energy efficiency cost term, lower economic costs, and lower environmental penalties, resulting in smaller corresponding objective values. All three types of objectives exhibit a consistent optimization direction in the descent search.

[0029] After dimensionless normalization of the three sub-objective functions, a multi-dimensional objective vector or matrix is ​​formed in the order of energy efficiency cost term, economic objective sub-function, and environmental objective sub-function. In the subsequent construction of the comprehensive expected objective function, a scalarized combination is performed using the prediction time-domain decay kernel and multi-objective preference weights. Specifically, the multi-dimensional objective function uses a weighted combination of the normalized values ​​of the energy efficiency cost term, economic objective sub-function, and environmental objective sub-function, instead of directly adding the original energy efficiency value (where higher energy efficiency is always better) to the cost-type objective. A Gaussian process regression model is constructed, using timestamp environmental conditions from historical operating data as input features and terminal heat load as a label for model training. The covariance function is a combination of the Matern kernel function and the white noise kernel function to characterize the smooth changes and random fluctuations of the load. After model training, real-time weather forecast data is input into the trained Gaussian process regression model. The expected load values ​​for each time step in the prediction time domain are obtained through forward calculation of the model and arranged in chronological order to generate a load prediction mean vector. The model's covariance calculation method is called to output the variance and covariance between each prediction time step, which constitute the load prediction covariance matrix.

[0030] In an optional embodiment, a multidimensional objective function representing energy efficiency, economy, and environment is constructed. Specifically, supply and return water temperatures, pipeline flow rates, and outdoor meteorological data are acquired, and the real-time thermal efficiency of each heat source is calculated in conjunction with fuel quality measurement parameters. The ratio of the total system heat supply to the total calorific value of consumed fuel is used as the energy efficiency objective sub-function. The weighted sum of fuel consumption costs, equipment depreciation costs, and start-up / shutdown penalty costs of each heat source is used as the economic objective sub-function. The greenhouse gas emissions and pollutant emissions generated by the operation of each heat source are multiplied by the corresponding environmental penalty coefficients and summed to form the environmental objective sub-function. The reciprocal of the energy efficiency objective sub-function, the economic objective sub-function, and the environmental objective sub-function are dimensionlessly normalized and then combined to form the multidimensional objective function.

[0031] When calculating the real-time thermal efficiency of each heat source, the supply water temperature is collected in real time by temperature sensors and electromagnetic flow meters deployed in the primary pipeline network. The preferred range is 60℃ to 95℃, for example, a supply water temperature of 85℃. The return water temperature is also collected in real time, with a preferred range of 40℃ to 60℃, for example, a return water temperature of 50℃. The pipeline flow rate is also collected in real time, for example, a pipeline flow rate of 1200 m³ / h. 3 The total real-time heating capacity is calculated using thermodynamic formulas, taking into account the ambient temperature collected by an outdoor weather station (e.g., -5°C). Fuel quality parameters include the lower heating value of natural gas, with an preferred range of 34 MJ / m³. 3 Up to 36 MJ / m 3 The lower heating value is taken as 35.5 MJ / m. 3The total calorific value of fuel consumed is calculated by combining the real-time consumption measured by the fuel consumption meter. The ratio of the total heat supply to the total calorific value of fuel consumed is the energy efficiency target sub-function. The function value of the energy efficiency target sub-function is between 0.85 and 1.05, for example, 0.95.

[0032] Before proceeding to the comprehensive optimization objective, the reciprocal of the energy efficiency objective sub-function is taken and normalized to obtain the energy efficiency cost term. For example, when the energy efficiency objective sub-function is 0.95, the energy efficiency cost term is: The values ​​obtained after minimax normalization are combined with economic and environmental objectives to form an optimization objective where smaller values ​​are better.

[0033] For the economic objective subfunction, fuel consumption cost is derived by multiplying the unit price of natural gas by the consumption rate, for example, a unit price of natural gas of 3.2 yuan per cubic meter. Equipment depreciation cost is amortized to the hour using the straight-line depreciation method, for example, 15 yuan per hour. Start-up and shutdown penalty cost is derived by multiplying a single start-up and shutdown cost constant by the number of start-ups and shutdowns per unit time, for example, a single start-up and shutdown cost constant of 500 yuan per start-up and shutdown. Fuel consumption cost, equipment depreciation cost, and start-up and shutdown penalty cost are weighted and summed in a weighted ratio, for example, 0.6:0.2:0.2. The environmental objective subfunction utilizes a Continuous Emission Monitoring System (CEMS) to obtain CO2 and NO2 data. x The real-time emission quality of SO2, where the environmental penalty coefficient for CO2 is taken as 0.05 yuan per kilogram, and NO... x The environmental penalty coefficient is set at 15 yuan per kilogram, and the environmental penalty coefficient for SO2 is set at 10 yuan per kilogram. These are multiplied separately and then summed to reflect the overall environmental cost; for example, the total environmental cost is 450 yuan per hour. In constructing the multidimensional objective function, the minimax normalization method is used to map the reciprocal of the energy efficiency objective sub-function, the economic objective sub-function, and the environmental objective sub-function to... The dimensionless interval is used to form a multidimensional objective function matrix, which is arranged according to the time nodes within the prediction time domain. The multi-objective preference weight vector, for example, has an energy efficiency weight of 0.4, an economic weight of 0.4, and an environmental weight of 0.2, which is used for scalar combination in the subsequent stage of constructing the comprehensive expected objective function.

[0034] In an optional embodiment, a Gaussian process regression model is used to process historical operating data and real-time weather forecast data to generate a mean vector of terminal load and a load prediction covariance matrix in the prediction time domain. Specifically, environmental meteorological sequences and load sequences within the historical operating cycle are extracted as training sets, and the hyperparameters of the Gaussian process regression model are optimized using the maximum likelihood estimation method. Real-time weather forecast data is input into the trained Gaussian process regression model to calculate the posterior probability distribution of terminal load at each time node in the prediction time domain. The mean of the posterior probability distribution is extracted as a mean vector, and the covariance between each prediction node is calculated based on the variance of the posterior probability distribution and the kernel function. A joint distribution matrix is ​​then constructed as the load prediction covariance matrix.

[0035] Extract environmental meteorological sequences from at least one complete heating season in the past, such as the last 3000 hours, including ambient temperatures ranging from -15°C to 10°C, relative humidity from 30% to 80%, wind speeds from 0 to 10 m / s, and solar radiation intensity, and synchronize them with the actual heat load sequence for the same period. This load range is determined based on the system's installed capacity and is set to 10MW to 38MW in the three heat source simulation scenarios described later. The construction dimension is... ,For example The training set includes historical features. When training the Gaussian process regression (GPR) model, a Matern kernel function superimposed with a white noise kernel function is preferred as the covariance function to describe the nonlinear mapping and random disturbance relationship between meteorological conditions and load. The squared exponential kernel function can be used as another alternative implementation. In the same embodiment, after selecting one main kernel function, consistency is maintained between the training and prediction phases. The white noise kernel function is used to describe measurement noise. The gradient of the marginal likelihood function is calculated using the maximum likelihood estimation method, and unconstrained optimization is performed. Hyperparameters such as the length scale (optimized value is 2.5), signal variance is 15, and noise variance is 0.5 are iteratively converged to the optimal solution.

[0036] The system acquires future forecast time domain data provided by meteorological departments, such as forecasting the next 24 hours with a time resolution of 1 hour, generating real-time weather forecast data for 24 discrete time points. This data is then fed into the trained GPR model as input vectors. Based on the conditional probability theorem, the posterior probability distribution of the output variable is calculated. The expected mean is extracted at each time point and arranged chronologically to form a dimension... The mean vector is calculated, for example, the mean at each node fluctuates between 20MW and 35MW, with a certain node having a mean of 28.5MW. The covariance between each prediction time node is calculated using the inner product of the prior kernel function matrix and the test data, and the variance is extracted as the diagonal element to construct a positive definite and symmetric matrix. The load forecast covariance matrix represents the temporal correlation and uncertainty level of the forecast load throughout the entire forecast time domain.

[0037] Reference Figure 2 This paper presents the load forecast mean and corresponding confidence intervals within the forecast time domain. The load forecast mean exhibits a smooth, nonlinear trend with fluctuations over time, reflecting the dynamic characteristics of terminal load influenced by meteorological conditions and random energy consumption behavior during the day. The confidence intervals, enveloping the load forecast mean, represent the uncertainty boundaries at each forecast time point. As the forecast time domain progresses, the confidence intervals gradually widen, indicating that the longer the forecast time, the higher the uncertainty. The representation based on the joint probability distribution enables the constructed multidimensional objective function to fully measure the potential impact of future load fluctuations on system scheduling, providing reliable input features for subsequent uncertainty-based adjustment of the time decay kernel, thereby improving the stability and adaptability of the expected optimization process.

[0038] S2. Construct the time-varying solution space and generate the expected function.

[0039] In an optional embodiment, based on the assessment of the available power output extreme value under the constraint of the heat source health status, the extreme distribution point that satisfies the total system load under the constraint of the available power output extreme value is used as the vertex. A simplex or set of simplexes with centroid coordinates is constructed as a time-varying solution space, and the available power output extreme value is used as the constraint boundary. Each point in the time-varying solution space is uniquely mapped to the power output distribution scheme. Based on the load prediction covariance matrix, the uncertainty scalar of the corresponding time node is generated. The attenuation kernel is adjusted to weight the components of the multidimensional objective function at each time point in the prediction time domain, and a comprehensive expected objective function representing the global expected performance is constructed.

[0040] An isolated forest algorithm is used to detect anomalies in the historical operating data of each heat source, extracting the maximum heating power under normal operating conditions as the baseline capacity value. A health degradation coefficient is calculated based on the current equipment operating years and maintenance records. Multiplying the baseline capacity value by the health degradation coefficient yields the extreme available output value of each heat source. When calculating extreme load distribution points at each time step based on the total system load provided by the mean vector, not only are empirical rules prioritizing load sharing by a single heat source or the minimum number of heat sources employed, but an active constraint enumeration method is further used to generate true extreme load distribution points. Specifically, for Taiwan heat source, establish power variable to Based on the principles of system load conservation and equipment output constraints, constraint equations for power variables are constructed as follows: ;

[0041] in, For the first The power variable of the heat source; The number of heat sources; This is the index number of the heat source; The total system load for the current control cycle or the first load in the forecast time domain. Forecast total load at each time point; For the first The maximum available output of the heat source at the corresponding time point or under the current control cycle.

[0042] enumerate For each combination of power variables at the lower or upper bounds, calculate the remaining power variable and retain the points that satisfy the upper and lower bound constraints as extreme distribution points. This avoids missing the true poles determined by the total load equality constraint and the upper and lower bounds.

[0043] After extracting the extreme allocation points and forming a matrix, the convex feasible region satisfying the total load balance and available power extreme value constraints is first obtained through the convex envelope algorithm. When the convex feasible region is composed of affine independent poles with one more degree of freedom, it is directly treated as a single simplex with barycentric coordinates. When the convex feasible region forms a complex convex polyhedron due to upper and lower limits, Delaunay partitioning or other convex polyhedron partitioning algorithms are used to divide the convex feasible region into multiple local simplex elements, and the number of each local simplex element is recorded. The convex envelope algorithm refers to calculating the minimum convex geometric boundary that can completely enclose a set of discrete points, thereby forming a convex feasible region representing the feasible allocation area. The unique mapping of each point in the time-varying solution space to the power allocation scheme means that the geometric points in the feasible region uniquely correspond to a power allocation scheme in the actual power coordinates. For the case of being divided into multiple local simplex elements, the barycentric coordinates need to be used together with the simplex element number. The barycentric coordinates are unique within the selected local simplex, and the barycentric coefficients spanning multiple simplex elements are not required to be globally unique. It should be noted that the simplex solution space of the barycentric coordinates can be a single simplex directly formed by affine independent poles with one more degree of freedom, or it can be a set of multiple local simplex elements formed by partitioning a complex convex feasible region. In the latter case, the barycentric coordinates and the local simplex element numbers together determine the unique actual output distribution scheme.

[0044] Arbitrary points within a selected local simplex are represented as linear combination weights of the vertices of the local simplex, with these weights serving as the local barycentric coordinates. For a convex feasible region composed of multiple local simplexes, points in the time-varying solution space are represented by both the local simplex number and the local barycentric coordinates. The main diagonal elements of the load prediction covariance matrix are extracted to form a prediction variance sequence for each time node within the prediction time domain. The prediction variance of each time node is divided by a preset variance benchmark to obtain the uncertainty scalar for that time node, eliminating the use of the integer matrix trace as the uncertainty scalar for a single time node. The integer matrix trace can be used as a monitoring quantity for the overall uncertainty in the entire prediction time domain but is not directly involved in the calculation of the decay kernel for each time node. An exponential decay function is constructed as the decay kernel, with the decay rate set as the product of the decay coefficient, the corresponding time node uncertainty scalar, and the normalized time position. This ensures that the decay rate is faster for time nodes with later prediction times and higher uncertainty. The multidimensional objective function value corresponding to each time step in the prediction time domain is multiplied by the exponential decay kernel function value of the corresponding time step, and the weighted objective values ​​of each time point are summed using a summation function to obtain the comprehensive expected objective function.

[0045] In an optional embodiment, a time-varying solution space is constructed using the extreme distribution point where each heat source combines to satisfy the total system load under the constraint of the available output extreme value as the vertex, with the barycentric coordinate single simplex or simplex set as the vertex, and the available output extreme value as the constraint boundary, and each point in the time-varying solution space uniquely maps the output distribution scheme. Specifically, the operating status, cumulative operating time, and fault alarm records of each heat source device are read, and the capacity attenuation coefficient of each device is calculated in combination with environmental conditions. The rated power of each heat source device is multiplied by the corresponding capacity attenuation coefficient to obtain the available output extreme value of each heat source at the current moment. In the multidimensional Euclidean space, the extreme distribution point that satisfies the current total load of the system under the constraint of the available output extreme value of each heat source is taken as the extreme distribution point. A convex feasible region is constructed based on the extreme distribution point. A single simplex or a set of simplex elements composed of multiple local simplex elements is generated by convex envelope calculation. The single simplex or the set of simplex elements is taken as the time-varying solution space. Any point inside the time-varying solution space is represented as a non-negative linear combination of vertices, and the sum of the combination coefficients is 1. The combination coefficients constitute the centroid coordinates and uniquely correspond to an output distribution scheme that satisfies the load balance. The available output extreme value constitutes the upper limit boundary of the limiting distribution ratio as the constraint boundary.

[0046] The equipment data acquisition and monitoring control system reads the underlying parameters of each heat source device. These parameters include fault record characteristic parameters such as a burner load rate of 75%, a cumulative full-load operating time of 30,000 hours, and the thickness of scale buildup on the heat exchanger. Using the Weibull distribution model or empirical degradation formula, the underlying parameters are combined with the current ambient temperature to calculate the capacity decay coefficient of each heat source device. The capacity decay coefficient is preferably between 0.8 and 1; for example, the decay coefficient of gas-fired boiler No. 1 due to scaling is 0.88. The rated heating power of each heat source is extracted; for example, the rated powers of heat sources No. 1 to No. 3 are 20MW, 15MW, and 10MW, respectively. The rated power is multiplied by the corresponding capacity decay coefficient to obtain the current limited available output limit; for example, the current limit for heat source No. 1 is 17.6MW, for heat source No. 2 it is 14.25MW, and for heat source No. 3 it is 9.8MW. The output limit constitutes the upper limit of multi-heat source scheduling.

[0047] In containing Taiwan heat source In Euclidean space, for example, the current total load demand is set to 35MW. First, it is determined whether the current total load demand is less than or equal to the sum of the extreme values ​​of available output from each heat source. If the current total load demand exceeds the sum of the extreme values ​​of available output from each heat source, it is determined that there is no feasible allocation point that satisfies the total load at the current moment, and a backup heat source access, peak shaving control, or load alarm process is triggered; the system does not enter simplex optimization. When the feasibility condition is met, an affine independent extreme value allocation point in multidimensional space is determined by finding extreme output conditions that satisfy the total load conservation constraint and where no heat source exceeds its available output extreme value. For example, in a system with three heat sources, an extreme value allocation point might be where heat source 1 has an output of 17.6MW, heat source 2 has an output of 14.25MW, and heat source 3 has an output of 3.15MW. In the specific enumeration, two of the heat sources are fixed at 0 or their corresponding available output extreme value boundaries. The output value of the third heat source is calculated, and it is determined whether the output value of the third heat source falls within the allowable range. Repeat the verification process for all boundary combinations to obtain a complete set of extreme allocation points.

[0048] A convex feasible region is constructed using the convex envelope algorithm based on the set of all identified extreme distribution points, thus creating a convex feasible region that satisfies the current total load constraint of the system. If the convex feasible region is a single two-dimensional triangular cross section, it is directly used as the simplex time-varying solution space. If the convex feasible region is a quadrilateral, pentagon, or higher-dimensional convex polyhedron, it is divided into multiple non-overlapping local simplex elements, and the local simplex number and local centroid coordinates are used together as optimization variables. The power distribution scheme is based on the actual power coordinates, so each geometric point uniquely corresponds to a power distribution scheme that satisfies load balance. The coordinates of each spatial point inside and on the boundary of a single local simplex can be linearly represented by the centroid coordinates of the corresponding local simplex vertex set, where the non-negative combination coefficients of the centroid coordinates add up to 1. For complex convex feasible regions formed by splicing multiple local simplexes, they are represented by the local simplex number and the corresponding centroid coordinates. This mapping ensures that each set of centroid coordinates within the simplex satisfies the load balance equation constraint of the system, while the external geometric boundary is constructed from the available output extrema, avoiding the generation of invalid allocation schemes.

[0049] Reference Figure 3 This diagram illustrates the feasible allocation region under a specific total load constraint. The extreme values ​​of available output from heat source 1 and heat source 2 constitute the physical upper bound of the system scheduling. The diagonal line connecting them represents the output combination trajectory of the two heat sources under the condition of satisfying total load conservation. Any point within the simplex subspace region represents an effective output allocation scheme that simultaneously satisfies total load balance and equipment extreme value constraints. Points located outside the simplex subspace region will trigger equipment overload constraints or fail to meet load demands. The simplex subspace confines the complex output allocation optimization process to a closed and safe geometric region, improving the feasibility and computational efficiency of subsequent searches.

[0050] In an optional embodiment, for the uncertainty scalar generated based on the load forecast covariance matrix at corresponding time points, the attenuation kernel is adjusted to weight the components of the multidimensional objective function at each time point in the prediction time domain, constructing a comprehensive expected objective function representing the global expected performance. Specifically, the main diagonal elements of the load forecast covariance matrix are extracted as the prediction variance at each time point in the prediction time domain; the prediction variance at each time point is divided by the variance benchmark value and normalized to obtain the uncertainty scalar at the corresponding time point; the uncertainty scalar is input into the exponential time decay function as the decay rate adjustment parameter, generating an asymmetric time decay kernel as the attenuation kernel; the expected values ​​of the multidimensional objective function at each discrete time point in the prediction time domain are weighted and summed using the attenuation kernel as the weight of each time point component, and the multi-objective preference weights are used to scalarize and combine the multidimensional objectives to obtain a comprehensive expected objective function in single-valued scalar form.

[0051] From the aforementioned construction, for example In the load forecast covariance matrix, the 24 diagonal elements are extracted from the main diagonal. These diagonal elements represent the marginal forecast variance of the terminal load at each time point from hour 1 to hour 24 within the forecast time domain; for example, the variance for hour 1 is 1.2 MW. 2 In the 24th hour, the variance increased to 8.5MW due to the longer forecast period. 2 The extracted variance values ​​are divided by a preset fixed variance benchmark value obtained from the statistical analysis of the overall load variance of historical heating seasons; the preferred range for the preset fixed variance benchmark value can be set to 5MW. 2 Up to 10MW 2 The example value is 8MW. 2 Normalization is performed using a linear scaling operation to obtain a dimensionless uncertainty scalar sequence for each time node, ranging from 0.15 to 1.06. This dimensionless uncertainty scalar sequence represents the difference in the model's confidence level for each prediction node. The trace of the load prediction covariance matrix is ​​the sum of all prediction variances. The entire matrix trace is used only to evaluate the overall uncertainty level across the entire prediction time domain and is not used as the decay rate input for a single time node, thus maintaining a one-to-one correspondence between the uncertainty scalar and the prediction time nodes.

[0052] Based on the theory of exponential decay, an exponential time decay function is constructed, which reflects the physical law that the prediction confidence gradually decreases over time. Substituting the calculated uncertainty scalar sequence into the exponential time decay function, the formula is as follows:

[0053] in, For the first The decay kernel value at each time point; The attenuation coefficient is determined by the set through historical verification; For the first An uncertainty scalar at each time point; To normalize the time position, , This represents the sequence number of the current prediction time point. To predict the total number of steps in the time domain.

[0054] The aforementioned adjustment mechanism results in smaller weights for nodes with later prediction times and higher uncertainties. The attenuation kernel matrix is ​​multiplied by the expected economic, energy efficiency, and environmental benefits of the multi-dimensional objective function at each time point within the prediction time domain, and then a weighted summation is performed. Combined with the multi-objective preference weight vector set by the heating manager (e.g., setting energy efficiency preference to 0.5, economic preference to 0.3, and environmental preference to 0.2), a scalar inner product is calculated on the weighted multi-dimensional objective function, transforming the multi-objective vector matrix into a continuous, smooth, single-valued scalar comprehensive expected objective function.

[0055] S3. Construct the metric tensor and search for the optimal assignment point.

[0056] In an optional embodiment, a Riemann metric tensor is constructed by projecting the inverse of the load prediction covariance matrix onto the time-varying solution space through the energy balance mapping matrix. The Riemann metric tensor is then used to adjust the original Euclidean gradient obtained by taking the partial derivative with respect to the comprehensive expectation objective function to obtain the Riemann natural gradient. The comprehensive expectation objective function is then subjected to a descent search and converges to the optimal expectation allocation point.

[0057] An energy balance mapping equation is constructed from the simplex centroid coordinate space to the actual heat source power space and the terminal heating response space. This energy balance mapping equation does not simply sum the output allocation points that already satisfy total load conservation into a single total load scalar. Instead, it describes the impact of free centroid coordinate disturbances on the effective terminal heating response, network heat transfer loss, equipment efficiency correction, and load balance residuals in the prediction time domain. Therefore, even if all candidate allocation points satisfy nominal total load conservation, different heat source combinations will still result in different terminal heating response trajectories due to differences in efficiency, transmission loss, thermal inertia, and delay parameters; the Jacobian matrix will not degenerate into a zero matrix.

[0058] The first derivative matrix of the energy balance mapping equation is obtained using the Jacobian matrix solver, generating the energy balance mapping matrix. The information matrix is ​​obtained by inverting the load forecast covariance matrix using an inverse function. A small regularization term, a tiny constant, is added to the load forecast covariance matrix before inversion. The product of the identity matrix, where The constant is a small constant greater than 0 to ensure numerical stability. The Riemann metric tensor in the time-varying solution space is obtained by multiplying the transpose of the energy balance mapping matrix by the information matrix and then by the energy balance mapping matrix itself. The partial derivative vector of the synthesized expectation objective function with respect to the simplex coordinate variables is obtained using the gradient solution function; this vector is the original Euclidean gradient. The Riemann natural gradient is obtained using the linear equations solution function, with the modified Riemann metric tensor (with added regularization terms) as the coefficient matrix and the original Euclidean gradient as the constant term vector; the Riemann natural gradient method eliminates the influence of parameter space distortion. The Riemann gradient descent algorithm is used on the simplex manifold to set the initial simplex coordinates and learning rate. In each iteration, the current barycentric coordinates are subtracted from the product of the learning rate and the Riemann natural gradient to obtain temporary updated coordinates; subsequently, simplex projection or retraction mapping operations are performed on the temporary updated coordinates so that they again satisfy the constraint that the barycentric coordinates are non-negative and the sum of the combined coefficients is 1. If the temporarily updated coordinates exceed the current local simplex, the temporarily updated coordinates are projected onto the boundary of the current local simplex, or repositioned to an adjacent local simplex cell based on the actual power coordinates and the local simplex number is updated. The learning rate employs backtracking search or an adaptive step size within a preset upper limit to avoid oscillations near the boundary caused by a fixed learning rate. The system stops iterating when the change in the comprehensive expected objective function over two consecutive iterations is less than a preset small threshold; at this point, the coordinates have converged to a first-order stable optimal expected allocation point that satisfies the current model and constraints.

[0059] In an optional embodiment, a Riemann metric tensor is constructed by projecting the inverse of the load prediction covariance matrix onto the time-varying solution space using the energy balance mapping matrix. The Riemann natural gradient is obtained by adjusting the original Euclidean gradient obtained by taking the partial derivative of the comprehensive expectation objective function using the Riemann metric tensor. Specifically, an energy balance mapping equation is constructed from the centroid coordinate sequence of each time node in the prediction time domain to the terminal load trajectory space of each time node in the prediction time domain, and the partial derivative of the energy balance mapping equation is obtained to obtain the Jacobian matrix as the energy balance mapping matrix. The information matrix is ​​obtained by inverting the load prediction covariance matrix. The product of the transpose of the energy balance mapping matrix, the information matrix, and the energy balance mapping matrix is ​​calculated to generate the Riemann metric tensor constructed on the time-varying solution space. The original Euclidean gradient is obtained by taking the first-order partial derivative of the comprehensive expectation objective function. The Riemann natural gradient is obtained by multiplying the inverse matrix of the Riemann metric tensor with the original Euclidean gradient.

[0060] Based on the first law of thermodynamics, an energy balance mapping equation is established from the centroid coordinate sequence in simplex space to the trajectory of the system's terminal heat load response. This equation describes the impact of different heat source output combinations on the terminal heat supply after corrections for equipment efficiency, pipeline losses, and transmission delays. The energy balance mapping equation incorporates parameters such as a pipeline heat loss rate (e.g., set to 3%) and a heat network transmission delay time (e.g., a hydraulic delay constant set to 1.5 hours). The output of the energy balance mapping equation is selected as the load balance residual at each time point within the prediction time domain, the effective terminal heat supply contribution, or the heat response after corrections for pipeline losses and equipment efficiency, rather than simply the sum of the outputs of each heat source. If the output uses the load balance residual, a load balance residual is constructed to accurately assess the deviation level between the output allocation scheme and the load prediction mean.

[0061] in, For the first Load balance residuals at each time point; The centroid coordinates are free variables; To take into account efficiency, losses, and delays in the effective heat supply to the terminal; To predict the average load.

[0062] By taking the first-order partial derivatives of the energy balance mapping equation in each free coordinate direction of the current runtime variation space, the Jacobian matrix representing the local geometric rate of change is derived, and the Jacobian matrix serves as the energy balance mapping matrix. In this invention, a rolling optimization method is adopted, meaning that in each control cycle, only one set of centroid coordinates at the current moment is used as the free decision variable, and subsequent coordinates in the prediction time domain are recursively generated by combining the current coordinates with the ramp rate and load prediction. Therefore, in a three-heat source system with 2 degrees of freedom and a prediction step size of 24, the dimension of the Jacobian matrix is... The Riemannian metric tensor has a dimension of If a full-predictive time-domain joint scheduling method is adopted, then an independent centroid coordinate free variable is set for each time node. The dimension of the free variable for the 3 heat sources 24-step system is... In this case, the Jacobian matrix and Riemannian metric tensor should be expanded accordingly to a matrix structure matching the 48-dimensional free variables. The two implementations do not mix dimensions.

[0063] For example, the distribution of load forecast uncertainty The load prediction covariance matrix is ​​inverted to extract the information matrix representing the deterministic probability distribution. The transpose of the energy balance mapping matrix, the information matrix, and the product of the energy balance mapping matrix are calculated sequentially. A quadratic transformation projects the statistical metrics in the high-dimensional terminal load uncertainty space onto the reduced-dimensional time-varying simplex solution space, generating a Riemann metric tensor matrix with symmetric positive semi-definite properties. To ensure the stability of the solution in the descent direction, based on the geometric principle that metric tensors in manifold optimization must be positive definite matrices, when the Riemann metric tensor matrix has insufficient rank or is close to singular, a regularization term is added to the Riemann metric tensor to construct a positive definite modified Riemann metric tensor:

[0064] in, To correct the Riemannian metric tensor; This is the energy balance mapping matrix; This is the transpose of the energy balance mapping matrix; For information matrix; It is a regular term; It is a small constant greater than 0; It is an identity matrix.

[0065] Using automatic differentiation techniques, partial derivatives are taken with respect to each independent barycentric coordinate variable for the single scalarized comprehensive expectation objective function, yielding the original Euclidean gradient vector in standard Euclidean space. For example, the gradient vector is... Based on the mapping relationship between the natural gradient and the Euclidean gradient in Riemannian geometry, the Riemannian natural gradient is derived and constructed:

[0066] in, This is the modified Riemannian metric tensor after adding regularization terms; For the Riemann natural gradient; This is the original Euclidean gradient vector.

[0067] By solving The Riemann natural gradient was obtained. This process linearly corrects the Euclidean gradient, which originally ignored the spatial curvature, to the Riemannian natural gradient, which conforms to the actual intrinsic geometry and probability distribution uncertainties of the Riemannian manifold. The correction process guides the optimization search path, reducing the oscillations caused by variable coupling and achieving stable convergence towards the desired optimal solution. The optimal expected allocation point refers to the optimization result obtained under given prediction models, objective functions, constraints, and convergence criteria. Without further proof of the convexity of the objective function and the constraint set, it is not limited to interpreting it as the mathematically optimal solution.

[0068] S4. Search for feasible allocation points and issue power setting values.

[0069] In an optional embodiment, within a subspace that satisfies the device's hard constraints, a feasible allocation point with the smallest Riemann distance to the optimal desired allocation point is searched, and the feasible allocation point is inversely mapped to the power setpoints of each heat source and then executed.

[0070] The minimum start-up power and maximum load increase / decrease rate of the heat source equipment are collected as hard constraints. For the minimum start-up and shutdown time limit, the current start-up / shutdown status of the equipment and the duration of continuous operation or shutdown are first read. If the equipment has not yet met the minimum operating time or minimum shutdown time, the equipment status is locked as "on" or "off" within the current control cycle, and the allowable output range of the equipment is adjusted accordingly. For example, the allowable output range of the shut-down equipment is locked to 0, and the allowable output range of the operating equipment is locked to be no less than the minimum stable output. After the status is locked, the pipeline flow limit, ramp rate limit, and available output limit of the remaining continuous output variables are transformed into linear or nonlinear inequality constraints in simplex space. If the current control cycle allows changes to the start-up and shutdown status, the optional start-up and shutdown combinations are first determined through state enumeration or mixed integer quadratic programming, and then the continuous output allocation problem is solved for each combination.

[0071] Hard constraints, along with available power extrema, are transformed into a system of polyhedral linear inequalities or nonlinear boundary constraints in simplex space, forming a subspace. If the feasible subspace obtained after cutting through the equipment's hard constraints is empty, a protection process is triggered, including backup heat source access, load reduction, constraint relaxation alarm, or maintaining safe power output from the previous control cycle. The system does not execute the conventional Riemann distance projection solution. The distance objective function is constructed by multiplying the transpose of the difference between the coordinates of the point to be solved and the coordinates of the optimal desired allocation point by the Riemann metric tensor, and then by the difference in coordinates. For continuous power allocation problems where start-stop states are fixed, an optimization module is used to solve the problem using a sequential quadratic programming algorithm or a quadratic programming solver. For scheduling problems involving start-stop state selection, they are not directly treated as ordinary continuous quadratic programming problems, but rather a two-stage solution method using state enumeration, mixed integer quadratic programming, or locking the state before continuous optimization is adopted. The coordinate solutions output by the optimization algorithm are used as feasible allocation points. The coordinate vector of the feasible allocation point is multiplied by the extreme allocation point matrix using the inverse transformation matrix of the barycentric coordinates to convert it into the actual power values ​​of each heat source, thus obtaining the power setpoints of each heat source. By calling the industrial control communication library based on the TCP / IP protocol, the power setpoints of each heat source are written in parallel to the register address corresponding to the underlying distributed control system, triggering the underlying closed-loop controller to perform power adjustment.

[0072] In an optional embodiment, within a subspace satisfying the hard constraints of the equipment, a feasible allocation point with the minimum Riemann distance to the optimal desired allocation point is searched. Specifically, the pipeline flow limit, equipment ramp rate limit, and minimum start-stop time limit of the multi-heat source system are obtained. For the minimum start-stop time limit, state locking or state enumeration is first performed based on the current start-stop state and duration of the equipment. Then, the locked or enumerated allowable output boundary, pipeline flow limit, and equipment ramp rate limit are transformed into multidimensional linear inequality constraint boundaries or nonlinear inequality constraint boundaries in the time-varying solution space as hard constraints of the equipment. The time-varying solution space is then cut using the inequality constraint boundaries to eliminate regions that violate the hard constraints of the equipment, resulting in a closed subspace. Within the subspace, a quadratic programming problem or a constrained continuous optimization problem is solved with the objective function of minimizing the quadratic distance between two points based on the Riemann metric tensor, to obtain the coordinate point with the minimum Riemann distance to the optimal desired allocation point as a feasible allocation point.

[0073] The system acquires underlying hard constraint parameters in real time, specifically covering the maximum flow rate that the main pipelines of the heating network can withstand, such as a flow limit of 1500 m³ / h; the extreme values ​​of the load ramp rate for each heat source device due to thermal stress limitations, such as a maximum allowable ramp rate of 0.5 MW per minute for boiler No. 1; and the minimum continuous operating time or downtime set to avoid equipment thermal fatigue, such as a minimum requirement of 120 minutes of continuous operation. The minimum start-up and shutdown times are not directly treated as linear inequalities for ordinary continuous variables, but rather a state-variable flag is first formed based on the current state of the equipment and the duration already elapsed. For equipment with immutable states, the start-up and shutdown states are fixed within the current optimization cycle, and the output upper and lower limits formed under the fixed states are then written into the continuous inequality constraints. Using a mapping function, boundary conditions spanning different dimensions are transformed into a multidimensional linear equation system and nonlinear inequalities containing state discrimination logic, uniformly projected onto a simplex time-varying solution space represented by barycentric coordinates, constructing a polyhedral hyperplane constraint boundary representing the hard constraints of the system equipment. By using the hyperplane cutting plane method, the generated multidimensional inequality constraint boundary is used to perform Boolean intersection and polyhedral trimming on the simplex time-varying solution space, eliminating regions that violate the constraints but meet the total load but will cause equipment overload or failure.

[0074] The remaining core region after the cutting process forms a closed feasible subspace that satisfies the hard constraints of the equipment. When the optimal expected allocation point is outside this feasible region, a quadratic Mahalanobis distance objective function based on the Riemann metric tensor is constructed within the feasible subspace, using the optimal expected allocation point as the anchor point. If all hard constraints have been transformed into continuous linear constraints and the distance objective function is quadratic, a quadratic programming solver is invoked to solve the problem; if nonlinear continuous constraints exist, the sequential quadratic programming algorithm (SQP) is invoked to solve the corresponding constrained continuous optimization problem, finding the feasible allocation point that is closest to the theoretical optimum and satisfies the hard constraints under the quadratic distance defined by the Riemann metric tensor. The feasible allocation point is then inversely mapped to a specific output power setting value, for example, instructing device 1 to output 15.5MW and device 2 to output 12MW, and this is sent to the actuator.

[0075] A simulation environment was constructed containing three heat source heating systems with rated capacities of 20MW, 15MW, and 10MW. Based on actual operating data and weather forecasts from a 30-day heating season, the load range was set to fluctuate between 10MW and 38MW. At each simulation moment, the terminal load was verified to not exceed the sum of the extreme values ​​of the available output of each heat source after capacity decay. If the predicted load at any moment exceeds the total available output, the system initiates a backup heat source access or peak shaving alarm process, and this moment is not included in the regular three-heat-source allocation optimization sample. Three sets of experiments were constructed for comparative analysis. The basic control group adopted a deterministic load forecasting and scheduling strategy with only economic cost as the single objective, fixed the time weight of the forecast time domain, and performed spatial optimization based on the traditional Euclidean gradient. The ablation experimental group added a multi-dimensional objective function under the conditions of the basic control group and incorporated an uncertainty scalar adjustment time decay kernel mechanism based on the load forecast covariance matrix, but the optimization process still retained the standard Euclidean gradient algorithm. The complete solution group adopts a comprehensive set of optimization strategies, including multidimensional objective functions, time decay kernels, Riemann natural gradient correction, and hard constraint boundary search based on Riemann distance.

[0076] Various scheduling comparison data were obtained after 720 hours of continuous simulation operation. The average daily comprehensive operating cost of the basic control group was 56,000 yuan, the average total thermal efficiency remained at 88.5%, the average daily environmental penalty cost for pollutants was 12,000 yuan, and the scheduling commands generated a total of 120 optimization oscillations due to heat source coupling. In the ablation experimental group, due to the addition of multiple targets, the average daily comprehensive cost decreased to 52,000 yuan, the total thermal efficiency increased to 91.2%, and the average daily environmental penalty cost decreased to 8,500 yuan, but the number of optimization oscillations was still 85. The average daily comprehensive cost of the complete scheme group decreased to 48,000 yuan, the total thermal efficiency stabilized at 94.5%, and the average daily environmental penalty cost decreased to 6,200 yuan. Due to the intervention of the Riemann metric tensor, the number of optimization oscillations decreased to 3, and the smoothness of the scheduling commands improved.

[0077] The synergistic effect of the multidimensional objective function and the time decay mechanism breaks through the limitations of traditional single-dimensional optimization, enabling the system to achieve a more stable balance between environmental benefits, equipment energy efficiency, and economic expenditure. Simultaneously, the evaluation model reduces the error impact caused by weather forecast variance. Furthermore, the application of the Riemannian metric tensor to extract the Riemannian natural gradient corrects the tendency to ignore the geometric correlation of variables when using Euclidean gradient descent. Combined with the hard-constrained subspace search rule constructed by polyhedral pruning, the system ensures that it does not exceed the extreme limits of the physical pipeline network and equipment in the pursuit of desired performance optimization, enhancing the overall scheduling stability and the safety of underlying equipment operation.

[0078] This invention also discloses a multi-heat source heating and hot water energy efficiency distribution system, including a processor and a memory. The memory stores computer program instructions, which, when executed by the processor, implement the multi-heat source heating and hot water energy efficiency distribution method according to the invention.

[0079] The system also includes other components well known to those skilled in the art, such as communication buses and communication interfaces, the settings and functions of which are known in the art and will not be described in detail here.

[0080] It should be noted that those skilled in the art can make various modifications and improvements without departing from the inventive concept, and these all fall within the scope of protection of this invention. Therefore, the scope of protection of this patent should be determined by the appended claims.

Claims

1. A method for multi-heat source heating and hot water energy efficiency distribution, characterized in that, include: S1. Obtain real-time operating parameters and environmental conditions of the multi-heat source system, and construct a multi-dimensional objective function representing energy efficiency, economy and environment; Using a Gaussian process regression model to process historical operating data and meteorological forecast data, a mean vector and a load forecast covariance matrix for terminal load prediction are generated; S2, based on the extreme values ​​of available output assessed by the health status of heat sources, and taking the extreme distribution points where each heat source combines under the constraints of the extreme values ​​of available output to satisfy the total load of the system as vertices, a single simplex or a set of simplexes with centroid coordinates is constructed as a time-varying solution space; based on the load forecast covariance matrix, an uncertainty scalar for the corresponding time node is generated, and the attenuation kernel is adjusted to weight the components of the multidimensional objective function at each time point to construct a comprehensive expected objective. S3. Project the inverse matrix of the load prediction covariance matrix onto the time-varying solution space using the energy balance mapping matrix to construct a Riemann metric tensor. Use the Riemann metric tensor to adjust the Euclidean gradient obtained by taking the partial derivative with respect to the comprehensive expectation objective function to obtain the Riemann natural gradient. Perform a descent search on the comprehensive expectation objective function to converge to the optimal expectation allocation point. S4. In the subspace that satisfies the hard constraints of the equipment, search for the feasible allocation point with the smallest Riemann distance to the optimal expectation allocation point. Inversely map the feasible allocation point to the power setpoints of each heat source and issue them for execution.

2. The multi-heat source heating and hot water energy efficiency distribution method according to claim 1, characterized in that, The construction of the multidimensional objective function representing energy efficiency, economy, and environment includes: using the ratio of the total heat supply of the system to the total calorific value of the fuel consumed as the energy efficiency sub-objective function; using the sum of the energy consumption of each heat source multiplied by the real-time time-of-use electricity price or natural gas price as the economic sub-objective function; using the sum of the energy consumption of each heat source multiplied by the corresponding energy carbon emission factor as the environmental sub-objective function; and combining the reciprocal of the energy efficiency sub-objective function, the economic sub-objective function, and the environmental sub-objective function after performing dimensionless normalization processing to form the multidimensional objective function.

3. The multi-heat source heating and hot water energy efficiency distribution method according to claim 1, characterized in that, The process of generating the mean vector and load prediction covariance matrix for terminal load forecasting includes: inputting real-time weather forecast data into a Gaussian process regression model that optimizes hyperparameters using the maximum likelihood estimation method, calculating the posterior probability distribution of terminal load at each time node in the prediction time domain; extracting the mean of the posterior probability distribution as the mean vector; calculating the covariance between each prediction node based on the variance of the posterior probability distribution and the kernel function; and constructing a joint distribution matrix as the load prediction covariance matrix.

4. The multi-heat source heating and hot water energy efficiency distribution method according to claim 3, characterized in that, The construction of a single simplex or a set of simplexes as the time-varying solution space includes: multiplying the rated power of each heat source device by the device capacity attenuation coefficient calculated in conjunction with environmental conditions to obtain the available output extreme value; constructing a convex feasible region by combining each heat source under the constraint of the available output extreme value to satisfy the extreme distribution points of the current total load of the system; and generating a single simplex or a set of simplexes composed of multiple local simplex elements as the time-varying solution space through convex envelope calculation, wherein any point inside the time-varying solution space is represented by a non-negative linear combination of vertices and the sum of the combination coefficients is 1.

5. The multi-heat source heating and hot water energy efficiency distribution method according to claim 1 or 2, characterized in that, The construction of the comprehensive expected objective function includes: extracting the main diagonal elements of the load forecast covariance matrix, dividing them by the variance benchmark value, and performing normalization to obtain the uncertainty scalar; inputting the uncertainty scalar into an exponential time decay function as a decay rate adjustment parameter to generate the decay kernel; using the decay kernel to perform a weighted summation of the expected values ​​of the multidimensional objective function at each discrete time point, and using multi-objective preference weights to perform a scalar combination of the multidimensional objective function to obtain the comprehensive expected objective function in single-valued scalar form.

6. The multi-heat source heating and hot water energy efficiency distribution method according to claim 1, characterized in that, The step of obtaining the Riemann natural gradient includes: constructing an energy balance mapping equation from the centroid coordinate sequence to the terminal load trajectory space; obtaining the Jacobian matrix as the energy balance mapping matrix by taking the partial derivative of the energy balance mapping equation; calculating the product of the transpose of the energy balance mapping matrix, the inverse of the load prediction covariance matrix, and the energy balance mapping matrix to generate the Riemann metric tensor constructed on the time-varying solution space; and multiplying the inverse matrix of the Riemann metric tensor with the Euclidean gradient to obtain the Riemann natural gradient.

7. The multi-heat source heating and hot water energy efficiency distribution method according to claim 1, characterized in that, The search for the feasible allocation point with the minimum Riemann distance to the optimal desired allocation point includes: obtaining the pipeline flow limit, equipment ramp rate limit, and minimum start-up / shutdown time limit of the multi-heat source system; performing state locking based on the current start-up / shutdown state and duration of the equipment; transforming the allowed output boundary, pipeline flow limit, and equipment ramp rate limit after state locking into multidimensional inequality constraint boundaries in the time-varying solution space as the hard constraints of the equipment; using the multidimensional inequality constraint boundaries to cut the time-varying solution space to obtain the closed subspace; and solving within the subspace with the objective function of minimizing the quadratic distance between two points based on the Riemann metric tensor.

8. The multi-heat source heating and hot water energy efficiency distribution method according to claim 6, characterized in that, The output of the energy balance mapping equation is selected as the load balance residual at each time node in the prediction time domain. The load balance residual is the difference between the effective heat supply of the terminal after considering efficiency, loss and delay and the mean vector.

9. The multi-heat source heating and hot water energy efficiency distribution method according to claim 6, characterized in that, After generating the Riemann metric tensor by multiplying the transpose of the energy balance mapping matrix, the inverse of the load prediction covariance matrix, and the energy balance mapping matrix, the method further includes: when the Riemann metric tensor has insufficient rank, adding a regularization term consisting of the product of a small constant and the identity matrix to the Riemann metric tensor to construct a positive definite modified Riemann metric tensor.

10. A multi-heat source heating and hot water energy efficiency distribution system, characterized in that, include: A processor and a memory, wherein the memory stores computer program instructions that, when executed by the processor, implement the multi-heat source heating and hot water energy efficiency distribution method according to any one of claims 1-9.

Images