A method for predicting xenon oscillations in nuclear reactors based on data assimilation and physical constraints

Abstract

This invention discloses a method for predicting xenon oscillations in nuclear reactors based on data assimilation and physical constraints, relating to the field of xenon oscillation prediction in nuclear reactors. It addresses the difficulties in xenon oscillation prediction in existing technologies, the limitation of physical model methods to specific reactor types, and the high computational cost of data-driven methods. This invention establishes a xenon oscillation surrogate model with physical constraints based on data assimilation. The surrogate model consists of two parts: a parameter search part and a xenon oscillation process curve prediction part. The parameter search part refers to the model searching for globally optimal parameters for multiple operating conditions under physical constraints through data assimilation methods. The xenon oscillation process curve prediction part refers to using the fitting relationship between parameters to predict the xenon oscillation process under unknown operating conditions. This invention is also applicable to the field of surrogate model technology for predicting xenon oscillation curves.

Description

Technical Field

[0001] This invention relates to the field of nuclear reactor xenon oscillation prediction technology, and specifically to a method for predicting nuclear reactor xenon oscillations based on data assimilation and physical constraints. Background Technology

[0002] Xenon oscillations generally occur in large thermal neutron reactors with high neutron flux rates (such as large pressurized water reactors). Radial xenon oscillations converge, and their stability increases with increasing core burnup; however, the stability of axial xenon oscillations decreases with increasing core burnup or core power. Therefore, axial xenon oscillations are a worthy research and discussion topic in large thermal neutron reactors. In practical pressurized water reactors, xenon oscillations can cause axial power deviation Δ. I Periodic oscillation. When the power deviation Δ I When the power deviation exceeds ±3% FP (full power), it will trigger a main control warning; when the power deviation exceeds ±5% FP, the cumulative timing within 12 consecutive hours must not exceed 1 hour. Furthermore, xenon oscillations can cause shifts in the reactor heat pipe position and alter the power density peak factor, leading to increased power levels and temperatures in localized areas. Without control, this can even cause fuel element meltdown. Xenon oscillations also cause alternating temperature field changes in the reactor core, exacerbating temperature stress variations in the core materials and causing premature material failure. Therefore, to ensure the safe and stable operation of the unit, xenon oscillations must be detected early and addressed with appropriate intervention measures.

[0003] Xenon oscillations involve complex nonlinear processes, including the interaction between neutron flux, iodine-135, and xenon-135, which are difficult to describe using simple linear models. During reactor operation, the concentration of xenon-135 cannot be directly measured; its distribution can only be inferred indirectly (e.g., by measuring neutron flux), increasing prediction uncertainty. Therefore, effective prediction methods are needed. Accurate prediction of xenon oscillations allows for proactive measures to suppress them and avoid potential safety risks. In short, accurate prediction and effective control of xenon oscillations are crucial for ensuring reactor safety and optimized operation.

[0004] Among them, the physical model method proposed in the prior art is based on two-group one-dimensional diffusion equations and the concentration change equations of iodine and xenon. It uses Fourier sinusoidal series expansion to determine the spatial distribution of xenon, iodine, and flux, and introduces an axial difference parameter to transform the equations into non-homogeneous equations. This method proposes an analytical solution for the time-varying axial power shape index (ASI) of the reactor during xenon oscillation. Given the cross-sectional data under the equilibrium xenon conditions and the corresponding physical constants, the initial reactor axial power shape index (ASI) and its corresponding derivative can be input to predict the xenon oscillation process.

[0005] However, the physical model method has drawbacks: it requires a sufficient amount of pre-measured data to fit a cosine function during balanced xenon initialization, and it does not allow changes to the core conditions to make these data valid. This method has poor applicability to various operating conditions and is generally only applicable to specific reactor types.

[0006] The proposed data-driven approach in the prior art combines a data assimilation method with a dynamic mode decomposition (DMD) hybrid dynamic prediction scheme to predict the whole-core power distribution caused by xenon oscillations in the HPR1000 reactor. This coupled scheme is divided into two stages: (1) Reconstruction stage: In the initial stage, only the IDB-DA method is used for physical field reconstruction. (2) Prediction and reconstruction stage: In the subsequent stage, the detector node power values ​​predicted by DMD and the IDB-DA method are combined to predict the whole-core power distribution.

[0007] However, the disadvantages of the data-driven method are: (1) The implementation of the data assimilation method is relatively complex, requiring the selection of appropriate background physics fields and hyperparameters such as regularization factors, and a large amount of experimental data and adjustments. When the number of detectors configured in the reactor core is limited, insufficient data will lead to an increase in prediction error. This method is a data-driven method, which is not subject to physical constraints, and the reliability of the prediction results is relatively poor. (2) Although the DMD method has high computational efficiency, the IDB-DA method may require high computational costs when dealing with large-scale physics fields.

[0008] (3) The performance of data assimilation is highly dependent on the quality of the observation data. If the observation data contains large noise or errors, it may affect the accuracy of the prediction results. Summary of the Invention

[0009] This invention addresses the problems of difficulty in predicting xenon oscillations in existing technologies, the limitation of physical model methods to specific reactor types, and the high computational cost of data-driven methods.

[0010] To solve the above-mentioned technical problems, the present invention is achieved through the following technical solution:

[0011] This invention proposes a method for predicting xenon oscillations in nuclear reactors based on data assimilation and physical constraints. The method includes the following steps:

[0012] S1. Select five groups of xenon transient operating conditions with the same reactor power level dropping to different degrees. Set the duration of a xenon transient step to 1 hour. Calculate the change of reactor axial power shape index (ASI) over time and the corresponding physical parameters for each group of operating conditions under 38 xenon transient steps, and output the corresponding results.

[0013] S2. Set up a proxy model and input the changes in the reactor axial power shape index (ASI) over time and the corresponding physical parameters for each operating condition in S1 into the proxy model.

[0014] The physical parameters include the reactor core height. H Decay constants of iodine and xenon nuclides, and diffusion coefficients of the fast and hot groups of the entire reactor. D 1. D 2. Absorption cross section Σ of fast and hot groups in the entire reactor a1 Σ a2 Σ fission cross section of the entire fast group and hot group of the reactor f1 Σ f2 The total removed section Σ R The number of neutrons produced during each fission in the reactor n and the reaction cross section of radionuclide xenon in thermal groups s Xe And calculate the average value of the above parameters after heap-wide averaging;

[0015] S3. Define the search parameters: iodine and xenon yields. c I , c Xe Average nuclear concentrations of iodine and xenon I , Xe The average neutron flux of the fast and hot groups of the entire reactor (I and II) and the reactor power coefficient. α p The range of values ​​for ;

[0016] S4. The parameter search part of the execution agent model is performed. The agent model performs a global search and dimensionality reduction of the parameters for the five sets of input working conditions, and performs local correction of the search parameters for each input working condition.

[0017] S5. Determine whether the relative standard deviation of the final output search parameters of the five working conditions is within 5%. If yes, proceed to S6; otherwise, repeat the S4 process.

[0018] S6, the proxy model combined with the final output search parameter values ​​from S4, plots the corresponding ASI- for each working condition. t Curve, and based on ASI under this operating condition. t The original curve undergoes spatiotemporal translation and scaling transformations.

[0019] S7. Determine the improvement rate of MSE for each working condition before and after spatiotemporal translation and scaling corrections. Is it above 10%?

[0020] S8. The xenon oscillation process curve prediction part of the execution agent model is used to establish the power transient value ΔP and E by combining the input operating conditions and the output parameters.Pd (0), α p and Xe And the time shift value ΔP between each pair of operating conditions t Spatial translation value Δ y With scaling value α The linear interpolation expression;

[0021] S9. Input the ΔP value of the working condition to be predicted, and obtain the predicted working condition through the fitting relationship obtained in S8. E Pd (0) Xe , α p Parameters and time shift value Δ t Spatial translation value Δ y With scaling value α The predicted value; input the predicted value into the ASI calculation model to obtain the ASI value under the predicted working condition. t curve.

[0022] Furthermore, a preferred embodiment is provided, wherein the method for determining the change of the reactor axial power shape index (ASI) over time and the corresponding physical parameters for each operating condition in S1 is as follows:

[0023]

[0024] In the formula, P B and P T These represent the power values ​​of the lower and upper halves of the reactor axis, respectively.

[0025] Furthermore, in a preferred embodiment, S2 further includes calculating the xenon transient operating condition determination curve. t The steps for determining the ASI value at time =0, and the xenon transient condition determination curve in t The ASI value at time =0 is named E P (0) and derivative E Pd (0), the derivative value is the change in ASI during the first xenon transient step, and its calculation expression is:

[0026]

[0027] Furthermore, a preferred embodiment is provided in which the search parameters for iodine and xenon yields are defined in S3. c I , c Xe Average nuclear concentrations of iodine and xenon I ,Xe The average neutron flux of the fast and hot groups of the entire reactor (I and II) and the reactor power coefficient. α p The value range does not exceed 0.01 to 100 times the actual measured value inside the reactor.

[0028] Furthermore, in a preferred embodiment, S4 also includes a surrogate model for the time- and space-varying variables of the average neutron flux 1 and 2 of the whole-reactor fast group and hot group, and the average nuclear concentrations of iodine and xenon. I , Xe The steps of introducing axial Fourier expansion and axial difference parameters,

[0029] Their expressions are as follows:

[0030]

[0031]

[0032] In the formula, i =1,2,3,4 represent 1, 2, 3, 4 respectively. I , Xe Four variables, n This represents the order of the Fourier expansion coefficients. b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions; E i This represents the axial difference function of the second-order Fourier expansion coefficients of the variable.

[0033] Furthermore, a preferred embodiment is provided, wherein the method for local correction of the search parameters for each input working condition in S4 is as follows: the parameters of the working condition are globally searched and locally corrected based on the differential evolution-SLSQP hybrid algorithm.

[0034] Furthermore, a preferred embodiment is provided, wherein the method for performing spatiotemporal translation and scaling transformation based on the original ASI-t curve under this operating condition in S6 is as follows:

[0035] The expression for the transformation function is defined as follows:

[0036]

[0037] in, These represent the time transformation values ​​of the curve, with positive values ​​indicating leftward shift and spatial transformation values, respectively. The default positive values ​​are for leftward shift and scaling transformations, with 1 as the base value. Each value in the curve is multiplied by this scaling value to obtain the scaled new curve. The transformed ASI value;

[0038]

[0039] Furthermore, a preferred embodiment is provided, wherein in S7, the improvement rate of MSE for each set of working conditions before and after spatiotemporal translation and scaling corrections is determined. The method to determine whether it is above 10% is as follows:

[0040]

[0041] If so, output the corrected time shift value Δ t Spatial translation value Δ y With scaling value If not, output the time shift value Δ. t Spatial translation value Δ y The result is 0, scaling value The result is 1.

[0042] Option 3: A computer device, including a memory and a processor, wherein the memory stores a computer program, and when the processor runs the computer program stored in the memory, the processor executes the method described in any one of Options 1.

[0043] Option 4: A computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the method described in any one of Options 1.

[0044] The advantages of this invention are:

[0045] The present invention provides a method for predicting xenon oscillations in nuclear reactors based on data assimilation and physical constraints. This method establishes a surrogate model suitable for predicting xenon oscillation curves, combining the advantages of physical modeling and data-driven methods. It exhibits strong adaptability to different reactor types and can obtain prediction results without extensive computation; thus, it enables the prediction of xenon oscillation evolution within a certain range. The method described in this invention combines the advantages of physical modeling and data-driven approaches, and the results obtained have certain reference value.

[0046] The surrogate model described in this invention automatically reduces the dimensionality of the required input parameters, thereby reducing the complexity of subsequent variable fitting and also reducing the amount of computation.

[0047] The surrogate model described in this invention has low computational cost and high computational efficiency. On average, it can obtain the search parameter results for five sets of working conditions in 2-3 minutes, and the results of the subsequent parameter fitting and prediction part can be obtained in less than 10 seconds.

[0048] The proxy model described in this invention can predict the evolution trend of xenon oscillation at any point within the operating range using several operating conditions. This can greatly simplify the calculation of core power capability analysis and provide a certain reference for power changes during the xenon oscillation process in the core.

[0049] This invention is also applicable to the field of surrogate model technology for predicting xenon oscillation curves. Attached Figure Description

[0050] Figure 1 This is a flowchart of the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints as described in Implementation Method 1.

[0051] Figure 2 The original ASI- for the five calculation conditions with an initial power level of 100% FP as described in Implementation Method Eleven. t Schematic diagram of the curve.

[0052] Figure 3 For the ASI-based prediction of the initial power level of 100% FP as described in Implementation Method Eleven t Curve comparison chart.

[0053] Among them, (a) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 100%FP to 65%FP; (b) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 100%FP to 55%FP; (c) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 100%FP to 45%FP; (d) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 100%FP to 35%FP; (e) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 100%FP to 25%FP; and (f) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 100%FP to 15%FP.

[0054] Figure 4 The original ASI- for the five calculation conditions with an initial power level of 90% FP as described in Implementation Method Eleven. t Schematic diagram of the curve.

[0055] Figure 5 For the ASI-based prediction of the initial power level of 90% FP as described in Implementation Method Eleven t Curve comparison chart.

[0056] Among them, (a) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 90%FP to 65%FP; (b) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 90%FP to 55%FP; (c) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 90%FP to 45%FP; (d) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 90%FP to 35%FP; (e) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 90%FP to 25%FP; and (f) is a comparison of the predicted curve and the actual curve when the initial power level suddenly drops from 90%FP to 15%FP. Detailed Implementation

[0057] To make the objectives, technical solutions, and advantages of the embodiments of this application clearer, the technical solutions of the embodiments of this application will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them.

[0058] Implementation Method 1: This implementation method provides a method for predicting xenon oscillations in nuclear reactors based on data assimilation and physical constraints. The method includes the following steps:

[0059] S1. Select five groups of xenon transient operating conditions with the same reactor power level dropping to different degrees. Set the duration of a xenon transient step to 1 hour. Calculate the change of reactor axial power shape index (ASI) over time and the corresponding physical parameters for each group of operating conditions under 38 xenon transient steps, and output the corresponding results.

[0060] S2. Set up a proxy model and input the changes in the reactor axial power shape index (ASI) over time and the corresponding physical parameters for each operating condition in S1 into the proxy model.

[0061] The physical parameters include the reactor core height. H Decay constants of iodine and xenon nuclides, and diffusion coefficients of the fast and hot groups of the entire reactor. D 1. D 2. Absorption cross section Σ of fast and hot groups in the entire reactor a1 Σ a2 Σ fission cross section of the entire fast group and hot group of the reactor f1 Σ f2 The total removed section Σ R The number of neutrons ν produced by each fission in the reactor and the reaction cross section of the radionuclide xenon in the thermal group. s Xe And calculate the average value of the above parameters after heap-wide averaging;

[0062] S3. Define the search parameters: iodine and xenon yields. cI , c Xe Average nuclear concentrations of iodine and xenon I , Xe The average neutron flux of the fast and hot groups of the entire reactor (I and II) and the reactor power coefficient. α p The range of values ​​for ;

[0063] S4. The parameter search part of the execution agent model is performed. The agent model performs a global search and dimensionality reduction of the parameters for the five sets of input working conditions, and performs local correction of the search parameters for each input working condition.

[0064] S5. Determine whether the relative standard deviation of the final output search parameters of the five working conditions is within 5%. If yes, proceed to S6; otherwise, repeat the S4 process.

[0065] S6, the proxy model combined with the final output search parameter values ​​from S4, plots the corresponding ASI- for each working condition. t Curve, and based on ASI under this operating condition. t The original curve undergoes spatiotemporal translation and scaling transformations.

[0066] S7. Determine the improvement rate of MSE for each working condition before and after spatiotemporal translation and scaling corrections. Is it above 10%?

[0067] S8. The xenon oscillation process curve prediction part of the execution agent model, combining the input operating conditions and the output parameters, establishes the power transient value ΔP and... E Pd (0), α p and Xe And the time shift value ΔP between each pair of operating conditions t Spatial translation value Δ y With scaling value α The linear interpolation expression;

[0068] S9. Input the ΔP value of the working condition to be predicted, and obtain the predicted working condition through the fitting relationship obtained in S8. E Pd (0) Xe , α p Parameters and time shift value Δ t Spatial translation value Δ y With scaling value α The predicted value; input the predicted value into the ASI calculation model to obtain the ASI value under the predicted working condition. t curve.

[0069] Implementation Method Two: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method One. The method for the change of the reactor axial power shape index (ASI) over time and the corresponding physical parameters for each operating condition in S1 is as follows:

[0070]

[0071] In the formula, P B and P T These represent the power values ​​of the lower and upper halves of the reactor axis, respectively.

[0072] Implementation Method 3: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method 2. S2 also includes the calculated xenon transient condition determination curve. t The steps for determining the ASI value at time =0, and the xenon transient condition determination curve in t The ASI value at time =0 is named E P (0) and derivative E Pd (0), the derivative value is the change in ASI during the first xenon transient step, and its calculation expression is:

[0073] .

[0074] Implementation Method Four: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method One. In S3, the search parameters iodine and xenon yields are defined. c I , c Xe Average nuclear concentrations of iodine and xenon I , Xe The average neutron flux of the fast and hot groups of the entire reactor (I and II) and the reactor power coefficient. α p The value range does not exceed 0.01 times to 100 times the actual measured value inside the reactor.

[0075] Implementation Method 5: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method 1. S4 also includes a surrogate model for the time- and space-varying variables: the average neutron fluxes 1 and 2 of the whole reactor fast and hot groups, and the average nuclear concentrations of iodine and xenon. I , Xe The steps of introducing axial Fourier expansion and axial difference parameters,

[0076] Their expressions are as follows:

[0077]

[0078]

[0079] In the formula, i =1,2,3,4 represent 1, 2, 3, 4 respectively. I , Xe Four variables, n This represents the order of the Fourier expansion coefficients. b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions; E i This represents the axial difference function of the second-order Fourier expansion coefficients of the variable.

[0080] Implementation Method Six: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method Five. The method for local correction of search parameters for each input operating condition in S4 is as follows: global search and local correction of the parameters of the operating condition are performed based on the differential evolution-SLSQP hybrid algorithm.

[0081] Implementation Method Seven: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method Two. In S6, based on the ASI under this operating condition... t The methods for performing spatiotemporal translation and scaling transformations on the original curve are as follows:

[0082] The expression for the transformation function is defined as follows:

[0083]

[0084] in, These represent the time transformation values ​​of the curve, with positive values ​​indicating leftward shift and spatial transformation values, respectively. The default positive values ​​are for leftward shift and scaling transformations, with 1 as the base value. Each value in the curve is multiplied by this scaling value to obtain the scaled new curve. The transformed ASI value;

[0085] .

[0086] Implementation Method Eight: This implementation method further defines the nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints described in Implementation Method Two. In S7, the improvement rate of MSE before and after spatiotemporal translation and scaling corrections for each operating condition is determined. The method to determine whether it is above 10% is as follows:

[0087]

[0088] If so, output the corrected time shift value Δ t Spatial translation value Δ y With scaling value If not, output the time shift value Δ. t Spatial translation value Δ y The result is 0, scaling value The result is 1.

[0089] Implementation Method Nine: A computer device, including a memory and a processor, wherein the memory stores a computer program, and when the processor runs the computer program stored in the memory, the processor executes the method described in any one of Implementation Methods One to Eight.

[0090] Implementation Method 10: A computer-readable storage medium storing a computer program that, when executed by a processor, implements the steps of the method described in any one of Implementation Methods 1 to 8.

[0091] Implementation Method Eleven: The embodiments presented in this implementation method are used to explain the above-described Implementation Methods One to Ten, and specifically include the following:

[0092] The purpose of this implementation is to establish a physically constrained xenon oscillation surrogate model based on data assimilation. The surrogate model consists of two parts: a parameter search part and a xenon oscillation process curve prediction part. The parameter search part refers to the model searching for globally optimal parameters for multiple operating conditions using a data assimilation method under physical constraints. The xenon oscillation process curve prediction part refers to using the fitting relationship between parameters to predict the xenon oscillation process under unknown operating conditions.

[0093] This implementation method specifically includes the following steps:

[0094] S1: Select five groups of xenon transient operating conditions with the same reactor power level but varying degrees of sudden drop. Set the xenon transient step duration to 1 hour. Based on a traditional nuclear design program, calculate the change of reactor axial power shape index (ASI) over time for each group of operating conditions under 38 xenon transient steps, along with the corresponding physical parameters, and output the results. The ASI calculation expression is:

[0095]

[0096] In the formula, P B and P T These represent the power values ​​of the lower and upper halves of the reactor axis, respectively.

[0097] S2: Input the reactor axial power shape index (ASI) over time and the corresponding physical parameters for each operating condition into the proxy model.

[0098] The required physical parameters include reactor core height. H Decay constants of iodine and xenon nuclides, and diffusion coefficients of the fast and hot groups of the entire reactor. D 1. D 2. Absorption cross section Σ of fast and hot groups in the entire reactor a1 Σ a2 Σ fission cross section of the entire fast group and hot group of the reactor f1 Σ f2 The total removed section Σ R The number of neutrons produced during each fission in the reactor n and the reaction cross section of radionuclide xenon in thermal groups s Xe .

[0099] These parameters are all averaged values ​​obtained after averaging across the entire stack. At the same time, it is also necessary to determine the curve based on the calculated xenon transient operating conditions. t The ASI value at time =0 (named) E P (0) and derivative E Pd (0), the derivative value is the change in ASI during the first xenon transient step, and its calculation expression is:

[0100]

[0101] S3: Define the dimensions and value range of the search parameters. The number of search parameters is the dimension of the search parameters, typically taken as the yield of iodine and xenon. c I , c Xe Average nuclear concentrations of iodine and xenon I , Xe The average neutron flux of the fast and hot groups of the entire reactor (I and II) and the reactor power coefficient. α p The upper and lower limits of the search parameter value range can be set to a maximum of 0.01 times and 100 times the actual measured value inside the reactor, respectively. The specific value range of the parameter needs to be adjusted based on the results of the model algorithm.

[0102] S4: Execute the parameter search part of the surrogate model. After running the surrogate model, it will combine the input ASI- t The curve and the range of search parameters are used for a global search, so that the search results obtained from the input can be plotted according to the formulas in the model to match the original ASI for the corresponding working condition. tCurves that are similar in shape can be used to reproduce the effect.

[0103] ASI of the proxy model t The calculation function is derived based on the two-group one-dimensional diffusion model of reactor physics and the physical model of iodine and xenon concentration changes, as shown in expression 3-6:

[0104]

[0105] The surrogate model considers time- and space-varying variables such as the average neutron flux 1 and 2 of the fast and hot groups of the entire reactor, and the average nuclear concentrations of iodine and xenon. I , Xe Introducing the axial Fourier expansion and axial difference parameter, their expressions are as follows:

[0106]

[0107]

[0108] In the formula, i =1,2,3,4 represent 1, 2, 3, 4 respectively. I , Xe Four variables, n This represents the order of the Fourier expansion coefficients. b i,2 (∞) represents the value of the second-order Fourier expansion coefficient of the variable under equilibrium conditions. E i This represents the axial difference function of the second-order Fourier expansion coefficients of the variable.

[0109] This leads to a homogeneous system of equations for the four axial difference parameters:

[0110]

[0111]

[0112] In the formula,

[0113]

[0114]

[0115]

[0116]

[0117]

[0118]

[0119]

[0120]

[0121]

[0122]

[0123]

[0124]

[0125]

[0126] Finally, the relationship between the axial power shape index (ASI) and time is obtained through Laplace transform:

[0127]

[0128] In the formula, E P This refers to the reactor axial power shape index (ASI), and has

[0129]

[0130]

[0131]

[0132]

[0133]

[0134]

[0135] Therefore, the ASI of the proxy model t The curve will be plotted based on equation (13).

[0136] Let the vector x represent the combination of seven search parameter values. y j,1 For the original ASI input t Curve value, y j,2 (x) represents the curve value obtained after inputting seven search parameter values, and MSE represents the original ASI- t The mean square error between the curve value and the calculated new curve value is expressed as follows:

[0137]

[0138]

[0139] In the formula, A This represents the number of samples, and in the model, it represents the number of xenon transient steps.

[0140] The differential evolution-SLSQP hybrid algorithm of the surrogate model searches for the optimal x-vector value that minimizes the mean squared error (MSE). The principle of the hybrid algorithm is as follows:

[0141] 1. Population initialization

[0142]

[0143] in, NP The population size is set to 50 in this model, but can be adjusted based on the model results. ln and one These represent the values ​​of the x vector, respectively. n The minimum and maximum values ​​within the range of values ​​for each parameter.

[0144] 2. Mutation Strategy

[0145] After initializing the population, a mutation strategy is used to induce mutation. The best / 1 mutation strategy is selected, and its expression is as follows:

[0146]

[0147]

[0148] in, G To determine the maximum number of generations, this model uses 150. For the value of the mutated individual, For the optimal individual value, F ( G ) is a scaling factor that varies with the maximum number of generations.

[0149] 3. Adaptive crossover operation

[0150] After mutation, perform a total cross between the mutated vector and the target vector, as shown in the following expression:

[0151]

[0152] in, CR This represents the crossover operator; in this model, it is set to 0.7. Indicates the first m Population, number n The search parameter, the first G +1 individual value of the population.

[0153] 4. Dynamic range shrinkage

[0154] During the process of population evolution, the parameter range also continuously shrinks, and its expression is as follows:

[0155]

[0156] in, and They represent the first G The generation n The average value of all populations under each search parameter. and These are the upper and lower limits of the new search parameter's value range, respectively.

[0157] 5. Individual local correction

[0158] The optimal variable vector x obtained through the aforementioned process best Apply neighborhood constraints to achieve better results for each operating condition. The constraints are as follows:

[0159]

[0160] The constraint equations are solved using the BFGS quasi-Newton method.

[0161] After inputting the physical parameters and search parameter ranges for the corresponding operating conditions, the model will perform a global search and local correction of the parameters for the five operating conditions based on a differential evolution-SLSQP hybrid algorithm. Under the algorithm's effect, parameters that are insensitive to the model (generally...) c I , c Xe , I (1) and (2) will yield equal results, achieving dimensionality reduction. In other words, as long as the parameters corresponding to these five working conditions have the same value, the original curve can be reproduced. The hybrid algorithm applies sensitive parameters (i.e., ...) to each input working condition. Xe , α p Local corrections are made to reduce the mean square error (MSE) and obtain search results that better match the original operating conditions.

[0162] S5: Determine whether the relative standard deviation of the final output search parameters of the five working conditions is within 5%. If yes, proceed to the next step; otherwise, repeat the S4 process.

[0163] S6: The surrogate model combines the search parameter values ​​output from S4 to plot the ASI under the corresponding working condition. t Curve, and based on ASI under this operating condition. t The original curve undergoes spatiotemporal translation and scaling transformations.

[0164] The expression for the transformation function is defined as follows:

[0165]

[0166] in, These represent the time transformation value (the default positive value is a leftward shift), the spatial transformation value (the default positive value is a leftward shift), and the scaling transformation value (based on 1, each value in the curve is multiplied by this scaling value to obtain the new scaled curve). This is the transformed ASI value.

[0167] The range of the three values ​​is as follows:

[0168]

[0169] The expression for the mean square error between the transformed curve value and the original curve value is:

[0170]

[0171] The new curve obtained after the above transformation is... t The initial value at =0 is E P The constraint of (0) and the effect of further reducing MSE.

[0172] S7: Determine the improvement rate of MSE for each working condition before and after spatiotemporal translation and scaling corrections. Is it above 10%? The calculation formula is as follows:

[0173]

[0174] If so, output the corrected time shift value Δ t Spatial translation value Δ y With scaling value If not, output the time shift value Δ. t Spatial translation value Δ y The result is 0, scaling value The result is 1.

[0175] S8: Execute the Xenon Oscillation Process Curve Prediction Section. To predict the ASI change over time under xenon transient conditions with the same reactor power level but different degrees of sudden drops, a new variable, the power transient value ΔP, needs to be introduced. This variable represents the magnitude of the sudden drop in reactor power level under this condition (e.g., if the reactor power level drops from 100% FP to 40% FP, then the value of ΔP is 60). Then, a relationship is established between the power transient value ΔP and... E Pd (0) Xe ,α p The relationship between ΔP and... E Pd (0), α p and Xe ΔP and Xe These relationships approximate linear, linear, and exponential changes, respectively, thus allowing the establishment of corresponding fitting formulas to calculate all relevant parameters. For the time shift value Δ output in S7... t Spatial translation value Δ y With scaling value α The prediction is approximately made by using linear fitting between pairs of operating conditions.

[0176] S9: Input the ΔP value of the working condition to be predicted, and obtain the predicted working condition through the fitting relationship obtained in S8. E Pd (0) Xe , α p Parameters and time shift value Δ t Spatial translation value Δ y With scaling value α The predicted value is input into the ASI calculation model to obtain the ASI value under the predicted operating condition. t curve.

[0177] Example 1: Five sets of xenon transient operating conditions were selected during the initial fuel loading of the HPR1000 reactor, with an initial power level of 100% FP and power transient changes ΔP of 30, 40, 60, 80, and 90. The ASI of the five sets of operating conditions was calculated. -t Curve and corresponding initial physical parameters

[0178] Table 1 Initial physical parameter values

[0179]

[0180] Using a proxy model to perform parameter search, the following results were obtained:

[0181] Table 2 Final values ​​of search parameters

[0182]

[0183] Table 3 Transformation Parameter Values

[0184]

[0185] The fitting relationships obtained from Table 2 are as follows:

[0186] E Pd(0) = 5.69526E-7 - 4.02001E-7 *ΔP

[0187] Xe = 7.32448E16 - 1.24041E18 * α p

[0188] Xe = 8.63005E13 * exp(ΔP / 15.09661) + 1.23629E16

[0189] Δ t , α With Δ y The power transient ΔP is obtained using linear interpolation.

[0190] By inputting ΔP = 35,45,55,65,75,85 and calculating the corresponding parameters, plot the ASI- t The prediction curve is shown in Figure 3.

[0191] Example 2: Five sets of xenon transient operating conditions were selected during the initial fuel loading of the HPR1000 reactor, with an initial power level of 90% FP and power transients ΔP values ​​of 30, 40, 60, 70, and 80. The ASI values ​​for these five operating conditions were calculated using a conventional nuclear design program. t Curves and corresponding initial physical parameters, such as Figure 4 As shown in Table 4.

[0192] Table 4 Initial physical parameter values

[0193]

[0194] Using a proxy model to perform parameter search, the following results were obtained:

[0195] Table 5 Final values ​​of search parameters

[0196]

[0197] Table 6 Transformation Parameter Values

[0198]

[0199] The fitting relationships obtained from Table 5 are as follows:

[0200] E Pd (0) = 6.79992E-7 - 4.02487E-7 *ΔP

[0201] Xe= 7.4254E16 - 1.2372E18 * α p

[0202] Xe = 1.22511E14 * exp(ΔP / 14.8749) + 1.05896E16

[0203] Δ t , α With Δ y The power transient ΔP is obtained using linear interpolation.

[0204] By inputting ΔP = 25,35,45,55,65,75 and calculating the corresponding parameters, plot the ASI- t The prediction curve is as follows Figure 5 As shown.

[0205] Those skilled in the art will understand that the above description is merely a preferred embodiment of the present invention, and the features described in the various embodiments and / or claims of this disclosure can be combined or combined in various ways, even if such combinations or combinations are not explicitly described in this disclosure. This is not intended to limit the present invention. Although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art can still modify the technical solutions described in the foregoing embodiments or make equivalent substitutions for some of the technical features. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the protection scope of the present invention.

[0206] Although preferred embodiments of the invention have been described, those skilled in the art, upon learning the basic inventive concept, can make other changes and modifications to these embodiments. Therefore, the appended claims are intended to be interpreted as including both the preferred embodiments and all changes and modifications falling within the scope of the invention. Clearly, those skilled in the art can make various alterations and modifications to the invention without departing from its spirit and scope. Thus, if these modifications and modifications of the invention fall within the scope of the claims and their equivalents, the invention is also intended to include these modifications and modifications.

Claims

1. A method for predicting xenon oscillations in a nuclear reactor based on data assimilation and physical constraints, the method comprising: The method comprises the following steps: S1, selecting five groups of xenon transient conditions with different degrees of sudden drop in reactor power level, setting a xenon transient step length of 1h, calculating the change of reactor axial power shape index ASI with time and the corresponding physical parameters of each group of conditions under 38 xenon transient steps, and outputting the corresponding results; S2, setting a surrogate model, inputting the change of reactor axial power shape index ASI with time and the corresponding physical parameters of each group of conditions in S1 into the surrogate model; said physical parameters include the reactor core height H , the decay constants of iodine and xenon nuclides, the diffusion coefficients of the fast and thermal groups throughout the reactor D 1、 D 2, the absorption cross sections of the fast and thermal groups throughout the reactor a1 , Σ a2 , the fission cross sections of the fast and thermal groups throughout the reactor f1 , Σ f2 , the removal cross sections throughout the reactor R , the number of neutrons produced per fission in the reactor, ν, and the reaction cross section of the xenon nuclide in the thermal group σ Xe ; and the average values obtained by averaging said parameters throughout the reactor are calculated; S3, defining search parameters yield of iodine and xenon γ I 、γ Xe , average nuclear concentration of iodine and xenon I 、 Xe full core fast group and thermal group average neutron flux 1 、 2 and power coefficient within the reactor α p value range; S4, executing the parameter search part of the surrogate model, the surrogate model performing global search and dimension reduction of parameters for the input five groups of conditions, and performing local correction of search parameters for each input condition; S5, judging whether the relative standard deviation of the final output search parameters of the five groups of conditions is within 5%, if yes, proceeding to S6, if not, repeating S4; S6, the agent model combines the search parameter values output in S4 to draw the ASI- t curve corresponding to each group of working conditions, and performs time-space translation transformation and scaling transformation based on the ASI- t curve under the working condition. S7, judging the promotion rate of MSE before and after the space-time translation transformation correction and the scaling correction of each group of working conditions whether above 10%; S8, the xenon oscillation process curve prediction part of the agent model, combined with the input working condition and the output parameter result, establishes the power transient value ΔP and E Pd (0), α p With Xe and the time translation value Δ t , the space translation value Δ y and the linear interpolation expression of the scaling value α of ΔP and E between two working conditions S9, input the ΔP value of the working condition to be predicted, and obtain the predicted value of the working condition to be predicted through the fitting relationship obtained in S8 E Pd (0)、 Xe 、 α p parameters and time translation value Δ t , space translation value Δ y and scaling value α of the predicted value; input the predicted value into the ASI calculation model to obtain the ASI curve under the predicted working condition. t ​ 2. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, The method for the change of reactor axial power shape index ASI with time and the corresponding physical parameters of each group of conditions in S1 is: wherein P B with P T respectively represent the power values of the lower and upper halves of the reactor axis.

3. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, The step of calculating the ASI value at the time of t = 0 of the xenon transient condition determination curve is also included in S2, and the ASI value at the time of t = 0 of the xenon transient condition determination curve is named as ASI (0). t t E P (0) and the derivative E Pd (0), the derivative value being the amount of change of ASI in the first xenon transient step time, and the calculation expression being:​​ 。 4. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, S3 defines the search parameters yield of iodine and xenon γ I 、γ Xe , average nuclear concentration of iodine and xenon I 、 Xe , average neutron flux of the full core fast group and thermal group 、 2 and the power coefficient in the reactor α p The value range of the above-mentioned parameters does not exceed 0.01 times to 100 times of the real measured value in the reactor core.

5. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, S4 also includes surrogate models for time- and space-varying variables, such as the average neutron fluxes 1 and 2 of the whole-reactor fast and hot groups, and the average nuclear concentrations of iodine and xenon. I , Xe The steps of introducing axial Fourier expansion and axial difference parameters, The expressions are respectively: wherein i = 1, 2, 3, 4, respectively, represent 1, 2, I , Xe four variables, n represents the order of the Fourier expansion coefficient, b i,2 (∞) represents the value of the second order Fourier expansion coefficient of the variable at equilibrium conditions; E i represents the axial difference function of the second order Fourier expansion coefficient of the variable.

6. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, The method for local correction of search parameters for each input condition in S4 is: based on the differential evolution-SLSQP hybrid algorithm, the parameters of the condition are globally searched and locally corrected.

7. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, The method for time-space translation transformation and scaling transformation based on the ASI-t original curve of the condition in S6 is: The expression of the transformation function is defined as: wherein, respectively represent time transform value of the curve, the default positive value is left shift, space transform value, the default positive value is left shift and scaling transform value, wherein 1 is the reference value, each value in the curve is multiplied by the scaling value to obtain the new curve after scaling, is the transformed ASI value; 。 8. The nuclear reactor xenon oscillation prediction method based on data assimilation and physical constraints of claim 1, wherein, S7, the improvement rate of MSE of each group of working conditions before and after the space-time translation transformation correction and the scaling correction is determined The method for determining whether the percentage is more than 10% is: If so, output the modified time-panning value Δ t , spatial-panning value Δ y , and scaling value ; if not, output the time-panning value Δ t , spatial-panning value Δ y , and scaling value whose result is 0.

9. Computer device comprising a memory and a processor, characterized in that The memory stores a computer program, and when the processor runs the computer program stored in the memory, the processor executes the method of any one of claims 1-8.

10. A computer readable storage medium characterized by, The computer readable storage medium stores a computer program, and the computer program is executed by the processor to realize the steps of the method of any one of claims 1-8.

Images