Optimization method for design parameters of double-layer film phase change energy storage greenhouse for biogas

Abstract

This invention provides a method for optimizing design parameters of a double-layer membrane phase change energy storage greenhouse for biogas, belonging to the field of greenhouse design optimization technology. This invention sets phase change temperature and membrane spacing as joint optimization variables and defines constraints to prevent spatial interference; it obtains a dataset consisting of annual heating energy consumption, life-cycle costs, and life-cycle carbon emissions through environmental boundary numerical simulation, comprehensively reflecting long-term, multi-dimensional operating conditions; it constructs a surrogate model using polynomial surface fitting to reduce the high-frequency computational cost of multi-physics fields; based on the Rayleigh-Stefan ratio, first-order partial derivatives, and normalized difference, it constructs evaluation functions for thermal mismatch factor, local topological sensitivity, and marginal benefit index, respectively, quantifying physical and economic characteristics; it integrates these evaluation functions into multi-objective optimization to construct penalized mutation and guided congestion operators, and outputs the optimal combination based on the Pareto solution set proximity, eliminating thermophysical distorted solutions and ensuring the value of engineering decision-making.

Description

Technical Field

[0001] This invention relates to the field of greenhouse design optimization technology, specifically to a method for optimizing design parameters of a biogas-using double-layer film-covered phase change energy storage solar greenhouse. Background Technology

[0002] Biogas fermentation requires maintaining a stable fermentation temperature during winter, and introducing a solar greenhouse with double-layer film covering and phase change energy storage is an effective approach. However, the design of its core parameters involves complex unsteady coupling of fluid mechanics and heat transfer. Current greenhouse design is increasingly adopting algorithmic optimization, but it faces the problem of high computational costs for multi-physics fields, and the modern agricultural evaluation system has evolved towards a multi-dimensional synergy of energy consumption, economy, and environmental protection. Conventional algorithms and physical models are only superficially integrated, and when solving complex multi-objective spaces, they suffer from low optimization efficiency and are prone to falling into local fallacies.

[0003] In the prior art, for example, CN118569104A discloses a method for optimizing greenhouse insulation parameters. This method generates multimodal data based on the insulation design scheme, inputs it into a trained temperature prediction model to obtain prediction results, and directly uses the prediction results as the objective function value. It then employs a conventional evolutionary algorithm for crossover, mutation, and iterative optimization to finally output an optimized result that meets the preset temperature.

[0004] However, the aforementioned existing technologies still have fundamental technical shortcomings in practical applications. First, the existing optimization and evaluation systems are limited to construction costs and insulation effects, failing to incorporate the dynamic costs throughout the entire life cycle of the project and the carbon emissions from construction to demolition into the joint optimization space, thus failing to reflect the comprehensive multi-dimensional performance of the greenhouse under long-term operation. Second, when using evolutionary algorithms for optimization, there is a lack of underlying physical mechanism guidance. The algorithm's crossover and mutation actions are only controlled by pure probability distributions, failing to introduce the thermophysical mismatch mechanism between natural convection and phase change material heat storage in the greenhouse space for dynamic intervention. Furthermore, the local topological gradient of the objective function surface is not considered, causing the algorithm to easily get stuck in high-gradient sensitive regions where even small variable fluctuations can lead to drastic performance degradation during iteration. This results in a large number of invalid solutions with poor anti-interference capabilities and thermophysical distortions, causing a significant waste of computing power. Finally, when the algorithm truncates and filters non-dominated solutions, the existing technology relies solely on the pure geometric crowding degree based on the Euclidean distance in multidimensional space. It lacks a quantitative consideration of the relative replacement rate between environmental reduction and cost of candidate solutions. This leads to the easy elimination of advantageous engineering solutions with high marginal benefits in areas of overcrowding, resulting in the final output set of solutions lacking accurate guidance value for actual engineering.

[0005] The information disclosed in the background section is only intended to enhance the understanding of the background of this disclosure, and therefore may include information that does not constitute prior art known to those skilled in the art. Summary of the Invention

[0006] The purpose of this invention is to provide a method for optimizing the design parameters of a biogas-use double-layer membrane phase change energy storage solar greenhouse, so as to solve the problems mentioned in the background art.

[0007] To achieve the above objectives, the present invention provides the following technical solution: A method for optimizing design parameters of a biogas-use double-layer membrane phase change energy storage greenhouse, comprising the following steps: Based on the greenhouse space structure model, the phase change temperature and the film covering distance are defined as joint optimization variables. The constraint range of the joint optimization variables is set by combining the preset target fermentation temperature and spatial interference constraints. Meteorological and soil parameters were extracted to calculate the initial soil temperature, which was then set together with the target fermentation temperature as the physical field boundary condition. Multiple sets of joint optimization variables were set according to the constraint range, and numerical simulation was used under the physical field boundary condition to obtain annual heating energy consumption, life cycle cost and life cycle carbon emission data, which constituted a multi-objective sample dataset. Based on a multi-objective sample dataset, a proxy model for annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions is constructed using polynomial surface fitting. Based on the population and physical field boundary conditions to be input, the ratio of the natural logarithm of the Rayleigh number of air convection to the Stefan number of phase change heat storage is calculated to construct the thermal mismatch factor evaluation function. The first-order partial derivatives of the surrogate model are extracted and the local maximum gradient is calculated to construct the local topological sensitivity evaluation function. The marginal substitution rate between the target normalized differences is calculated to construct the marginal benefit index evaluation function. A multi-objective optimization algorithm is used to iterate and optimize the surrogate model. Within the constraints of the joint optimization variables, an initial population is randomly generated. During iteration, the three evaluation functions are called. A penalty mutation operator is constructed by introducing the thermal mismatch factor and local topological sensitivity to generate a subpopulation to be screened. Then, a directional crowding operator is constructed by introducing the marginal benefit index to screen and generate a subpopulation. The process is repeated until the termination condition is met, and the Pareto solution set is output. Calculate the proximity of each solution in the Pareto solution set, and output the joint optimization variable combination with the highest proximity as the optimal joint optimization variable.

[0008] Furthermore, based on the preset target fermentation temperature, a temperature lower than the target fermentation temperature is set as the lower limit of the phase transition temperature, and a temperature not lower than the target fermentation temperature is set as the upper limit of the phase transition temperature, thus setting a constraint range for the phase transition temperature; Based on the building structure, the effective static air layer thickness threshold is set as the lower limit, and the maximum critical distance to prevent physical contact between the internal membrane and the building's external envelope is set as the upper limit. Specifically, the membrane spacing is the distance between the internal membrane on the shady side of the biogas digester and the insulation layer on the shady side.

[0009] Furthermore, the average surface temperature, annual periodic surface temperature, and annual temperature fluctuation period of the project site are extracted as meteorological parameters, and the soil thermal conductivity is extracted as a soil physical property parameter. Based on the principle of unsteady periodic heat conduction, a semi-infinite flat wall periodic heat conduction analytical solution is constructed to calculate the initial soil temperature at the preset target depth. The initial soil temperature at the target depth and the target fermentation temperature are set together as the physical field boundary conditions. Multiple sets of joint optimization variables are substituted into the established physical field model. The target fermentation temperature is used as the internal heat source boundary, and the initial soil temperature and meteorological parameters are used as the external environment boundary. Fluid dynamics and heat transfer transient coupling simulation is performed using COMSOL Multiphysics to calculate the cumulative dynamic heating load required to maintain the target fermentation temperature. Based on the principles of heat conservation and electrothermal conversion efficiency, the ratio of the cumulative dynamic heating load to the conversion parameter is calculated using the preset electrothermal constant and the energy efficiency ratio of the heating equipment as conversion parameters to obtain the annual heating energy consumption.

[0010] Furthermore, based on the principle of the time value of money comprehensive discount model in the P1-P2 economic model, the present value factor of total expenditure within the economic analysis period and the ratio of total expenditure to initial investment within the economic analysis period are calculated. The specific calculation logic is as follows: Based on the preset economic analysis period, energy price growth rate, and market discount rate, when the energy price growth rate and the market discount rate are equal, the ratio of the economic analysis period to the energy price growth rate plus one is calculated as the present value factor of the total expenditure within the economic analysis period. When the energy price growth rate and the market discount rate are not equal, the ratio of the energy price growth rate plus one to the market discount rate plus one is calculated as the economic analysis period raised to the power of the economic analysis period. The result of subtracting the power from one is then divided by the difference between the market discount rate and the energy price growth rate to obtain the present value factor of the total expenditure within the economic analysis period. Based on the preset initial investment and down payment ratio, loan term, actual repayment period, loan interest rate, ratio of annual maintenance costs to initial investment, and ratio of resale price to initial investment, calculate the difference between the initial investment and down payment ratio and multiply it by the present value factor calculated using the market discount rate for the actual repayment period, and divide it by the present value factor calculated using the loan term and loan interest rate to obtain the loan cost present value ratio. Calculate the product of the ratio of annual maintenance costs to initial investment and the present value factor calculated using the market discount rate for the economic analysis period to obtain the maintenance cost present value ratio. Calculate the quotient of the ratio of resale price to initial investment raised to the power of the market discount rate plus one for the economic analysis period to obtain the residual value present value ratio. Add the initial investment and down payment ratio, the loan cost present value ratio, and the maintenance cost present value ratio, and subtract the residual value present value ratio to calculate the ratio of total expenditure to initial investment within the economic analysis period. In calculating heating operation costs, the annual heating energy consumption is multiplied by the preset electricity price to obtain the annual heating operation costs. In the initial investment cost calculation for the renovation, the product of the sum of the film covering construction cost and the film covering unit price and the film covering area is calculated; the product of the sum of the galvanized steel pipe construction cost and the galvanized steel pipe unit price and the galvanized steel pipe length is calculated; the product of the plastering construction cost and the plastering area is calculated; and the product of the phase change energy storage material unit price and the phase change energy storage material mass is calculated. The sum of the four product results in the initial investment cost calculation for the renovation is then obtained. Based on the above calculation results, in the calculation of the total life cycle cost, the annual heating operation cost is multiplied by the present value factor of the total expenditure within the economic analysis period, and the initial investment cost of the renovation is multiplied by the ratio of the total expenditure within the economic analysis period to the initial investment. The sum of the two products in the calculation of the total life cycle cost is then used to obtain the total life cycle cost.

[0011] Furthermore, the annual heating energy consumption is multiplied by the preset economic lifespan and the grid carbon emission factor to calculate the carbon emissions during the operation phase. Extract the bill of materials from the greenhouse space construction model. Based on the preset carbon emission factors of various building materials, multiply the carbon emission factors of various building materials by their corresponding bill of materials and sum them up to calculate the carbon emissions during the production stage. Based on the preset carbon emission factors and average transportation distances of various building materials, the carbon emission factors and average transportation distances of various building materials are multiplied together with their corresponding bill of materials and summed to calculate the carbon emissions during the transportation stage. The carbon emissions during the construction phase are estimated using a proportional method, with the carbon emissions during the construction phase calculated as 4% of the total carbon emissions during the production and transportation phases, and the carbon emissions during the demolition phase calculated as 10% of the total carbon emissions during the construction phase. The carbon emissions during the operation phase, production phase, transportation phase, construction phase, and demolition phase are summed to obtain the carbon emissions data for the entire life cycle. A multi-objective sample dataset is constructed based on annual heating energy consumption, total life-cycle costs, and total life-cycle carbon emissions.

[0012] Furthermore, based on the multi-objective sample dataset, the combination of joint optimization variables in each group of the multi-objective sample dataset is used as the independent variable, and the corresponding annual heating energy consumption, life cycle cost and life cycle carbon emission data are used as the dependent variables. Based on the principle of multivariate higher-order nonlinear regression, a sixth-order bivariate polynomial surface fitting is performed to obtain the constant term, the higher-order main effect regression coefficients of each independent variable and the multivariate higher-order cross-coupling effect regression coefficients. Using the obtained regression coefficients, three sixth-order polynomial surrogate models are constructed, with phase change temperature and coating spacing as inputs and annual heating energy consumption, life cycle cost and life cycle carbon emission as outputs.

[0013] Furthermore, the entire combination of candidate joint optimization variables to be input is treated as a population, and a single joint optimization variable in the combination of candidate joint optimization variables to be input is treated as an individual in the population; Based on each individual population and physical field boundary conditions, according to the preset gravitational acceleration and the thermal expansion coefficient, kinematic viscosity and thermal diffusivity of the air between the double-layer membrane, and according to the internal target fermentation temperature and the external outdoor weather temperature, the absolute value of the difference between the two is calculated as the internal and external temperature difference. The gravitational acceleration, the thermal expansion coefficient of the air between the double-layer membrane, the internal and external temperature difference and the cube of the membrane spacing to be input are multiplied together and then divided by the product of the kinematic viscosity and thermal diffusivity of the air between the double-layer membrane to calculate the Rayleigh number characterizing the natural convection intensity of the air between the double-layer membrane. According to the preset specific heat capacity and latent heat enthalpy of the phase change material, the absolute value of the difference between the phase change temperature and the internal target fermentation temperature is calculated and multiplied by the specific heat capacity of the phase change material and then divided by the latent heat enthalpy to calculate the Stefan number characterizing the heat storage and release characteristics of the phase change material. The natural logarithm of the Rayleigh number is calculated and divided by the Stefan number to construct the thermal mismatch factor evaluation function. For three sixth-order polynomial surrogate models of annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions, the first-order partial derivatives of phase transition temperature and film spacing for each individual in the population are obtained. The two first-order partial derivatives corresponding to each surrogate model are squared, summed, and the square root is taken to obtain the spatial gradient magnitude of each surrogate model at the current variable position. The maximum value of the spatial gradient magnitudes of the three sixth-order polynomial surrogate models is selected to construct a local topological sensitivity evaluation function. The arithmetic mean of the predicted values ​​output by three sixth-order polynomial surrogate models for all individuals in the population is calculated. The population mean is obtained for the dimensions of annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions. The normalized difference of each individual relative to the population mean in the dimensions of annual heating energy consumption and life-cycle carbon emissions is calculated. The absolute values ​​of the two are added to obtain the comprehensive environmental benefit component. The normalized difference of each individual relative to the population mean in the dimension of life-cycle cost is calculated. The absolute value of the difference is taken and a preset small constant is added to obtain the cost component. The comprehensive environmental benefit component is divided by the cost component to construct the marginal benefit index evaluation function.

[0014] Furthermore, within the constraints of the joint optimization variables, a preset number of joint optimization variable combinations are randomly generated as the initial population. When using the NSGA-II algorithm for iterative solution, the generated initial population is used as the parent population to call the thermal mismatch factor evaluation function, the local topological sensitivity evaluation function, and the marginal benefit index evaluation function. When the next generation of offspring population is generated, the generated next generation of offspring population is used as the new parent population to continue calling the three evaluation functions for optimization iteration. When the NSGA-II algorithm reaches the preset maximum number of iterations, it terminates and outputs the Pareto solution set. In the NSGA-II algorithm crossover and mutation stage, based on the preset polynomial mutation distribution index and penalty sensitivity, the heat mismatch factor, local topological sensitivity and the preset penalty sensitivity of each individual in the current parent population are multiplied together and the result is negative. The exponential function value of this negative value with the natural constant as the base is calculated. This exponential function value is multiplied by the preset polynomial mutation distribution index to obtain the dynamically corrected variable length. Based on the dynamically corrected variable length, the mutation operation is performed on the current population individuals to generate the next generation of offspring population to be screened. During the elite retention phase of the NSGA-II algorithm, the generated next-generation offspring population to be screened is merged with the current parent population to form a hybrid population with double the size. In the hybrid population, two individuals adjacent to the current calculated population are found. The differences between these two adjacent individuals in each target dimension are calculated, normalized, and summed to obtain the geometric crowding degree of the current calculated population. Combining the preset guidance amplification coefficient, the preset guidance amplification coefficient is multiplied by the marginal benefit index of the current population individual and then incremented by one to obtain the weighting factor. The geometric crowding degree of the current calculated population individual is multiplied by the weighting factor to obtain the guidance crowding degree. When the NSGA-II algorithm constructs the next-generation offspring population, during the screening operation of the algorithm's next-generation population construction, individuals from the hybrid population are preferentially retained in descending order of guidance crowding degree to form a next-generation offspring population with half the size.

[0015] Furthermore, using the annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions corresponding to the Pareto solution set as negative evaluation indicators, the objective weights of the three evaluation indicators are calculated using the entropy weight-TOPSIS method. Combining the objective weights, the Euclidean distance between the Pareto solution set and the positive and negative ideal solutions is calculated using the approximation ideal ranking method. The comprehensive closeness of each solution in each Pareto solution set is obtained, and the solution with the largest comprehensive closeness is extracted. The joint optimization variable combination corresponding to this solution is output as the optimal joint optimization variable.

[0016] Compared with the prior art, the beneficial effects of the present invention are: This invention, based on a greenhouse spatial construction model, defines phase change temperature and film spacing as joint optimization variables, setting a constraint range under a preset target fermentation temperature and spatial interference limitations. The initial soil temperature is calculated by extracting meteorological and soil parameters, and this temperature, along with the target fermentation temperature, is set as the physical field boundary condition. Numerical simulations are used to obtain annual heating energy consumption, life-cycle cost, and life-cycle carbon emission data, thus constructing a multi-objective sample dataset and building a surrogate model. Compared to existing technologies that only consider construction costs and direct insulation effects, this design incorporates life-cycle dynamic costs and life-cycle carbon emissions into the joint optimization space, comprehensively quantifying the multi-dimensional integrated performance of the greenhouse under long-term operation. More importantly, a heat mismatch factor evaluation function is constructed based on the ratio of the Rayleigh number of air convection to the Stefan number of phase change heat storage; a local topological sensitivity evaluation function is constructed by extracting the first-order partial derivatives of the surrogate model to calculate the local maximum gradient; and a marginal benefit index evaluation function is constructed by calculating the marginal substitution rate between the target normalized differences. The establishment of these three functions provides precise thermodynamic perception and economic guidance for subsequent optimization, enabling the underlying logic to identify the degree of distortion of physical characteristics and to distinguish high-performance engineering solutions. This invention also employs a multi-objective optimization algorithm to iteratively optimize the surrogate model, calling the three evaluation functions mentioned above during iteration. Specifically, a thermal mismatch factor and local topological sensitivity are introduced to construct a penalized mutation operator to generate a population of offspring to be screened, overcoming the deficiency in existing evolutionary algorithms where crossover and mutation actions are only controlled by pure probability distributions. By intervening through the thermophysical mismatch mechanism between natural convection and phase change material heat storage in greenhouse spaces, and considering the local topological gradient of the objective function surface, the algorithm is effectively prevented from getting stuck in high gradient sensitive regions during iteration, reducing the generation of thermophysically distorted and invalid solutions, and saving computational power. Furthermore, a marginal benefit index is introduced to construct a directional congestion operator to screen and generate an offspring population, changing the existing method of relying solely on pure geometric congestion during truncation screening. By quantitatively considering the relative replacement rate between environmental reduction and cost for candidate solutions, the incorrect elimination of advantageous solutions with high marginal benefits in overcrowded regions is prevented, ensuring that the final Pareto solution set and the optimal joint optimization variables selected based on proximity have higher practical engineering guidance value. Attached Figure Description

[0017] Figure 1 This is a schematic diagram of the overall method flow of the present invention; Figure 2 This is a schematic cross-sectional view of the greenhouse and biogas fermentation tank of the present invention. Figure 3 This is a three-dimensional structural diagram of the greenhouse and biogas fermentation tank of the present invention; Figure 4 This is a schematic diagram illustrating the influence of the film spacing and phase change temperature on annual heating energy consumption in this invention. Figure 5 This is a schematic diagram illustrating the impact of coating spacing and phase change temperature on the total life cycle cost of this invention. Figure 6 This is a schematic diagram illustrating the impact of the coating spacing and phase change temperature on carbon emissions throughout the entire life cycle of this invention. Detailed Implementation

[0018] To make the objectives, technical solutions, and advantages of this invention clearer, the invention will be further described in detail below with reference to specific embodiments.

[0019] It should be noted that, unless otherwise defined, the technical or scientific terms used in this invention should have the ordinary meaning understood by one of ordinary skill in the art to which this invention pertains. The terms "first," "second," and similar terms used in this invention do not indicate any order, quantity, or importance, but are merely used to distinguish different components. Terms such as "comprising" or "including" mean that the element or object preceding the word encompasses the elements or objects listed following the word and their equivalents, without excluding other elements or objects. Terms such as "connected" or "linked" are not limited to physical or mechanical connections, but can include electrical connections, whether direct or indirect. Terms such as "upper," "lower," "left," and "right" are used only to indicate relative positional relationships; when the absolute position of the described object changes, the relative positional relationship may also change accordingly.

[0020] Example: Please see Figures 1 to 3 The present invention provides a technical solution: A method for optimizing design parameters of a biogas-use double-layer membrane phase change energy storage greenhouse, comprising the following steps: Step 1: Based on the greenhouse space construction model, the phase change temperature and the film covering distance are defined as joint optimization variables. Combined with the preset target fermentation temperature and spatial interference constraints, the constraint range of the joint optimization variables is set. In this embodiment, based on the preset target fermentation temperature, a temperature lower than the target fermentation temperature is set as the lower limit of the phase change temperature, and a temperature not lower than the target fermentation temperature is set as the upper limit of the phase change temperature, thus setting a constraint range for the phase change temperature; Based on the building structure, the effective static air layer thickness threshold is set as the lower limit, and the maximum critical distance to prevent physical contact between the internal membrane and the building's external envelope is set as the upper limit. Specifically, the membrane spacing is the distance between the internal membrane on the shady side of the biogas digester and the shady side insulation layer. Based on the architectural structure of a double-layer membrane-covered phase change energy storage greenhouse, a phase change plaster layer is applied to the inner surface of the greenhouse partition wall adjacent to the biogas digester, with the inner membrane covering the outer surface of the digester. According to this structure, the phase change temperature and the membrane spacing are defined as joint optimization variables. In practical engineering applications, the constraint range of these joint optimization variables needs to be flexibly set according to different fermentation process requirements. In a specific embodiment of this invention, when the biogas fermentation adopts a mesophilic fermentation process, i.e., the target fermentation temperature requirement is 35℃, to ensure that the phase change plaster layer can fully cross the phase change point for heat storage and release during the winter-spring transition season, this embodiment specifically sets the constraint range of the phase change temperature to 18℃ to 36℃. Simultaneously, considering the common spatial dimensions of double-layer membrane-covered greenhouses, to ensure maximum thermal resistance of the air layer without structural interference, this embodiment sets the specific constraint range of the membrane spacing to 0.2 meters to 1 meter. It should be understood that when the fermentation process is changed to high-temperature fermentation, such as 55℃, those skilled in the art can adjust the upper and lower limits of the phase change temperature according to the constraint rules described in this invention.

[0021] Step 2: Extract meteorological and soil parameters to calculate the initial soil temperature and set it together with the target fermentation temperature as the physical field boundary condition. Set multiple sets of joint optimization variables according to the constraint range and use numerical simulation under the physical field boundary condition to obtain annual heating energy consumption, life cycle cost and life cycle carbon emission data, forming a multi-objective sample dataset. In this embodiment, the average surface temperature, annual periodic surface temperature and annual temperature fluctuation period of the project site are extracted as meteorological parameters, and the soil thermal conductivity is extracted as soil physical property parameters. Based on the principle of unsteady periodic heat conduction, a semi-infinite flat wall periodic heat conduction analytical solution is constructed to calculate and obtain the initial soil temperature at the preset target depth. The initial soil temperature at the target depth and the target fermentation temperature are set together as the physical field boundary conditions. Multiple sets of joint optimization variables are substituted into the established physical field model. The target fermentation temperature is used as the internal heat source boundary, and the initial soil temperature and meteorological parameters are used as the external environment boundary. Fluid dynamics and heat transfer transient coupling simulation is performed using COMSOL Multiphysics to calculate the cumulative dynamic heating load required to maintain the target fermentation temperature. Based on the principles of heat conservation and electrothermal conversion efficiency, the ratio of the cumulative dynamic heating load to the conversion parameter is calculated using the preset electrothermal constant and the energy efficiency ratio of the heating equipment as conversion parameters to obtain the annual heating energy consumption.

[0022] When establishing the thermophysical field model of a biogas fermentation greenhouse, due to the deep thermodynamic coupling between the greenhouse structure and the foundation soil, it is necessary to extract meteorological parameters and soil physical property parameters of the project site to calculate the initial soil temperature at the target depth. This calculation method, based on the principle of unsteady periodic heat conduction to construct an analytical solution for the periodic heat conduction of a semi-infinite flat wall, can accurately characterize the temperature delay and attenuation characteristics of the soil at different depths with seasonal changes. This avoids the huge transient calculation errors caused by arbitrarily assuming the initial temperature in traditional numerical simulations, and provides the model with a dynamic equilibrium start-up condition that conforms to the laws of real natural climate. Temperature and the preset target fermentation temperature are set together as the physical field boundary conditions, which perfectly construct the real heat transfer driving force field between the heat source inside the greenhouse and the cold source in the external environment. Under this rigorous boundary condition, the cumulative dynamic heating load required to maintain the target fermentation temperature can be accurately reflected through transient coupling simulation of fluid mechanics and heat transfer. Furthermore, based on the principles of heat conservation and electrothermal conversion efficiency, this load is divided by a conversion parameter that includes the electrothermal constant and the equipment energy efficiency ratio to obtain the annual heating energy consumption. This process greatly ensures the authenticity of the underlying sample data on which subsequent fitting and optimization depend and the necessity of engineering decisions. In this calculation logic, the annual heating energy consumption is used as... As the core dependent variable, it essentially reflects the total external energy input required to withstand external cold loads and maintain a constant temperature throughout the fermentation process under a specific building configuration. Its magnitude directly constitutes the physical foundation for evaluating the long-term operational economics and carbon emission environmental costs of an engineering project. The evolution of this dependent variable profoundly depends on the combined effect of the two independent variables: phase change temperature and membrane spacing. The membrane spacing directly determines the thickness of the air layer between the two membranes and dominates the overall resistance to unsteady convective heat transfer within the membrane. The phase change temperature determines whether the phase change plaster layer can accurately cross the phase change point during the winter-spring weather transition to fully release latent heat. Flow buffering efficiency, from the perspective of mathematical correlation analysis between variables, when the phase change temperature gradually approaches the internal target fermentation temperature and the film spacing approaches the critical thickness that can just suppress convection, the thermophysical potential of the greenhouse enclosure structure is maximized. At this time, the dependent variable annual heating energy consumption shows a significant negative correlation and decreases as the independent variable approaches this ideal parameter range. Conversely, if the independent variable deviates from the optimal thermodynamic coupling range, resulting in violent natural convection heat dissipation between the double film or the phase change material becoming an ineffective sensible heat carrier due to temperature mismatch, the dependent variable annual heating energy consumption will show a sharp nonlinear positive correlation increase as the degree of deviation of the independent variable increases.

[0023] Based on the principle of the time value of money comprehensive discount model in the P1-P2 economic model, the present value factor of total expenditure within the economic analysis period and the ratio of total expenditure to initial investment within the economic analysis period are calculated. The specific calculation logic is as follows: Based on the preset economic analysis period, energy price growth rate, and market discount rate, when the energy price growth rate and the market discount rate are equal, the ratio of the economic analysis period to the energy price growth rate plus one is calculated as the present value factor of the total expenditure within the economic analysis period. When the energy price growth rate and the market discount rate are not equal, the ratio of the energy price growth rate plus one to the market discount rate plus one is calculated as the economic analysis period raised to the power of the economic analysis period. The result of subtracting the power from one is then divided by the difference between the market discount rate and the energy price growth rate to obtain the present value factor of the total expenditure within the economic analysis period. Based on the preset initial investment and down payment ratio, loan term, actual repayment period, loan interest rate, ratio of annual maintenance costs to initial investment, and ratio of resale price to initial investment, calculate the difference between the initial investment and down payment ratio and multiply it by the present value factor calculated using the market discount rate for the actual repayment period, and divide it by the present value factor calculated using the loan term and loan interest rate to obtain the loan cost present value ratio. Calculate the product of the ratio of annual maintenance costs to initial investment and the present value factor calculated using the market discount rate for the economic analysis period to obtain the maintenance cost present value ratio. Calculate the quotient of the ratio of resale price to initial investment raised to the power of the market discount rate plus one for the economic analysis period to obtain the residual value present value ratio. Add the initial investment and down payment ratio, the loan cost present value ratio, and the maintenance cost present value ratio, and subtract the residual value present value ratio to calculate the ratio of total expenditure to initial investment within the economic analysis period. In calculating heating operation costs, the annual heating energy consumption is multiplied by the preset electricity price to obtain the annual heating operation costs. In the initial investment cost calculation for the renovation, the product of the sum of the film covering construction cost and the film covering unit price and the film covering area is calculated; the product of the sum of the galvanized steel pipe construction cost and the galvanized steel pipe unit price and the galvanized steel pipe length is calculated; the product of the plastering construction cost and the plastering area is calculated; and the product of the phase change energy storage material unit price and the phase change energy storage material mass is calculated. The sum of the four product results in the initial investment cost calculation for the renovation is then obtained. Based on the above calculation results, in the calculation of the total life cycle cost, the annual heating operation cost is multiplied by the present value factor of the total expenditure within the economic analysis period, and the initial investment cost of the renovation is multiplied by the ratio of the total expenditure within the economic analysis period to the initial investment. The sum of the two products in the calculation of the total life cycle cost is then used to obtain the total life cycle cost.

[0024] After completing the thermophysical field simulation and obtaining the annual heating energy consumption, the P1-P2 economic model needs to be introduced to conduct full life-cycle cost accounting for the biogas greenhouse project. This calculation logic, based on the principle of comprehensive discounting of the time value of money, can accurately calculate complex dynamic realities such as energy price fluctuations, loan interest rate discounting, and equipment maintenance depreciation over a long period. This overcomes the shortcomings of traditional static investment returns, which cannot accurately measure the long-term operating capital burden, and provides a quantitative benchmark with practical guiding significance for engineering decisions. In this economic evaluation calculation, the full life-cycle cost, as the dependent variable, essentially reflects the total financial cost incurred during the entire service life of the greenhouse building, from construction to final demolition. Its technical effect lies in providing rigorous engineering cost constraints for multi-objective intelligent optimization, preventing the algorithm from evolving into a distorted parameter combination with optimal thermodynamic performance but extremely high cost and impracticality. The total life-cycle cost is directly controlled by two independent variables: phase change temperature and film spacing. The film spacing directly determines the amount of supporting steel pipes used and the area of ​​flexible film coverage, while the phase change temperature determines the selection and scale of use of specific phase change plastering materials. These bill of materials requirements derived directly from the independent variables constitute the basic economic parameter of initial investment cost for the renovation. At the same time, the change in the overall thermal resistance of the greenhouse driven by the independent variables determines the scale of annual heating energy consumption, which is then directly mapped to the long-term annual heating operation cost by multiplying it by the electricity price. From the analysis of the mathematical correlation between specific variables, as the film spacing increases or the phase change material related structures increase, the initial investment cost of the renovation shows a direct positive correlation increase. However, the strengthening of this heat storage boundary will significantly reduce the intrusion of external cold load, resulting in a strong negative correlation decrease in the annual heating energy consumption required to maintain the constant temperature of fermentation and the corresponding annual heating operation cost. The P1-P2 model precisely calculates the ratio of the present value factor of operating costs to the present value of investment, and integrates this positive correlation of short-term initial investment and negative correlation of long-term operating costs across time. Finally, it calculates the nonlinear numerical characteristics of the life cycle cost of the dependent variable, so that the change of the independent variable is essentially transformed into an engineering trade-off process of finding the lowest comprehensive cost between short-term construction costs and long-term operation energy saving.

[0025] The annual heating energy consumption is multiplied by the preset economic lifespan and the grid carbon emission factor to calculate the carbon emissions during the operation phase. Extract the bill of materials from the greenhouse space construction model. Based on the preset carbon emission factors of various building materials, multiply the carbon emission factors of various building materials by their corresponding bill of materials and sum them up to calculate the carbon emissions during the production stage. Based on the preset carbon emission factors and average transportation distances of various building materials, the carbon emission factors and average transportation distances of various building materials are multiplied together with their corresponding bill of materials and summed to calculate the carbon emissions during the transportation stage. The carbon emissions during the construction phase are estimated using a proportional method, with the carbon emissions during the construction phase calculated as 4% of the total carbon emissions during the production and transportation phases, and the carbon emissions during the demolition phase calculated as 10% of the total carbon emissions during the construction phase. The carbon emissions during the operation phase, production phase, transportation phase, construction phase, and demolition phase are summed to obtain the carbon emissions data for the entire life cycle. A multi-objective sample dataset is constructed based on annual heating energy consumption, total life-cycle costs, and total life-cycle carbon emissions.

[0026] After establishing economic evaluation benchmarks, further in-depth calculations of the life-cycle emission costs of biogas greenhouse projects are needed. This carbon emission accounting logic, based on the whole-process evaluation theory, can comprehensively capture the physical emission details of the entire chain from building material acquisition to disposal, thereby avoiding the emission transfer problem caused by the excessive accumulation of high-energy-consuming building materials in pursuit of energy conservation during operation. This provides an objective and rigorous quantitative scale for evaluating the green attributes of engineering technical solutions. In this accounting process, the life-cycle carbon emission data, as the dependent variable, essentially reflects the total greenhouse gas emission equivalent of the greenhouse building during its complete service life. Its technical effect is to set strict environmental limit constraints for multi-dimensional intelligent optimization, enabling the NSGA-II algorithm to fully consider long-term ecological costs when exploring multi-objective extreme values. The magnitude of the life-cycle carbon emission data is composed of the emission components from multiple stages such as operation, production, transportation, construction, and demolition, and its depth is controlled by the variation trajectory of the two independent variables: phase change temperature and membrane spacing. Since the physical extension of the cladding spacing directly determines the length of the supporting steel pipe and the usable area of ​​the flexible cladding, and the selection of the phase change temperature directly determines the specific ratio and laying quality of the phase change energy storage material, these changes in the bill of materials quantity directly driven by the phase change temperature and cladding spacing constitute the direct multiplier basis for calculating carbon emissions in the production and transportation stages by combining corresponding factors. This, in turn, leads to the estimated carbon emissions for the construction and demolition stages according to a fixed proportional coefficient. From a deeper mathematical correlation analysis, with a moderate expansion of the cladding spacing or an increase in the amount of phase change material used, the massive consumption of building materials leads to a significant positive correlation increase in carbon emissions during production, transportation, and subsequent construction and demolition stages based on proportional estimation. However, this enhanced building insulation and thermophysical energy storage properties significantly resist the intrusion of external cold loads, resulting in a significant decrease in the annual heating energy consumption required to maintain a constant internal fermentation temperature. Consequently, the carbon emissions during the operation stage, obtained by multiplying this energy consumption by the grid emission factor, show a strong negative correlation and a decreasing trend. The carbon accounting step in multi-objective optimization precisely integrates the positive correlation of emissions increase during the physical phase and the negative correlation of emissions decrease during the operational phase, thus finely characterizing the nonlinear numerical features of carbon emission data throughout the entire life cycle in the multi-dimensional optimization space. This transforms the changes in phase change temperature and membrane spacing into an engineering calculation process that seeks the optimal balance between high consumption of building materials in the short term and low emissions during long-term operation.

[0027] Step 3: Based on the multi-objective sample dataset, construct a proxy model for annual heating energy consumption, total life cycle cost, and total life cycle carbon emissions using polynomial surface fitting; In this embodiment, based on a multi-objective sample dataset, the combination of joint optimization variables in each group of the multi-objective sample dataset is used as the independent variable, and the corresponding annual heating energy consumption, life-cycle cost and life-cycle carbon emission data are used as the dependent variables. Based on the principle of multivariate higher-order nonlinear regression, a sixth-order bivariate polynomial surface fitting is performed to obtain the constant term, the higher-order main effect regression coefficients of each independent variable and the multivariate higher-order cross-coupling effect regression coefficients. Using the obtained regression coefficients, three sixth-order polynomial surrogate models are constructed, with phase change temperature and coating spacing as inputs and annual heating energy consumption, life-cycle cost and life-cycle carbon emission as outputs.

[0028] After acquiring a multi-objective sample dataset covering annual heating energy consumption, life-cycle costs, and life-cycle carbon emissions, the computational demands of real-world transient coupling simulations and life-cycle economic and environmental assessments are enormous. Directly embedding this costly assessment process into subsequent multi-objective optimization algorithms requiring tens of thousands of iterations would lead to a complete loss of control over the engineering computation timescale and even prevent convergence. Therefore, it is necessary to construct a corresponding mathematical proxy model based on the multi-objective sample dataset using polynomial surface fitting techniques. This transformation from low-level physical simulation to high-dimensional mathematical mapping reduces the time-consuming solution of multi-physics equations to explicit algebraic polynomial operations that can be completed in a very short time. This removes the computational burden from the algorithm's exploration of a massive optimization space, making high-frequency population propagation and large-scale global optimization feasible in practical engineering applications. In this fitting process, the annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions output by the three sixth-order polynomial surrogate models are established as independent dependent variables. Essentially, these variables reflect the greenhouse's energy input requirements in the thermophysical dimension, the comprehensive economic costs in the financial and time dimension, and the scale of greenhouse gas emissions in the full-process evaluation dimension, given any set of design parameters. The technical effect of this explicit expression is that it provides a continuously differentiable function basis for subsequent algorithms to construct mutation and screening operators, enabling real-time gradient calculation and marginal benefit assessment. The high-precision approximation of annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions on the fitting surface deeply depends on the interaction mapping between the two independent variables: phase change temperature and mulch spacing. The reason why phase change temperature and mulch spacing can determine the trends of the above three evaluation indicators is because there is a highly complex nonlinear thermodynamic coupling effect between the air layer thermal resistance evolution dominated by mulch spacing and the latent heat buffering efficiency controlled by phase change temperature. This underlying physical and economic coupling characteristic is precisely quantified and solidified in the surrogate model through the high-order main effect regression coefficients and cross-coupling effect regression coefficients in the polynomial surface fitting. To explain the connotation of this algebraic expression from the perspective of specific mathematical correlations, when the phase change temperature and the film-coating spacing change, for example, the film-coating spacing gradually increases within its effective threshold or the phase change temperature is optimized towards the range that matches the fermentation process, the improvement of physical insulation performance will trigger strong nonlinear feedback in the surrogate model. This will cause the annual heating energy consumption and the carbon emissions throughout the entire life cycle, which are dominated by the operation phase, to show an accelerated negative correlation decay surface. At the same time, affected by the increase in the physical consumption of building materials such as phase change materials and film-coated steel frames, the life cycle cost will show a complex wave-shaped or even saddle-shaped mapping relationship in different parameter ranges, which is formed by the game between the positive correlation of the initial investment increase and the negative correlation of the operating cost decrease. The sixth-order bivariate polynomial surrogate model, with its sufficiently flexible high-dimensional algebraic surface features, can capture and reconstruct the deep nonlinear mapping law of mixed positive and negative correlation between the phase change temperature and the film-coating spacing and the above three comprehensive evaluation indicators.

[0029] Step 4: Based on the population and physical field boundary conditions to be input, calculate the ratio of the natural logarithm of the Rayleigh number of air convection to the Stefan number of phase change heat storage to construct the thermal mismatch factor evaluation function, extract the first-order partial derivatives of the surrogate model and calculate the local maximum gradient to construct the local topological sensitivity evaluation function, and calculate the marginal substitution rate between the target normalized differences to construct the marginal benefit index evaluation function. In this embodiment, the entire combination of candidate joint optimization variables to be input is regarded as a population, and a single joint optimization variable in the combination of candidate joint optimization variables to be input is regarded as an individual in the population. Based on each individual population and physical field boundary conditions, according to the preset gravitational acceleration and the thermal expansion coefficient, kinematic viscosity and thermal diffusivity of the air between the double-layer membrane, and according to the internal target fermentation temperature and the external outdoor weather temperature, the absolute value of the difference between the two is calculated as the internal and external temperature difference. The gravitational acceleration, the thermal expansion coefficient of the air between the double-layer membrane, the internal and external temperature difference and the cube of the membrane spacing to be input are multiplied together and then divided by the product of the kinematic viscosity and thermal diffusivity of the air between the double-layer membrane to calculate the Rayleigh number characterizing the natural convection intensity of the air between the double-layer membrane. According to the preset specific heat capacity and latent heat enthalpy of the phase change material, the absolute value of the difference between the phase change temperature and the internal target fermentation temperature is calculated and multiplied by the specific heat capacity of the phase change material and then divided by the latent heat enthalpy to calculate the Stefan number characterizing the heat storage and release characteristics of the phase change material. The natural logarithm of the Rayleigh number is calculated and divided by the Stefan number to construct the thermal mismatch factor evaluation function. The constructed thermal mismatch factor evaluation function is as follows: in, Indicates the thermal mismatch factor. To represent the calculation of the natural logarithm, This represents the acceleration due to gravity, taken as 9.8 m / s². 2 , The coefficient of thermal expansion of the air between the two layers of film is obtained by referring to a standard table of air thermal properties. For an ideal gas, it is usually taken as the reciprocal of the absolute temperature. This represents the internal and external temperature difference, obtained by calculating the absolute value of the difference between the preset internal target fermentation temperature and the external outdoor ambient temperature. Indicates the intercellular distance within a population. The kinematic viscosity of the air between the two layers of coating is obtained by referring to a standard air thermal property table. The thermal diffusivity of the air between the two layers of coating is indicated by consulting a standard air thermal property table. This indicates the specific heat capacity of the phase change material, based on the selected phase change material provided by the material supplier. This indicates the phase transition temperature within an individual in a population. This indicates the preset internal target fermentation temperature. The latent heat enthalpy of the phase change material is provided by the material supplier. To provide rigorous guidance on the heat transfer physics mechanism for the NSGA-II algorithm optimization in step five, and to prevent the algorithm from blindly exploring purely mathematical probability distributions during the optimization process, it is necessary to first construct a heat mismatch factor evaluation function at this stage. The construction logic of this function is deeply rooted in the most critical thermodynamic contradiction in greenhouse engineering, namely the dynamic game between the natural convection heat loss of air between the double-layer covering and the latent heat buffering capacity of the phase change plaster layer. By introducing the Rayleigh number, which characterizes the intensity of fluid convection, and the Stefan number, which characterizes the heat storage and release characteristics of the phase change material, these two complex fluid dynamics and thermophysical properties are quantitatively integrated into a dimensionless physical evaluation benchmark. The rationale and necessity of constructing this calculation formula in advance lies in the fact that the NSGA-II algorithm, as a pure algebraic optimization tool, cannot perceive the thermodynamic distortions that engineering parameters may cause at the physical level. Without the real-time constraints of this underlying physical mechanism in future iterations, it is easy to generate invalid parameter solutions that lead to space flow control or complete failure of phase change energy storage. This evaluation function can act as a mathematical operator to help directly identify and avoid regions of thermophysical deterioration, thereby significantly improving the optimization effectiveness of computing resources in complex solution spaces. In this mathematical framework, the thermal mismatch factor, as a specific dependent variable, deeply reflects the relative imbalance between the structural heat dissipation risk and the material thermal storage efficiency of the greenhouse under specific parameter configurations. The technical effect of this physical quantification is that it can directly serve as a physical scale for dynamically correcting the variable asynchrony length of the penalty in the subsequent NSGA-II algorithm, enabling the evolution of individual populations to have a profound thermodynamic perception. The numerical evolution of this dependent variable is completely controlled by the variation trajectories of the two independent variables, namely the film-coating spacing and the phase change temperature. The reason why the independent variables can have a decisive influence on the internal trend of the dependent variable is that the physical scale of the film-coating spacing directly dominates the critical point of convection excitation of the boundary layer air, while the deviation between the phase change temperature and the target fermentation temperature directly defines whether the phase change material is in a state of efficient latent heat release or an inefficient state of sensible heat accumulation.From a precise mathematical correlation perspective, as the independent variable of the membrane spacing increases, the air between the double membranes gradually overcomes viscous resistance, causing a sharp increase in the Rayleigh number, which characterizes convection. Since it is in the numerator of the logarithmic operation, this results in a highly significant positive correlation increase in the dependent variable, the heat mismatch factor, physically representing a nonlinear increase in the risk of convective heat loss. At the same time, when the independent variable of the phase change temperature gradually deviates from the preset internal target fermentation temperature in the calculation, the increase in the absolute value of the difference will directly cause the Stradivarius number, which characterizes the heat storage characteristics, to increase accordingly. Since the Stradivarius number is strictly set in the denominator of the formula, this deviation in temperature distribution will cause the dependent variable, the heat mismatch factor, to show a clear inversely proportional negative correlation decrease. Utilizing this unique positive and negative correlation interaction mapping, the evaluation function successfully transforms the small physical deviation of the independent variable into a numerical feedback that is highly sensitive to the dependent variable, laying a solid mathematical foundation for the subsequent NSGA-II algorithm to tighten the parameter variation range in advance before the heat transfer state tends to collapse.

[0030] For three sixth-order polynomial surrogate models of annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions, the first-order partial derivatives of phase transition temperature and film spacing for each individual in the population are obtained. The two first-order partial derivatives corresponding to each surrogate model are squared, summed, and the square root is taken to obtain the spatial gradient magnitude of each surrogate model at the current variable position. The maximum value of the spatial gradient magnitudes of the three sixth-order polynomial surrogate models is selected to construct a local topological sensitivity evaluation function. The constructed local topology sensitivity evaluation function is as follows: in, This indicates local topological sensitivity, and `max` indicates the operation of taking the maximum value. Indicates the target sample identifier. This indicates the annual heating energy consumption. Indicates total life-cycle cost. Indicates carbon emissions over the entire life cycle. This represents the first-order partial derivative operator. This represents the sixth-order polynomial proxy model corresponding to the target sample identifier. Based on the proxy model of annual heating energy consumption, life-cycle cost and life-cycle carbon emissions, a large number of fragile solutions with excellent theoretical prediction performance are hidden in its polynomial surface. In the construction and long-term operation of greenhouses, even a small construction error in the film spacing or a small fluctuation in the actual phase change temperature can lead to a precipitous decline in its thermodynamic and economic performance. In order to avoid this engineering hazard, it is necessary to construct a local topological sensitivity evaluation function through rigorous mathematical differentiation to quantify the strength of anti-interference ability, so as to ensure that the final selected combination of engineering parameters has extremely high parameter tolerance and operational stability. In this evaluation function, local topological sensitivity, as the dependent variable to be solved, specifically reflects the maximum spatial gradient modulus that the three surrogate models—annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions—can generate in response to minor physical perturbations of phase change temperature and coating spacing at the current parameter coordinate point. Its engineering implication is that it represents the maximum risk of sudden performance changes and the steepness of parameter boundaries when a specific design scheme encounters fluctuations in underlying variables. This precise quantification provides a numerical benchmark for the optimization process to directly identify and effectively avoid high-risk sensitive areas. The numerical evolution of local topological sensitivity is directly controlled by the specific spatial stationary points of the two independent variables, phase change temperature and coating spacing, on the polynomial surface of the surrogate model. Since the numerical changes in phase change temperature and coating spacing and their interaction directly trigger transient slope changes in higher-order principal effects and cross-coupling effects on the surface, the trajectory of these two independent variables fundamentally dominates the sharp fluctuations in local topological sensitivity. Analyzing the mathematical correlation characteristics in the formula derivation, when the phase transition temperature and coating spacing approach the physical boundary, constraint limit, or critical inflection point where there is a drastic thermodynamic change, the sum of squares of the partial derivatives of the surrogate model in that coordinate will expand sharply, directly causing the local topological sensitivity to show an extremely strong positive correlation surge. This surge in value intuitively and sensitively exposes the extreme instability of the region. Conversely, when the phase transition temperature and coating spacing gradually slide into the wide, flat, and highly resistant optimal solution basin region, the partial derivatives of the phase transition temperature and coating spacing rapidly converge to zero, and the local topological sensitivity will show a significant in the same direction of decay. This evaluation function uses this precise positive linkage feedback between the spatial displacement of the phase transition temperature and coating spacing and the local topological sensitivity to transform the invisible topological steepness risk hidden in the prediction surface of the surrogate model into a numerical indicator that can be directly read.

[0031] The arithmetic mean of the predicted values ​​of all individuals in the population output by three sixth-order polynomial surrogate models is calculated respectively. The population mean is obtained in the dimensions of annual heating energy consumption, life cycle cost and life cycle carbon emission. The normalized difference of each individual in the population relative to the population mean in the dimensions of annual heating energy consumption and life cycle carbon emission is calculated. The absolute values ​​of the two are added to obtain the comprehensive environmental benefit component. The normalized difference of each individual in the population relative to the population mean in the dimension of life cycle cost is calculated. The absolute value of the difference is taken and a preset small constant is added to obtain the cost component. The comprehensive environmental benefit component is divided by the cost component to construct the marginal benefit index evaluation function. The constructed marginal benefit index evaluation function is as follows: in, Indicates the marginal benefit index. This represents the normalized difference between an individual in the population and the population mean in terms of annual heating energy consumption. This represents the normalized difference between an individual in a population and the population mean in terms of life-cycle carbon emissions. This represents the normalized difference between an individual in the population and the population mean in terms of life-cycle cost. This represents a small constant that prevents the denominator from being zero, and its value is 10. -6 .

[0032] In the non-dominated solution truncation and screening stage of multi-objective optimization, to overcome the defect that the traditional pure geometric distance congestion calculation method is prone to erroneously eliminating solutions with extremely high engineering value advantages, it is necessary to further construct a marginal benefit index evaluation function. The rationality and necessity of this calculation lies in the fact that when candidate solutions are densely congested in the multi-dimensional objective space, pure mathematical distance measurement cannot measure the true cost-effectiveness of different parameter schemes in terms of environmental reduction and economic cost. Therefore, it is necessary to introduce the concept of marginal substitution rate from economics to inject precise engineering value guidance into the algorithm screening. In this calculation framework, the marginal benefit index, as the dependent variable to be solved, essentially reflects the relative substitution rate between the sum of the normalized differences of annual heating energy consumption and life-cycle carbon emissions and the normalized difference of life-cycle costs for each individual in the population relative to the population mean. The technical effect of this precise quantification is that it directly constitutes the core weighting factor for the subsequent construction of the guiding congestion operator, ensuring that the optimization iteration can accurately retain those high-quality parameter combinations that exchange huge environmental benefits for minimal economic costs in the dense local space of the population. The numerical evolution of the marginal benefit index, the dependent variable, is ultimately controlled by the two most basic independent variables: phase change temperature and film spacing. The reason why the independent variables can profoundly affect the various components within the dependent variable is that the setting of the phase change temperature and the physical extension of the film spacing directly determine the prediction output of each proxy model, thereby affecting the degree of deviation of each engineering target dimension from the population mean. Analyzing the ingenious correlation mapping logic in the formula, when the interlayer spacing is appropriately increased or the phase change temperature approaches the optimal thermodynamic coupling point, the significant improvement in physical insulation and heat storage performance will cause the annual heating energy consumption and life-cycle carbon emissions to be significantly lower than the population average. This performance directly leads to an extremely strong positive correlation surge in the comprehensive environmental benefits component, which is the numerator of the formula. However, at the same time, the increase in greenhouse building materials will inevitably lead to an increase in life-cycle costs, which will be higher than the population average. This causes the cost component, which is the denominator of the formula, to also show a positive correlation increase. The marginal benefit index evaluation function divides this synchronously increasing environmental numerator with the cost denominator, transforming the spatial displacement of phase change temperature and interlayer spacing into a nonlinear ratio feedback that measures how many units of environmental benefits can be brought about by each additional unit of economic cost. This gives the computational model a powerful computational insight to directly capture and retain the engineering solution with the highest cost-effectiveness on the complex and crowded Pareto front.

[0033] Step 5: Use a multi-objective optimization algorithm to iterate and optimize the surrogate model. Randomly generate an initial population within the constraints of the joint optimization variables. During iteration, call the three evaluation functions, introduce the heat mismatch factor and local topological sensitivity to construct a penalty mutation operator to generate a subpopulation to be screened, and then introduce the marginal benefit index to construct a directional crowding operator to screen and generate a subpopulation. Iterate and update until the termination condition is met, and output the Pareto solution set. In this embodiment, a preset number of joint optimization variable combinations are randomly generated within the constraints of the joint optimization variables as an initial population. When using the NSGA-II algorithm for iterative solution, the generated initial population is used as the parent population to call the heat mismatch factor evaluation function, the local topological sensitivity evaluation function, and the marginal benefit index evaluation function. When the next generation of offspring population is generated, the generated next generation of offspring population is used as the new parent population to continue calling the three evaluation functions for optimization iteration. When the NSGA-II algorithm reaches the preset maximum number of iterations, it terminates and outputs the Pareto solution set. In the NSGA-II algorithm crossover and mutation stage, based on the preset polynomial mutation distribution index and penalty sensitivity, the heat mismatch factor, local topological sensitivity and the preset penalty sensitivity of each individual in the current parent population are multiplied together and the result is negative. The exponential function value of this negative value with the natural constant as the base is calculated. This exponential function value is multiplied by the preset polynomial mutation distribution index to obtain the dynamically corrected variable length. Based on the dynamically corrected variable length, the mutation operation is performed on the current population individuals to generate the next generation of offspring population to be screened. In the crossover and mutation phase of the NSGA-II algorithm, to overcome the limitations of conventional evolutionary algorithms that are only controlled by pure probability distributions and lack guidance from underlying physical mechanisms, a rigorous mathematical penalty mechanism is needed to calculate the dynamically corrected variable length. The dynamically corrected variable length, as the core dependent variable of this mutation operation, essentially reflects the physical search range that individuals in the population are allowed to explore when reproducing the next generation in the current multidimensional parameter space. This technique, which imbues the optimization process with spatial constraints, precisely restrains candidate solutions that are in adverse thermodynamic states or on the edge of steep terrain, preventing them from further generating severely distorted and invalid thermophysical solutions that deviate significantly from practical engineering principles during blind probability mutation, thereby greatly saving valuable computing resources. The precise evolution of the dynamically corrected variable asynchrony length is deeply controlled by two key independent variables: the thermal mismatch factor and the local topological sensitivity. The reason why the thermal mismatch factor and the local topological sensitivity can jointly dominate the trend of the dynamically corrected variable asynchrony length is that the former represents the degree of physical imbalance between natural convection heat loss and latent heat buffering in the greenhouse, while the latter reveals the fragile limit of the anti-interference ability on the surrogate model prediction surface. These two risk parameters that characterize engineering disadvantages directly constitute the natural physical and geometric basis for imposing variation penalties. From a rigorous mathematical correlation and the underlying operational logic of the penalty operator, when the current parent population individual participating in the iteration exhibits an extremely high risk of thermophysical mismatch or its coordinates reside in a sensitive region where the objective function undergoes a drastic change, the values ​​of the thermal mismatch factor and local topological sensitivity corresponding to that individual will be at a high level. Under the set mathematical constraints, after multiplying these two continuously increasing independent variables with the preset penalty sensitivity and taking the negative value, the calculated exponential function value with the natural constant as the base will show a strong exponential decay. This directly leads to the extremely significant negative correlation and contraction trend of the dynamically corrected variable length obtained after multiplying the severely decayed exponential function value with the preset polynomial variation distribution exponent. This precise mapping design based on the natural exponential decay law successfully transforms the engineering risks represented by the thermal mismatch factor and local topological sensitivity into extremely severe numerical suppression of the dynamically corrected variable asynchrony length. This allows the NSGA-II algorithm to spontaneously impose physical exploration constraints on high-risk individuals when generating the next generation of offspring to be screened, ensuring that the entire parameter evolution process always steadily advances along a safe optimization space with high anti-interference capability and in accordance with the real heat transfer law.

[0034] During the elite retention phase of the NSGA-II algorithm, the generated next-generation offspring population to be screened is merged with the current parent population to form a mixed population with double the size. In the mixed population, two individuals adjacent to the current calculated population are found. The differences between these two individuals in each target dimension are calculated, normalized, and summed to obtain the geometric crowding degree of the current calculated population. Combined with a preset guidance amplification coefficient, the preset guidance amplification coefficient is multiplied by the marginal benefit index of the current population and then incremented by one to obtain the weighting factor. The geometric crowding degree of the current calculated population is multiplied by the weighting factor to obtain the guidance crowding degree. When the NSGA-II algorithm constructs the next-generation offspring population, during the screening operation of the algorithm constructing the next-generation offspring population, individuals from the mixed population are retained first in descending order of guidance crowding degree to form a next-generation offspring population with half the size. In the elite retention phase of the NSGA-II algorithm, when the size of the merged mixed population exceeds the next generation quota and truncation is required, traditional optimization logic often relies solely on pure multidimensional spatial distance to maintain population diversity. This purely geometric measurement method easily leads to the algorithm incorrectly eliminating solutions with extremely high engineering exchange cost-effectiveness in areas of excessively dense solution space. To overcome this fundamental flaw, guided crowding is calculated to inject a real engineering value judgment benchmark into the survival of the fittest in the population. In this precise selection calculation, guided crowding is the core dependent variable being solved. Its essence reflects the comprehensive survival weight of the geometrical distribution sparsity of individuals in the multi-objective parameter space and their marginal substitution rate between environmental reduction and cost. The technical effect of constructing this dependent variable is that it completely reshapes the critical ranking criterion of non-dominated solutions, making the final selected next generation population not only spatially uniform but also possessing extremely high decision-making value in practical engineering applications. The numerical evolution of orientation congestion is deeply controlled by two key independent variables: geometric congestion and marginal benefit index. Geometric congestion and the marginal benefit index jointly dominate the final value of orientation congestion because the former accurately depicts the spatial congestion state at a purely mathematical level by calculating the normalized difference between adjacent individuals across various evaluation dimensions, while the latter directly quantifies the real environmental benefits gained at a specific coordinate point by sacrificing unit economic costs. This deep integration of the two precisely compensates for the logical loophole that a single geometric scale cannot identify the merits and demerits of physical and economic attributes. Analyzing from a rigorous mathematical correlation mechanism, under the pre-set algebraic mapping rules, orientation congestion is obtained by directly multiplying geometric congestion by a weighting factor constructed based on the marginal benefit index. When a candidate individual is located in an extremely sparse boundary or isolated region in the solution space, its larger geometric congestion directly causes a significant positive correlation increase in orientation congestion, thus naturally ensuring the diversity of global exploration by the population. More importantly, when a candidate has an extremely high marginal benefit index, multiplying this index by the guidance amplification coefficient and adding one will generate a weighting factor with a very large value. With the mathematical amplification structure of the multiplication, this expanded weighting factor will have an extremely strong positive pull on the guidance congestion degree. Relying on this positive correlation value surge driven by high marginal benefits, even if the high-quality engineering solution is unfortunately located in an extremely congested geometric coordinate region, its final calculated guidance congestion degree can still achieve a strong numerical reversal, thereby avoiding the indiscriminate elimination of conventional truncation operations and ensuring that parameter combinations that can exchange small investments for large energy-saving and emission-reduction benefits can be stably inherited to the next generation.

[0035] Step 6: Calculate the proximity of each solution in the Pareto solution set, and output the joint optimization variable combination with the highest proximity as the optimal joint optimization variable; In this embodiment, the annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions corresponding to the Pareto solution set are used as negative evaluation indicators. The entropy weight-TOPSIS method is used to calculate the objective weights of the three evaluation indicators. Combined with the objective weights, the Euclidean distance between the Pareto solution set and the positive and negative ideal solutions is calculated by the approximation ideal ranking method. The comprehensive closeness of each solution in each Pareto solution set is obtained. The solution with the largest comprehensive closeness is extracted, and the joint optimization variable combination corresponding to the solution is output as the optimal joint optimization variable.

[0036] In this method, the standard entropy weight-TOPSIS approach is used. Specifically, annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions are all used as negative evaluation indicators to construct an initial evaluation matrix. First, the range standardization of each negative indicator in the initial evaluation matrix is ​​performed. Based on the standardized data, the objective weight of each negative evaluation indicator is calculated using information entropy theory, and a weighted decision matrix is ​​constructed accordingly. The Euclidean distance between each candidate solution and the positive and negative ideal solutions in the weighted decision matrix is ​​then calculated.

[0037] This design parameter optimization selected five coating spacings: 0.2m, 0.4m, 0.6m, 0.8m, and 1.0m, and ten phase change temperatures: 18℃, 20℃, 22℃, 24℃, 26℃, 28℃, 30℃, 32℃, 34℃, and 36℃. A total of 50 sets of joint optimization variables were used, with a comparison of single-layer coating without phase change and a coating spacing of 0.4m without phase change. Annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions for a total of 52 sets of joint optimization variables were calculated, as detailed in Table 1. Table 1: Reference Table of Related Impacts of Joint Optimization Variables A surrogate model of annual heating energy consumption as a function of phase transition temperature and coating spacing was obtained by performing sixth-order polynomial surface fitting using MATLAB. The fitting effect is good. The specific sixth-order polynomial surrogate model for annual heating energy consumption is as follows: A surrogate model for life-cycle cost as a function of phase transition temperature and coating spacing was obtained by fitting a sixth-order polynomial surface using MATLAB. The fitting effect is good. The specific sixth-order polynomial surrogate model for the total life cycle cost is as follows: A surrogate model for life-cycle carbon emissions as a function of phase transition temperature and coating spacing was obtained by fitting a sixth-order polynomial surface using MATLAB. The fitting effect is good. The specific life-cycle carbon emission sixth-order polynomial surrogate model is as follows: The above formulas are derived from 50 calculation results (excluding the two control schemes) that have been processed and then fitted into MATLAB using a sixth-order polynomial surface. These models represent the annual heating energy consumption, total life-cycle cost, and total life-cycle carbon emissions as a function of phase transition temperature and membrane spacing. The influence curves are shown in the figures below. Figure 4 , Figure 5 and Figure 6 .

[0038] The above formulas are all dimensionless calculations. The formulas are derived from software simulations based on a large amount of collected data to obtain the most recent real-world results. The preset parameters in the formulas are set by those skilled in the art according to the actual situation.

[0039] The above embodiments can be implemented, in whole or in part, by software, hardware, firmware, or any other combination thereof. When implemented in software, the above embodiments can be implemented, in whole or in part, as a computer program product. Those skilled in the art will recognize that the units and algorithm steps of the various examples described in conjunction with the embodiments disclosed herein can be implemented by electronic hardware, or a combination of computer software and electronic hardware. Whether these functions are implemented in hardware or software depends on the specific application and design constraints of the technical solution.

[0040] The units described as separate components may or may not be physically separate. The components shown as units may or may not be physical units; they may be located in one place or distributed across multiple network units. Some or all of the units can be selected to achieve the purpose of this embodiment, depending on actual needs.

[0041] The above description is merely a specific embodiment of this application, but the scope of protection of this application is not limited thereto. Any changes or substitutions that can be easily conceived by those skilled in the art within the scope of the technology disclosed in this application should be included within the scope of protection of this application.

Claims

1. A method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas, characterized in that, The specific steps include: Step 1: Based on the greenhouse space construction model, the phase change temperature and the film covering distance are defined as joint optimization variables. Combined with the preset target fermentation temperature and spatial interference constraints, the constraint range of the joint optimization variables is set. Step 2: Extract meteorological and soil parameters to calculate the initial soil temperature and set it together with the target fermentation temperature as the physical field boundary condition. Set multiple sets of joint optimization variables according to the constraint range and use numerical simulation under the physical field boundary condition to obtain annual heating energy consumption, life cycle cost and life cycle carbon emission data, forming a multi-objective sample dataset. Step 3: Based on the multi-objective sample dataset, construct a proxy model for annual heating energy consumption, total life cycle cost, and total life cycle carbon emissions using polynomial surface fitting; Step 4: Based on the population and physical field boundary conditions to be input, calculate the ratio of the natural logarithm of the Rayleigh number of air convection to the Stefan number of phase change heat storage to construct the thermal mismatch factor evaluation function, extract the first-order partial derivatives of the surrogate model and calculate the local maximum gradient to construct the local topological sensitivity evaluation function, and calculate the marginal substitution rate between the target normalized differences to construct the marginal benefit index evaluation function. Step 5: Use a multi-objective optimization algorithm to iterate and optimize the surrogate model. Randomly generate an initial population within the constraints of the joint optimization variables. During iteration, call the three evaluation functions, introduce the heat mismatch factor and local topological sensitivity to construct a penalty mutation operator to generate a subpopulation to be screened, and then introduce the marginal benefit index to construct a directional crowding operator to screen and generate a subpopulation. Iterate and update until the termination condition is met, and output the Pareto solution set. Step 6: Calculate the proximity of each solution in the Pareto solution set, and output the joint optimization variable combination with the highest proximity as the optimal joint optimization variable.

2. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 1, characterized in that: Based on the preset target fermentation temperature, a temperature lower than the target fermentation temperature is set as the lower limit of the phase transition temperature, and a temperature not lower than the target fermentation temperature is set as the upper limit of the phase transition temperature, thus setting the constraint range of the phase transition temperature. Based on the building structure, the effective static air layer thickness threshold is set as the lower limit, and the maximum critical distance to prevent physical contact between the internal membrane and the building's external envelope is set as the upper limit. Specifically, the membrane spacing is the distance between the internal membrane on the shady side of the biogas digester and the insulation layer on the shady side.

3. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 1, characterized in that: The average surface temperature, annual periodic surface temperature, and annual temperature fluctuation period of the project site are extracted as meteorological parameters, and the soil thermal conductivity is extracted as a soil physical property parameter. Based on the principle of unsteady periodic heat conduction, a semi-infinite flat wall periodic heat conduction analytical solution is constructed to calculate the initial soil temperature at the preset target depth. The initial soil temperature at the target depth and the target fermentation temperature are set together as the physical field boundary conditions. Multiple sets of joint optimization variables are substituted into the established physical field model. The target fermentation temperature is used as the internal heat source boundary, and the initial soil temperature and meteorological parameters are used as the external environment boundary. Fluid dynamics and heat transfer transient coupling simulation is performed using COMSOL Multiphysics to calculate the cumulative dynamic heating load required to maintain the target fermentation temperature. Based on the principles of heat conservation and electrothermal conversion efficiency, the ratio of the cumulative dynamic heating load to the conversion parameter is calculated using the preset electrothermal constant and the energy efficiency ratio of the heating equipment as conversion parameters to obtain the annual heating energy consumption.

4. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 3, characterized in that: Based on the principle of the time value of money comprehensive discount model in the P1-P2 economic model, the present value factor of total expenditure within the economic analysis period and the ratio of total expenditure to initial investment within the economic analysis period are calculated. The specific calculation logic is as follows: Based on the preset economic analysis period, energy price growth rate, and market discount rate, when the energy price growth rate and the market discount rate are equal, the ratio of the economic analysis period to the energy price growth rate plus one is calculated as the present value factor of the total expenditure within the economic analysis period. When the energy price growth rate and the market discount rate are not equal, the ratio of the energy price growth rate plus one to the market discount rate plus one is calculated as the economic analysis period raised to the power of the economic analysis period. The result of subtracting the power from one is then divided by the difference between the market discount rate and the energy price growth rate to obtain the present value factor of the total expenditure within the economic analysis period. Based on the preset initial investment and down payment ratio, loan term, actual repayment period, loan interest rate, ratio of annual maintenance costs to initial investment, and ratio of resale price to initial investment, calculate the difference between the initial investment and down payment ratio and multiply it by the present value factor calculated using the market discount rate for the actual repayment period, and divide it by the present value factor calculated using the loan term and loan interest rate to obtain the loan cost present value ratio. Calculate the product of the ratio of annual maintenance costs to initial investment and the present value factor calculated using the market discount rate for the economic analysis period to obtain the maintenance cost present value ratio. Calculate the quotient of the ratio of resale price to initial investment raised to the power of the market discount rate plus one for the economic analysis period to obtain the residual value present value ratio. Add the initial investment and down payment ratio, the loan cost present value ratio, and the maintenance cost present value ratio, and subtract the residual value present value ratio to calculate the ratio of total expenditure to initial investment within the economic analysis period. In calculating heating operation costs, the annual heating energy consumption is multiplied by the preset electricity price to obtain the annual heating operation costs. In the initial investment cost calculation for the renovation, the product of the sum of the film covering construction cost and the film covering unit price and the film covering area is calculated; the product of the sum of the galvanized steel pipe construction cost and the galvanized steel pipe unit price and the galvanized steel pipe length is calculated; the product of the plastering construction cost and the plastering area is calculated; and the product of the phase change energy storage material unit price and the phase change energy storage material mass is calculated. The sum of the four product results in the initial investment cost calculation for the renovation is then obtained. Based on the above calculation results, in the calculation of the total life cycle cost, the annual heating operation cost is multiplied by the present value factor of the total expenditure within the economic analysis period, and the initial investment cost of the renovation is multiplied by the ratio of the total expenditure within the economic analysis period to the initial investment. The sum of the two products in the calculation of the total life cycle cost is then used to obtain the total life cycle cost.

5. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 4, characterized in that: The annual heating energy consumption is multiplied by the preset economic lifespan and the grid carbon emission factor to calculate the carbon emissions during the operation phase. Extract the bill of materials from the greenhouse space construction model. Based on the preset carbon emission factors of various building materials, multiply the carbon emission factors of various building materials by their corresponding bill of materials and sum them up to calculate the carbon emissions during the production stage. Based on the preset carbon emission factors and average transportation distances of various building materials, the carbon emission factors and average transportation distances of various building materials are multiplied together with their corresponding bill of materials and summed to calculate the carbon emissions during the transportation stage. The carbon emissions during the construction phase are estimated using a proportional method, with the carbon emissions during the construction phase calculated as 4% of the total carbon emissions during the production and transportation phases, and the carbon emissions during the demolition phase calculated as 10% of the total carbon emissions during the construction phase. The carbon emissions during the operation phase, production phase, transportation phase, construction phase, and demolition phase are summed to obtain the carbon emissions data for the entire life cycle. A multi-objective sample dataset is constructed based on annual heating energy consumption, total life-cycle costs, and total life-cycle carbon emissions.

6. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 1, characterized in that: Based on a multi-objective sample dataset, the combined optimization variables of each group in the multi-objective sample dataset are used as independent variables, and the corresponding annual heating energy consumption, life-cycle cost, and life-cycle carbon emission data are used as dependent variables. Based on the principle of multivariate higher-order nonlinear regression, a sixth-order bivariate polynomial surface fitting is performed to obtain the constant term, the higher-order main effect regression coefficients of each independent variable, and the multivariate higher-order cross-coupling effect regression coefficients. Using the obtained regression coefficients, three sixth-order polynomial surrogate models are constructed, with phase change temperature and coating spacing as inputs and annual heating energy consumption, life-cycle cost, and life-cycle carbon emission as outputs.

7. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 1, characterized in that: The entire combination of candidate joint optimization variables to be input is treated as a population, and a single joint optimization variable in the combination of candidate joint optimization variables to be input is treated as an individual in the population. Based on each individual population and physical field boundary conditions, according to the preset gravitational acceleration and the thermal expansion coefficient, kinematic viscosity and thermal diffusivity of the air between the double-layer membrane, and according to the internal target fermentation temperature and the external outdoor weather temperature, the absolute value of the difference between the two is calculated as the internal and external temperature difference. The gravitational acceleration, the thermal expansion coefficient of the air between the double-layer membrane, the internal and external temperature difference and the cube of the membrane spacing to be input are multiplied together and then divided by the product of the kinematic viscosity and thermal diffusivity of the air between the double-layer membrane to calculate the Rayleigh number characterizing the natural convection intensity of the air between the double-layer membrane. According to the preset specific heat capacity and latent heat enthalpy of the phase change material, the absolute value of the difference between the phase change temperature and the internal target fermentation temperature is calculated and multiplied by the specific heat capacity of the phase change material and then divided by the latent heat enthalpy to calculate the Stefan number characterizing the heat storage and release characteristics of the phase change material. The natural logarithm of the Rayleigh number is calculated and divided by the Stefan number to construct the thermal mismatch factor evaluation function. For three sixth-order polynomial surrogate models of annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions, the first-order partial derivatives of phase transition temperature and film spacing for each individual in the population are obtained. The two first-order partial derivatives corresponding to each surrogate model are squared, summed, and the square root is taken to obtain the spatial gradient magnitude of each surrogate model at the current variable position. The maximum value of the spatial gradient magnitudes of the three sixth-order polynomial surrogate models is selected to construct a local topological sensitivity evaluation function. The arithmetic mean of the predicted values ​​output by three sixth-order polynomial surrogate models for all individuals in the population is calculated. The population mean is obtained for the dimensions of annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions. The normalized difference of each individual relative to the population mean in the dimensions of annual heating energy consumption and life-cycle carbon emissions is calculated. The absolute values ​​of the two are added to obtain the comprehensive environmental benefit component. The normalized difference of each individual relative to the population mean in the dimension of life-cycle cost is calculated. The absolute value of the difference is taken and a preset small constant is added to obtain the cost component. The comprehensive environmental benefit component is divided by the cost component to construct the marginal benefit index evaluation function.

8. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 7, characterized in that: Within the constraints of the joint optimization variables, a preset number of joint optimization variable combinations are randomly generated as the initial population. When using the NSGA-II algorithm for iterative solution, the generated initial population is used as the parent population to call the thermal mismatch factor evaluation function, the local topological sensitivity evaluation function, and the marginal benefit index evaluation function. When the next generation of offspring population is generated, the generated next generation of offspring population is used as the new parent population to continue calling the three evaluation functions for optimization iteration. When the NSGA-II algorithm reaches the preset maximum number of iterations, it terminates and outputs the Pareto solution set. In the NSGA-II algorithm crossover and mutation stage, based on the preset polynomial mutation distribution index and penalty sensitivity, the heat mismatch factor, local topological sensitivity and the preset penalty sensitivity of each individual in the current parent population are multiplied together and the result is negative. The exponential function value of this negative value with the natural constant as the base is calculated. This exponential function value is multiplied by the preset polynomial mutation distribution index to obtain the dynamically corrected variable length. Based on the dynamically corrected variable length, the mutation operation is performed on the current population individuals to generate the next generation of offspring population to be screened. During the elite retention phase of the NSGA-II algorithm, the generated next-generation offspring population to be screened is merged with the current parent population to form a hybrid population with double the size. In the hybrid population, two individuals adjacent to the current calculated population are found. The differences between these two adjacent individuals in each target dimension are calculated, normalized, and summed to obtain the geometric crowding degree of the current calculated population. Combining the preset guidance amplification coefficient, the preset guidance amplification coefficient is multiplied by the marginal benefit index of the current population individual and then incremented by one to obtain the weighting factor. The geometric crowding degree of the current calculated population individual is multiplied by the weighting factor to obtain the guidance crowding degree. When the NSGA-II algorithm constructs the next-generation offspring population, during the screening operation of the algorithm's next-generation population construction, individuals from the hybrid population are preferentially retained in descending order of guidance crowding degree to form a next-generation offspring population with half the size.

9. The method for optimizing design parameters of a double-layer membrane phase change energy storage solar greenhouse for biogas as described in claim 8, characterized in that: Using the annual heating energy consumption, life-cycle cost, and life-cycle carbon emissions corresponding to the Pareto solution set as negative evaluation indicators, the objective weights of the three evaluation indicators are calculated using the entropy weight-TOPSIS method. Combining the objective weights, the Euclidean distance between the Pareto solution set and the positive and negative ideal solutions is calculated using the approximation ideal ranking method. The comprehensive closeness of each solution in each Pareto solution set is obtained, and the solution with the largest comprehensive closeness is extracted. The joint optimization variable combination corresponding to this solution is output as the optimal joint optimization variable.

Images