A semiconductor vertical furnace multi-process software control method and system

Abstract

The application relates to the technical field of semiconductor manufacturing equipment control, and discloses a semiconductor vertical furnace multi-process software control method and system. The method acquires historical temperature, gas flow and cavity pressure data and extracts relevant time sequence features; cross-physical quantity coupling analysis is performed when disturbance propagation is judged, and a disturbance propagation path is determined; real-time heat exchange efficiency is calculated according to the disturbance propagation path, and future temperature distribution is predicted; when the difference between the predicted temperature distribution and a target temperature distribution is greater than a temperature deviation threshold, control variables are reversely solved, and each temperature zone heater power and a mass flow controller are coordinately adjusted. The application can reduce compensation lag and improve the temperature control precision of the vertical furnace.

Description

Technical Field

[0001] This invention relates to the field of semiconductor manufacturing equipment control technology, and in particular to a multi-process software control method and system for semiconductor vertical furnaces. Background Technology

[0002] Currently, vertical furnaces are crucial heat treatment equipment in semiconductor manufacturing, primarily used for processes such as wafer oxidation, diffusion, doping, and annealing. During furnace operation, a SCADA system is required to collect and control parameters such as temperature in different zones within the furnace tube, process gas flow rate, and chamber pressure in real time, ensuring that the wafers are processed under preset process conditions.

[0003] In a prior art, a vertical furnace is typically divided into multiple temperature zones along its axial direction, and the heater power of each zone is independently controlled to ensure that the actual temperature of each zone is close to its corresponding target temperature. However, the temperature zones inside a vertical furnace are not independent of each other; heat transfer occurs between adjacent zones through the furnace wall, thermal radiation, and gas flow. When the heater power of one zone changes, the generated heat may be transferred to other zones, causing changes in the actual temperature of adjacent zones.

[0004] However, existing independent zone control logic combined with traditional SCADA monitoring architecture only adjusts based on the current temperature deviation of each zone, making it difficult to identify in advance the propagation process of disturbances caused by changes in parameters such as gas flow rate and chamber pressure. When the vertical furnace undergoes process switching, the control system usually only compensates after a significant temperature deviation has occurred, which easily leads to compensation lag and a longer time required for the temperature field to reach a stable state.

[0005] Therefore, existing technologies suffer from insufficient temperature control accuracy in vertical furnaces due to the difficulty in identifying the dynamic coupling relationship between temperature, gas flow rate, and cavity pressure. Summary of the Invention

[0006] This invention provides a multi-process software control method and system for semiconductor vertical furnaces to solve the technical problem of insufficient temperature control accuracy in vertical furnaces caused by the difficulty in identifying the dynamic coupling relationship between temperature, gas flow rate and cavity pressure in the prior art.

[0007] In a first aspect, to solve the above-mentioned technical problems, the present invention provides a multi-process software control method for a semiconductor vertical furnace, the method being executed by a computer and comprising:

[0008] Historical temperature data, historical gas flow data, and historical cavity pressure data of each temperature zone of the semiconductor vertical furnace during its historical operating cycle are obtained, and long-term and short-term time series analysis is performed to extract correlation time series features to characterize the dynamic correlation of multiple physical quantities.

[0009] When it is determined that disturbance propagation has occurred based on the associated time series characteristics, cross-physical quantity coupling analysis is performed on the associated time series characteristics to determine the disturbance propagation path;

[0010] Based on the disturbance propagation path, the heat transfer efficiency is calculated using the real-time gas flow rate data and real-time cavity pressure data at the current moment to obtain the real-time heat exchange efficiency.

[0011] Based on the real-time heat exchange efficiency and the real-time gas flow rate data, the temperature field at future times is predicted, and the predicted temperature distribution is obtained.

[0012] If the difference between the predicted temperature distribution and the preset target temperature distribution is greater than the preset temperature deviation threshold, then the control quantity used to coordinate the adjustment of the heater power and mass flow controller of each temperature zone of the semiconductor vertical furnace is solved in reverse based on the difference, and coordinated adjustment is performed based on the control quantity.

[0013] In a second aspect, the present invention provides a multi-process software control system for a semiconductor vertical furnace, including a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the above-described multi-process software control method for a semiconductor vertical furnace.

[0014] Compared with the prior art, the present invention has the following beneficial effects:

[0015] (1) This invention performs long-term and short-term time series analysis on historical temperature data, historical gas flow data and historical cavity pressure data, extracts correlation time series features to characterize the dynamic correlation of multiple physical quantities, and performs cross-physical quantity coupling analysis based on the correlation time series features to determine the disturbance propagation path. This enables the identification of the dynamic coupling relationship and disturbance propagation process between temperature, gas flow and cavity pressure, providing a data basis for subsequent real-time heat exchange efficiency calculation and temperature field prediction.

[0016] (2) Based on the disturbance propagation path, the present invention calculates the heat transfer efficiency of real-time gas flow data and real-time cavity pressure data to obtain the real-time heat exchange efficiency, and predicts the temperature distribution at future times based on the real-time heat exchange efficiency and the real-time gas flow data, thereby obtaining the temperature change trend before the temperature deviation further expands, reducing the compensation lag caused by adjusting only based on the current temperature deviation, and shortening the time required for the temperature field to reach a stable state.

[0017] (3) When the difference between the predicted temperature distribution and the target temperature distribution is greater than the temperature deviation threshold, the present invention determines the compensation heat flux density based on the difference, and determines the power adjustment amount of the heater in each temperature zone and the flow adjustment amount of the mass flow controller through reverse heat transfer, so as to coordinate the adjustment of the power of the heater in each temperature zone and the flow rate of the process gas, thereby improving the temperature control accuracy of the semiconductor vertical furnace. Attached Figure Description

[0018] Figure 1 This is a flowchart illustrating a multi-process software control method for a semiconductor vertical furnace provided in an embodiment of the present invention;

[0019] Figure 2 This is a schematic diagram of the partial correlation matrix provided in an embodiment of the present invention;

[0020] Figure 3 This is a schematic diagram of the predicted temperature distribution provided in an embodiment of the present invention. Detailed Implementation

[0021] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0022] Reference Figure 1 The first embodiment of the present invention provides a multi-process software control method for semiconductor vertical furnaces, comprising the following steps:

[0023] S1. Acquire historical temperature data, historical gas flow data, and historical cavity pressure data of each temperature zone of the semiconductor vertical furnace during its historical operating cycle, and perform long-term and short-term time series analysis to extract correlation time series features used to characterize the dynamic correlation of multiple physical quantities.

[0024] S2, When it is determined that disturbance propagation has occurred based on the associated time series characteristics, cross-physical quantity coupling analysis is performed on the associated time series characteristics to determine the disturbance propagation path;

[0025] S3. Based on the disturbance propagation path, calculate the heat transfer efficiency using the real-time gas flow rate data and real-time cavity pressure data at the current moment to obtain the real-time heat exchange efficiency.

[0026] S4. Based on the real-time heat exchange efficiency and the real-time gas flow data, predict the temperature field at future moments to obtain the predicted temperature distribution.

[0027] S5. If the difference between the predicted temperature distribution and the preset target temperature distribution is greater than the preset temperature deviation threshold, then the control quantity used to coordinate the adjustment of the heater power and mass flow controller of each temperature zone of the semiconductor vertical furnace is solved in reverse according to the difference, and coordinated adjustment is performed according to the control quantity.

[0028] In step S1, historical temperature data, historical gas flow data, and historical cavity pressure data of each temperature zone of the semiconductor vertical furnace during the historical operating cycle are obtained, and long-term and short-term time series analysis is performed to extract correlation time series features used to characterize the dynamic correlation of multiple physical quantities.

[0029] The historical temperature data consists of time-series temperature data collected by temperature sensors located in different temperature zones of the semiconductor vertical furnace during historical operating cycles; the historical gas flow data consists of time-series process gas flow data collected by flow sensors or mass flow controllers during corresponding historical operating cycles; and the historical chamber pressure data consists of time-series chamber pressure data collected by pressure sensors during corresponding historical operating cycles. These historical temperature, gas flow, and chamber pressure data can be acquired through the semiconductor vertical furnace's monitoring and data acquisition system and correlated according to corresponding timestamps and historical operating cycles to extract correlated time-series features.

[0030] In one implementation, the extraction of correlation time-series features for characterizing the dynamic correlation of multiple physical quantities includes:

[0031] The historical temperature data, historical gas flow data, and historical cavity pressure data are standardized, and the standardized data are segmented according to a preset sliding time window to obtain a multidimensional time series data sequence.

[0032] The multidimensional time-series data sequence is input into a preset long short-term memory network to generate the associated time-series features.

[0033] Specifically, historical temperature data, historical gas flow data, and historical cavity pressure data are time-aligned according to the acquisition timestamps corresponding to each historical data point. The sampling frequency of historical temperature data and historical cavity pressure data is set to 10 sampling points per second; when the original sampling frequency of historical gas flow data is 5 sampling points per second, linear interpolation is used to convert the historical gas flow data to 10 sampling points per second, so that the three types of historical data have the same sampling time interval.

[0034] Outlier removal and smoothing are performed on the time-aligned data from each sensor channel. In one implementation, an isolated forest algorithm is used to detect outlier sample values ​​in each sensor channel, and the detected outlier sample values ​​are replaced with interpolated results of adjacent normal sample values. A moving average algorithm with a window length of 5 sample points is used to smooth the time-series data after outlier removal. In this embodiment, the number of trees in the isolated forest algorithm is set to 100, and the contamination rate is set to 0.05.

[0035] It should be noted that the standardization process uses a single sensor channel as the statistical unit. Specifically, the temperature sensor channels corresponding to different temperature zones, the gas flow sensor channels corresponding to different gas branches, and the pressure sensor channels corresponding to different chamber pressure monitoring points are each treated as independent sensor channels.

[0036] For each sensing channel, based on all valid sampled data obtained after time alignment, outlier removal, and smoothing within multiple historical operating cycles used to train the Long Short-Term Memory network, the mean and standard deviation for that sensing channel are calculated separately. Therefore, different temperature zones, different gas branches, and different chamber pressure monitoring points each have their own corresponding mean and standard deviation, and multiple sensing channels under the same physical quantity type are not combined for statistical analysis.

[0037] Multiple historical operating cycles include a heating phase, a isothermal phase, a process reaction phase, and a cooling phase. Before training the Long Short-Term Memory (LSTM) network, the mean and standard deviation for each sensor channel are calculated and saved based on the training dataset. During training, the saved mean and standard deviation are used to standardize the training data. During real-time inference, the saved mean and standard deviation from the training phase are used again to standardize the real-time data of the corresponding sensor channels, without recalculating the mean and standard deviation based on the data within the current sliding time window. This ensures that the data from the same sensor channel in the training and real-time inference phases have a consistent standardization scale.

[0038] The standardized data from each sensor channel are synchronously segmented according to timestamps, so that the same multidimensional time-series data sequence contains temperature data of each temperature zone, flow data of each gas branch, and pressure data of each cavity within the same time range, thereby maintaining the temporal correspondence between different physical quantities.

[0039] Assume a semiconductor vertical furnace includes Temperature nodes in each temperature zone Each gas branch flow node and If there are 1 cavity pressure node, then the total number of physical quantity nodes is 1.

[0040]

[0041] The sliding time window length is set to 60 seconds, and the sliding step size is set to 10 seconds. Since the data from each sensor channel is standardized to 10 sampling points per second, the multidimensional time-series data sequence corresponding to each sliding time window includes 600 time steps, and its data dimension is... .

[0042] The first dimension of the multidimensional time series data sequence is the time axis, and the second dimension is the physical quantity node axis. The physical quantity nodes are arranged in a fixed order according to the temperature nodes of each temperature zone, the flow rate nodes of each gas branch, and the pressure nodes of each cavity. Different sliding time windows use the same arrangement order of physical quantity nodes.

[0043] In this embodiment, the preset long short-term memory network model includes a parameter-sharing long short-term memory network encoding layer, a one-dimensional convolutional layer, a global average pooling layer, and a prediction output layer used during the training phase.

[0044] The multidimensional time series data sequence is split according to the physical quantity node axis. Each node has a time series sequence, and each node time series sequence includes 600 time steps. The node time series sequences corresponding to each physical quantity node are respectively input into the same parameter-shared Long Short-Term Memory network coding layer for time series encoding.

[0045] The Long Short-Term Memory (LSTM) network encoding layer employs a single-layer unidirectional LSM network with an input dimension of 1 and a hidden state dimension of 128. For each physical quantity node, the hidden state vector output by the LSM network encoding layer at the 600th time step is selected as the independent temporal feature corresponding to the physical quantity node. Thus, each physical quantity node corresponds to a 128-dimensional hidden state vector.

[0046] The 128-dimensional hidden state vectors corresponding to each physical quantity node are arranged in a fixed order according to the physical quantity node axes to obtain the node-level temporal feature matrix. The node-level temporal feature matrix has the following number of rows. The number of columns is 128, meaning the dimension of the node-level time series feature matrix is ​​[value missing]. Each row uniquely corresponds to a temperature node in a temperature zone, a gas branch flow node, or a cavity pressure node, and each column corresponds to a hidden feature dimension output by the Long Short-Term Memory network encoding layer.

[0047] The node-level temporal feature matrix is ​​input into a one-dimensional convolutional layer for cross-physical quantity feature extraction. The one-dimensional convolutional layer performs convolution operations along the physical quantity node axis, with the number of input channels set to 128, the number of output channels set to 64, the kernel size set to 3, the stride set to 1, the padding size set to 1, and a linear rectified function used to perform non-linear processing on the convolution results.

[0048] Global average pooling is performed along the physical quantity node axes of the output of the one-dimensional convolutional layer to obtain a fixed 64-dimensional multi-physical quantity dynamic correlation feature. The multi-physical quantity dynamic correlation feature is used to characterize the overall dynamic correlation state between temperature, gas flow rate and cavity pressure within the current sliding time window.

[0049] The temperature values ​​of each temperature zone, the flow rate values ​​of each gas branch, and the pressure values ​​of each cavity at the next sampling time after each sliding time window are used as training targets. The 64-dimensional multi-physical quantity dynamic correlation features are mapped to a model containing these parameters through the prediction output layer. The prediction vector is a single-layer fully connected structure with an input dimension of 64 and an output dimension of N. No additional hidden layers or non-linear activation functions are used; the prediction vector is generated directly through linear mapping.

[0050] The training data consists of temperature, gas flow, and chamber pressure data collected over multiple historical operating cycles of the semiconductor vertical furnace, covering the heating, isothermal, reaction, and cooling stages. The training data is divided according to historical operating cycles, with 80% of the data from those cycles used as the training set and 20% as the validation set. Data from the same historical operating cycle is not simultaneously allocated to the training and validation sets. During model training, the training data covers at least 20 complete historical operating cycles, generating a total of at least 5000 sliding window samples. The sample sizes for the heating, isothermal, reaction, and cooling stages are naturally distributed according to the actual process time ratio of the semiconductor vertical furnace; no manual equalization sampling is performed for each stage.

[0051] The mean squared error between the predicted vector and the actual sampled vector at the next sampling time is used as the loss function for model training, and the adaptive moment estimation algorithm is used to update the model parameters. The initial learning rate for model training is set to 0.001, the batch size is set to 32, and the maximum number of training epochs is set to 100.

[0052] After each training epoch, the validation loss is calculated using the validation set. If the validation loss does not decrease after 10 consecutive training epochs, model training is stopped, and the model parameters of the long short-term memory network encoding layer and the one-dimensional convolutional layer corresponding to the lowest validation loss are saved.

[0053] After the model is trained, the prediction output layer is no longer used in the real-time inference process. Only the long short-term memory network encoding layer, one-dimensional convolutional layer and global average pooling layer are retained to generate node-level temporal feature matrix and 64-dimensional multi-physical quantity dynamic correlation features based on the real-time temperature data, real-time gas flow data and real-time cavity pressure data in the current sliding time window.

[0054] It should be noted that the associated time-series features include a node-level time-series feature matrix and 64-dimensional multi-physical quantity dynamic association features. The node-level time-series feature matrix is ​​used to retain the fixed mapping relationship between each physical quantity node and its corresponding hidden feature; the 64-dimensional multi-physical quantity dynamic association features are used to determine whether the preset disturbance propagation conditions are met. Specifically, the preset disturbance propagation conditions are determined according to the method in step S2: taking the current analysis time window and the previous 9 consecutive analysis time windows, a total of 10 analysis time windows, as a judgment period, the feature variance of each of the 64 feature dimensions is calculated along the time axis of the judgment period and the mean is taken to obtain the feature fluctuation variance; when the feature fluctuation variance is greater than the preset disturbance threshold, it is determined that the disturbance propagation conditions are met.

[0055] When it is determined that the disturbance propagation has occurred, the standardized time series data corresponding to each physical quantity node within the same sliding time window used when generating the associated time series features are used as the observation data for subsequent partial correlation analysis. The temperature nodes of each temperature zone, the flow nodes of each gas branch, and the pressure nodes of each cavity are used as the analysis variables for partial correlation analysis, instead of using the hidden feature dimensions output by the long short-term memory network as independent physical quantity nodes.

[0056] Therefore, a fixed correspondence is established between the time axis and physical quantity node axis in the multidimensional time series data sequence, the physical quantity node axis and hidden feature dimension in the node-level time series feature matrix, and the physical quantity nodes used in the subsequent cross-physical quantity coupling analysis.

[0057] In step S2, when it is determined that disturbance propagation has occurred based on the associated time-series characteristics, a cross-physical quantity coupling analysis is performed on the associated time-series characteristics to determine the disturbance propagation path, including:

[0058] A cross-physical quantity coupling analysis was performed on the associated time-series features to obtain the partial correlation matrix;

[0059] Extract elements whose absolute values ​​of partial correlation coefficients are greater than a preset correlation threshold from the partial correlation matrix, and construct a topological connection graph of physical quantity interaction;

[0060] The out-degree of each physical quantity node in the topology connection graph is counted, and the physical quantity node with the largest out-degree is determined as the disturbance source node;

[0061] The disturbance propagation path is determined based on the time lag between the disturbance source node and other physical quantity nodes, as well as the topology connection diagram.

[0062] The determination of whether disturbance propagation has occurred is based on the associated time-series characteristics corresponding to the continuous analysis time windows. Specifically, a preset number of continuous analysis time windows are arranged in chronological order to form a sequence of dynamic correlation characteristics of multiple physical quantities. For each feature dimension in the 64-dimensional dynamic correlation characteristics of multiple physical quantities, the variance of the feature value is calculated along the time axis formed by the continuous analysis time windows, obtaining the feature variances corresponding to the 64 feature dimensions respectively. The average of the 64 feature variances is used to obtain the feature fluctuation variance, which characterizes the overall fluctuation of the correlation time-series characteristics. The feature fluctuation variance is compared with a preset disturbance threshold. When the feature fluctuation variance is greater than the disturbance threshold, the correlation time-series characteristics are determined to meet the preset disturbance propagation conditions. The preset number is 10, that is, the current analysis time window and the previous 9 continuous analysis time windows, totaling 10 analysis time windows, are used as continuous analysis time windows.

[0063] The disturbance threshold is pre-calibrated based on historical operating data corresponding to normal process states and historical operating data corresponding to marked disturbance states. Specifically, the characteristic fluctuation variance distributions under normal process states and disturbance states are statistically analyzed separately, and a critical value that can distinguish between the two states is selected as the disturbance threshold. In this embodiment, the disturbance threshold is set to 5 and fixed before system operation. Specifically, the empirical cumulative distribution function of the characteristic fluctuation variances corresponding to normal process states and disturbance states is calculated separately, and the critical value is determined by maximizing the Youden exponent. The normal process state samples used for calibration have no less than 100 analysis time windows, and the disturbance state samples have no less than 20 analysis time windows corresponding to disturbance events.

[0064] After determining that disturbance propagation has occurred, standardized time-series data corresponding to each physical quantity node within the same analysis time window used to generate the associated time-series features are obtained, and partial correlation analysis is performed on the standardized time-series data. The time-series data corresponding to each physical quantity node are time-aligned at 10 sampling points per second, and each physical quantity node corresponds to 600 sample values ​​within a 60-second analysis time window.

[0065] In this system, each sampling time point is considered an observation sample, and each physical quantity node is considered an analysis variable, thus forming a data matrix with sampling time points as rows and physical quantity nodes as columns. When performing partial correlation analysis on any two physical quantity nodes, these two nodes are treated as the variables to be analyzed, while other nodes in the data matrix—temperature nodes in the temperature zone, gas branch flow rate nodes, and cavity pressure nodes—are used as control variables to eliminate the influence of other physical quantity nodes on the correlation between the two nodes. The partial correlation coefficient for any two physical quantity nodes i and j is calculated using the formula r. ij·rest =-Ω ij / √(Ω ii ×Ω jj), where Ω is the inverse of the covariance matrix of the standardized time series data for all N physical quantity nodes, Ω ij Ω represents the element in the i-th row and j-th column of the inverse matrix. ii and Ω jj These represent the i-th and j-th diagonal elements in the inverse matrix, respectively.

[0066] By performing pairwise calculations on each physical quantity node, a partial correlation matrix is ​​obtained. The rows and columns of the partial correlation matrix are arranged in a fixed order according to the temperature nodes of each temperature zone, the flow rate nodes of each gas branch, and the pressure nodes of each cavity. Each element in the matrix represents the degree of direct correlation between the corresponding two physical quantity nodes after excluding the influence of other physical quantity nodes.

[0067] The absolute values ​​of each partial correlation coefficient in the partial correlation matrix are compared with a preset correlation threshold. Elements whose absolute values ​​of partial correlation coefficients are not greater than the preset correlation threshold are set to zero, while elements whose absolute values ​​of partial correlation coefficients are greater than the preset correlation threshold are retained. Based on the retained non-zero elements, candidate connection relationships are determined between the corresponding physical quantity nodes, and a candidate topology connection graph is constructed.

[0068] The preset correlation threshold is pre-calibrated based on the distribution of partial correlation coefficients corresponding to normal operation data and labeled disturbance data, and is fixed and saved before system operation. In this embodiment, the preset correlation threshold is set to 0.6. It should be noted that at this time, the edges in the candidate topology connection graph only indicate the direct coupling relationship between the corresponding two physical quantity nodes, and do not yet indicate the direction of disturbance propagation. The preset correlation threshold is pre-calibrated based on the distribution of the absolute values ​​of the partial correlation coefficients corresponding to normal operation data and labeled disturbance data, and the critical value is determined by the selection criterion to maximize the accuracy of disturbance propagation path identification. The samples used for calibration include no less than 50 sets of partial correlation matrices under normal operation conditions and no less than 20 sets of partial correlation matrices under labeled disturbance conditions.

[0069] Cross-correlation analysis was performed on each pair of physical quantity nodes with candidate connectivity in the candidate topology connection graph. The cross-correlation analysis used standardized time-series data within the same analysis time window as the partial correlation analysis, with a sampling interval of 0.1 seconds and a lag search range of -30 seconds to +30 seconds.

[0070] Within the lag search range, the peak position with the largest absolute value of the cross-correlation coefficient is selected as the lag time difference between the corresponding two physical quantity nodes. When there are multiple identical maximum peaks, the peak with the smallest absolute lag time is selected; when the absolute lag times of multiple peaks are still the same, the peak with the same starting time sequence as the actual change of the two physical quantity nodes is selected.

[0071] For the standardized time-series data x of physical quantity nodes i and j i (t) and x jThe cross-correlation coefficient (t) under lag τ is defined as ρ. ij (τ)=Σ t (x i (t)-x̄ i (x) j (t+τ)-x̄ j ) / √(Σ t (x i (t)-x̄ i )²×Σ t (x j (t+τ)-x̄ j )²).

[0072] Let the cross-correlation function between the first physical quantity node and the second physical quantity node represent the correlation between the data of the first physical quantity node at the current moment and the data of the second physical quantity node at the lagging moment. When the lag time difference corresponding to the peak value is positive, it indicates that the change of the first physical quantity node precedes that of the second physical quantity node, and a directed edge is established from the first physical quantity node to the second physical quantity node; when the lag time difference is negative, a directed edge is established from the second physical quantity node to the first physical quantity node; when the lag time difference is zero, the corresponding nodes are marked as synchronous coupling.

[0073] In this embodiment, the out-degree of each physical quantity node in the topological connection graph is counted. Specifically, after determining the propagation direction of candidate connection relationships based on cross-correlation analysis, the number of directed edges sent from each physical quantity node to other physical quantity nodes in the formed directed topological connection graph is counted.

[0074] The physical quantity node with the largest out-degree is identified as the disturbance source node. When there are multiple physical quantity nodes with the same out-degree, the sum of the absolute values ​​of the partial correlation coefficients of all outgoing edges corresponding to each node is calculated, and the physical quantity node with the largest sum of the absolute values ​​of the partial correlation coefficients is identified as the disturbance source node.

[0075] Starting with the disturbance source node, the system searches for other reachable physical quantity nodes along the directed edges of the directed topology graph, and determines the disturbance propagation order based on the lag time difference corresponding to each directed edge. When a physical quantity node has only one outgoing edge, the corresponding nodes are connected sequentially to form a serial disturbance propagation path; when a physical quantity node has multiple outgoing edges, each outgoing edge is treated as a different parallel propagation branch; when different propagation branches continue to point to other physical quantity nodes, a multi-branch disturbance propagation path is formed.

[0076] If the occurrence of disturbance propagation is not determined based on the associated time series characteristics, the cross-physical quantity coupling analysis within the current analysis time window is not performed, the disturbance propagation path and real-time heat exchange efficiency are not updated, and the current control parameters are maintained to enter the next analysis time window to re-determine whether disturbance propagation has occurred.

[0077] It should be noted that the disturbance propagation path is used to represent the direction, sequence, and time relationship of the disturbance propagation from the disturbance source node to other physical quantity nodes, and is used to determine the gas flow rate data, cavity pressure data, and their corresponding furnace tube regions required for subsequent heat transfer efficiency calculations. Therefore, step S2 determines the disturbance propagation path through disturbance judgment, partial correlation analysis, topology connection diagram construction, disturbance source node determination, and lag time analysis.

[0078] like Figure 2 As shown in the figure, this diagram illustrates the distribution of partial correlation matrices obtained after cross-physical quantity coupling analysis of various physical quantity nodes in a semiconductor vertical furnace. A comparison clearly shows that in the coupling analysis matrix containing 8 physical quantity nodes (5 temperature zones (T1–T5), 2 gas branch flow rate nodes (F1–F2), and 1 cavity pressure node (P), the diagonal region is deep red (correlation coefficient 1.0), corresponding to the autocorrelation relationship of each node. The off-diagonal regions between temperature zones T1–T5 exhibit varying degrees of orange-red (correlation coefficient 1.0). The correlation coefficients (0.25–0.75) indicate significant thermal conduction coupling between different temperature zones. A distinct red cluster (correlation coefficient 0.50–0.75) is observed between the gas flow nodes F1–F2 and the cavity pressure node P, indicating a close dynamic correlation between gas flow and cavity pressure. After cross-correlation analysis and lag time difference calculation, the cavity pressure node P is marked as the disturbance source node (marked with a red pentagram), exhibiting the highest out-degree, indicating that changes at node P precede those at other nodes and propagate towards T1–T5 and F1–F2. Through this cross-physical quantity coupling analysis, the indirect coupling effects previously mixed in with simple correlations are effectively eliminated, and the direct correlation between each physical quantity node is clearly quantified after excluding the influence of other nodes. The direction of disturbance propagation is precisely determined by the lag time difference. This result intuitively demonstrates that the method described in this scheme for determining the disturbance propagation path based on the partial correlation matrix and topological connection graph can accurately identify the disturbance source node and its propagation direction in a dynamic coupling environment of multiple physical quantities, providing a solid data foundation for subsequent calculation of real-time heat exchange efficiency and prediction of future temperature fields based on the disturbance propagation path.

[0079] In step S3, based on the disturbance propagation path, the heat transfer efficiency is calculated using the real-time gas flow rate data and real-time cavity pressure data at the current moment to obtain the real-time heat exchange efficiency, including:

[0080] Based on the disturbance propagation path, the flow fluctuation amplitude is extracted from the real-time gas flow data, and the pressure gradient change is extracted from the real-time cavity pressure data;

[0081] Convection heat transfer characteristics are constructed based on the flow rate fluctuation amplitude and the pressure gradient change, and the heat transfer delay in the gas flow process is determined based on the convection heat transfer characteristics. A thermal resistance mapping diagram is generated based on the heat transfer delay.

[0082] If at least one region in the thermal resistance mapping map satisfies a preset thermal resistance threshold, then the dynamic heat flux density corresponding to the region is analyzed.

[0083] The real-time heat exchange efficiency is determined based on the ratio between the dynamic heat flux density and the total input heat.

[0084] In one implementation, the step of resolving the dynamic heat flux density corresponding to the region if at least one region in the thermal resistance map satisfies a preset thermal resistance threshold further includes:

[0085] If the preset thermal resistance threshold is not met in the thermal resistance mapping, the real-time heat exchange efficiency corresponding to the thermal resistance mapping is determined based on the average value of the dynamic heat flux density corresponding to each region.

[0086] It should be noted that, based on the disturbance propagation path, the gas flow channel, cavity pressure monitoring point, and corresponding furnace tube space region involved in the heat transfer efficiency calculation are determined. A sliding window algorithm is used to calculate the difference between the maximum and minimum gas flow rates in the real-time gas flow time series, obtaining the flow fluctuation amplitude. A spatial difference algorithm is used to calculate the pressure difference between adjacent pressure monitoring points, and the pressure gradient change is determined based on the pressure difference and the distance between the corresponding pressure monitoring points.

[0087] It should be noted that, since the flow rate fluctuation amplitude and pressure gradient change have different physical dimensions, they are both dimensionless using the average flow rate fluctuation amplitude and the average pressure gradient change over multiple normal historical operating cycles under the corresponding process stage. The dimensionless flow rate fluctuation amplitude and pressure gradient change are then concatenated in a fixed order to obtain the convective heat transfer characteristics. These characteristics characterize the combined impact of current gas flow velocity changes and pressure loss along the flow path on the heat transfer process inside the furnace tube. In this embodiment, for each process stage, data from at least 10 normal historical operating cycles are selected to calculate the average flow rate fluctuation amplitude μ. ΔQ and the average value of the pressure gradient change μ ΔP Dimensionless processing is achieved by dividing by the corresponding average value, i.e., ΔQ*=ΔQ / μ ΔQ ΔP*=ΔP / μ ΔP ΔQ refers to the flow rate fluctuation amplitude, which is the difference between the maximum and minimum real-time gas volumetric flow rate within a sliding time window. The dimensionless flow rate fluctuation amplitude and pressure gradient change are spliced ​​together in the order of first the flow rate fluctuation amplitude and then the pressure gradient change to form a two-dimensional convective heat transfer characteristic vector.

[0088] Based on the pressure gradient change and the spatial distance corresponding to the disturbance propagation path, the pressure loss along the disturbance propagation path is determined. Based on the real-time gas flow rate and pressure loss, the energy dissipation rate of the gas flowing along the disturbance propagation path is determined. The energy dissipation rate is then distributed according to the furnace tube space regions traversed by the disturbance propagation path, obtaining the energy dissipation distribution corresponding to each space region. Simultaneously, based on the real-time gas flow rate, gas flow cross-sectional area, and the spatial distance corresponding to the disturbance propagation path, the heat transfer delay required for the gas to flow from the region corresponding to the disturbance source node to other regions is determined.

[0089] The energy dissipation rate is determined by the product of pressure loss and real-time gas volumetric flow rate, and the calculation formula is Ẇ diss =ΔP loss ×Q, where ΔP loss The pressure loss is determined based on the spatial distance of the pressure gradient change and the disturbance propagation path, where Q is the real-time gas volumetric flow rate. The unit of energy dissipation rate is watts. The energy dissipation rate is allocated according to the proportion of pressure loss in each spatial region, i.e., the energy dissipation rate allocated to the k-th spatial region is Ẇ. diss,k =Ẇ diss ×(ΔP loss,k / Σ i ΔP loss,i ), where ΔP loss,k Let be the pressure loss corresponding to the k-th spatial region. The allocated energy dissipation distribution is used to determine the distribution ratio of the total input heat among the spatial regions, thereby determining the heat flow rate transferred to each spatial region in the subsequent thermal resistance calculation. The formula for calculating the heat transfer delay is τ. delay =L / (Q / A), where L is the spatial distance corresponding to the disturbance propagation path, Q is the real-time gas volumetric flow rate in m³ / s, and A is the gas flow cross-sectional area.

[0090] The heat transfer delay is not directly used as the thermal resistance value, but rather for time alignment of the input heat change corresponding to the disturbance source node and the temperature response corresponding to each spatial region. After time alignment, the thermal resistance value of each spatial region is determined based on the ratio between the temperature change of each spatial region at the corresponding time and the heat flow rate transferred to that spatial region. The unit of thermal resistance value is Kelvin per watt, used to represent the degree of impediment to heat transfer in the corresponding spatial region. The thermal resistance values ​​are arranged according to the position of each spatial region in the grid coordinate system inside the furnace tube to generate a thermal resistance mapping map. The heat flow rate transferred to each spatial region is q. k =P total ×(1-Ẇ diss,k / Σ i Ẇ diss,i ×α loss ) / N region , where Ptotal α represents the total output power of the heaters in each temperature zone during the current control cycle. loss N is the preset energy dissipation ratio coefficient. region This represents the number of spatial regions traversed by the disturbance propagation path. In this embodiment, α loss Set to 1×10⁻ 4 .

[0091] The thermal resistance values ​​corresponding to each spatial region in the thermal resistance mapping diagram are compared with a preset thermal resistance threshold. When the thermal resistance value of any spatial region is greater than or equal to the preset thermal resistance threshold, the corresponding spatial region is identified as the target region. In this embodiment, the preset thermal resistance threshold is determined based on the statistical results of the thermal resistance values ​​of each spatial region within a normal historical operating cycle and is set to 0.15 Kelvin per watt, so that spatial regions exceeding the normal heat transfer resistance range enter the dynamic heat flux density analysis process. In this embodiment, the 95th percentile of the thermal resistance values ​​of each spatial region within at least 50 normal historical operating cycles is used as the preset thermal resistance threshold, i.e., the thermal resistance threshold is 0.15 Kelvin per watt, so that spatial regions with thermal resistance values ​​higher than the upper limit of the normal operating range are identified as target regions that need to undergo dynamic heat flux density analysis.

[0092] Transient heat conduction analysis is performed on the target regions. Specifically, the furnace tube temperature distribution acquired at the current moment is used as the initial temperature field, the heat flow input from the heaters in each temperature zone to the outer wall of the furnace tube is used as the outer boundary condition, and the convective heat transfer between the inner wall of the furnace tube and the process gas is used as the inner boundary condition. Based on the thermal conductivity, density, and specific heat capacity of the furnace tube material, as well as the current gas temperature, cavity pressure, gas specific heat capacity, and gas thermal conductivity, a transient heat conduction model corresponding to each target region is established. Combining the grid size, region area, and temperature changes between adjacent grid nodes of each target region, the transient heat conduction model is discretized and solved to obtain the dynamic heat flux density of each target region within the current control cycle.

[0093] The governing equation for the transient heat conduction model is ρc p ∂T / ∂t=k(∂²T / ∂x²+∂²T / ∂y²), where ρ is the density of the furnace tube material, and c p Let be the specific heat capacity of the furnace tube material, T be the temperature, t be the time, and k be the thermal conductivity of the furnace tube material. The outer boundary condition is the heat flow boundary from the heaters in each temperature zone to the outer wall of the furnace tube, i.e., -k∂T / ∂n| outer =q heater The internal boundary condition is the convective heat transfer boundary between the inner wall of the furnace tube and the process gas, i.e., -k∂T / ∂n| inner =h gas (T wall -T gas ), where h gas T is the convective heat transfer coefficient of the process gas.wall is the inner wall temperature of the furnace tube, T gas is the process gas temperature. The finite difference method is used for discrete solution, the grid size is uniformly divided into 10 grid nodes along the wall thickness direction of the furnace tube, and the time step is set to 0.1 second. In this embodiment, the furnace tube material is quartz, whose thermal conductivity varies with temperature, and the corrected value at the working temperature is directly adopted, for example, about 2.14 W / (m·K) at 800°C, the density is 2200 kg / m³, and the specific heat capacity is 750 J / (kg·K). The specific heat capacity and thermal conductivity of the process gas are determined according to the type and temperature of the current process gas: for nitrogen, the specific heat capacity at 800°C is about 1120 J / (kg·K), and the thermal conductivity is about 0.057 W / (m·K).

[0094] Convective heat transfer coefficient h gas is determined according to the Nusselt number correlation formula Nu=0.023·Re 0 .8·Pr 0 .4, where Re is the Reynolds number and Pr is the Prandtl number, and h gas =Nu·k gas / D h is converted to obtain, k gas is the thermal conductivity of the process gas, D h is the hydraulic diameter of the furnace tube. For the typical process conditions of a vertical furnace (nitrogen flow rate 5 to 20 SLM, furnace tube inner diameter 200 mm to 300 mm, temperature 600°C to 1200°C), the Reynolds number is usually in the range of 10 4 to 10 5 , which meets the turbulent flow application conditions of the Dittus-Boelter correlation (Re>10000, 0.7<Pr<160). This two-dimensional transient heat conduction equation is applicable to the heat transfer analysis in the wall thickness direction (x direction) of the furnace tube and the gas flow direction (y direction). The temperature gradient in the axial direction (z direction) is smaller than that in the radial direction, and the recursive calculation of MPC already includes the independent evolution of each axial temperature node, so the local transient heat conduction analysis of each target area can be simplified to a two-dimensional model.

[0095] When there is no area meeting the preset thermal resistance threshold in the thermal resistance map, the dynamic heat flux density corresponding to each area in the thermal resistance map is determined respectively, the average value of the dynamic heat flux density corresponding to each area is calculated, and the average value is taken as the overall dynamic heat flux density inside the furnace tube.

[0096] It should be noted that, for the current control cycle, the dynamic heat flux density corresponding to each spatial region is multiplied by the effective heat transfer area of ​​the corresponding region to obtain the heat flux rate of each spatial region. The heat flux rate is then integrated over time within the current control cycle, and the integral results for each spatial region are summed to obtain the effective heat exchange within the current control cycle. The output power of the heaters in each temperature zone is then integrated over time within the current control cycle to obtain the total input heat within the current control cycle. The real-time heat exchange efficiency is determined based on the ratio between the effective heat exchange and the total input heat.

[0097] Real-time heat exchange efficiency is used to characterize the proportion of heat input from the heaters in each temperature zone that actually participates in the heat exchange inside the furnace tube under the current gas flow rate and cavity pressure conditions, and serves as the input for predicting the temperature field at future moments.

[0098] In step S4, based on the real-time heat exchange efficiency and the real-time gas flow rate data, the temperature field at future times is predicted to obtain the predicted temperature distribution, including:

[0099] Based on the real-time heat exchange efficiency and the real-time gas flow rate data, determine the state transition matrix and control input sequence corresponding to the temperature field;

[0100] Based on the state transition matrix and the control input sequence, a model predictive control algorithm is used to generate the initial temperature field evolution trajectory within a future preset time window;

[0101] The spatial temperature gradient and prediction error covariance are determined based on the initial temperature field evolution trajectory, and the optimal control increment is solved based on the spatial temperature gradient when the prediction error covariance meets the preset convergence condition.

[0102] The initial temperature field evolution trajectory is updated based on the optimal control increment to obtain the predicted temperature distribution.

[0103] Specifically, the actual temperature values ​​corresponding to each temperature monitoring node within the furnace tube at the current moment are arranged into a temperature state vector according to the spatial order of the nodes, and this temperature state vector is used as the current state of the temperature field dynamic model. Let n be the number of temperature nodes within the furnace tube used to characterize the temperature field, and let the number of heaters in each temperature zone be... The number of gas flow control channels is The temperature state vector includes n temperature state components, and the control input vector includes... Each heater power component and Each gas flow component.

[0104] In one implementation, a subspace identification algorithm is used to determine the dynamic model of the temperature field based on temperature status data, heater power data for each temperature zone, gas flow data, and corresponding real-time heat exchange efficiency over multiple historical operating cycles. The real-time heat exchange efficiency is used as a scheduling parameter of the dynamic temperature field model to correct the state transition matrix, but it is not considered an independent control component in the control input vector. In this embodiment, the subspace identification algorithm employs a numerical subspace state-space system identification algorithm, determining the order and parameters of the state-space model based on the QR decomposition and singular value decomposition of the historical input data matrix and the historical output data matrix. The historical input data matrix consists of heater power data and gas flow data at each sampling time for at least 30 historical operating cycles, and the historical output data matrix consists of temperature status data at the corresponding sampling time. The model order is determined using the Akaike information criterion, selecting the order with the smallest AIC value as the final model order.

[0105] Specifically, a reference heat exchange efficiency, a reference state transition matrix, a state transition correction matrix corresponding to changes in heat exchange efficiency, and a control input matrix are pre-determined based on historical operating data. The temperature state vector at each historical sampling moment is used as the model state, the heater power and gas flow rate at the corresponding sampling moment are used as the model control input, and the real-time heat exchange efficiency at the corresponding sampling moment is used as the scheduling variable for the state transition matrix, thus establishing a historical input data matrix and a historical output data matrix. Based on the historical input data matrix and the historical output data matrix, the dynamic relationship between temperature state and changes in heater power, gas flow rate, and real-time heat exchange efficiency is identified, resulting in the reference state transition matrix, the state transition correction matrix, and the control input matrix.

[0106] When there is a difference between the real-time heat exchange efficiency and the reference heat exchange efficiency for the current control cycle, the difference is multiplied by the state transition correction matrix, and the resulting product is added to the reference state transition matrix to obtain the state transition matrix for the current control cycle. This linear correction relationship holds true when the difference between the real-time heat exchange efficiency and the reference heat exchange efficiency is no greater than 0.3. When the difference exceeds 0.3, the difference is truncated to 0.3 before linear correction. The reference heat exchange efficiency is set as the median of the real-time heat exchange efficiency over multiple historical operating cycles.

[0107] In this embodiment of the invention, the temperature field evolution process is represented by the following discrete state-space model:

[0108]

[0109] in, Indicates the first The temperature state vector at each sampling time has a dimension of . ; This represents the state transition matrix after correction based on the current real-time heat exchange efficiency, and its dimension is... ; This represents the control input matrix, and its dimension is... ; This represents the control input vector, and its dimension is... Control the front of the input vector Each component corresponds to the heater power in each temperature zone, then Each component corresponds to the real-time gas flow rate of each gas flow control channel. This represents the output vector of the temperature field. This represents the state output matrix. When each temperature state component can be directly obtained through its corresponding temperature monitoring node, the state output matrix is ​​set to... An identity matrix of order 1, which makes the temperature field output vector correspond to the temperature state vector.

[0110] The heater power data and gas flow rate data are dimensionless according to their respective rated operating ranges to form the control input vector. Based on the control input matrix, the heater power components and gas flow rate components with different dimensions are converted into temperature change contributions to each temperature state component. The sampling period of the state-space model is set to 0.5 seconds, and the prediction time domain of the model predictive control algorithm is set to 10 sampling periods to generate the temperature field evolution trajectory within the next 5 seconds.

[0111] The prediction time domain is set to 10 sampling periods (5 seconds) to match the time constant of the furnace tube temperature response: with a 0.5-second sampling period, 5 seconds can cover the main rise phase of the temperature response, while keeping the computational load within an acceptable range for real-time solutions. The control time domain of the MPC is equal to the prediction time domain, i.e., the control time domain is also set to 10 sampling periods.

[0112] Based on the current temperature state vector, the corrected state transition matrix, and the current control input vector, the temperature state over the next 10 sampling periods is recursively calculated to obtain the initial temperature field evolution trajectory. Temperature nodes within the furnace tube are arranged in axial coordinate order, and the spatial distance between adjacent temperature nodes is determined by the installation position of the corresponding temperature sensors. For temperature nodes located inside the furnace tube, the spatial temperature gradient is determined using a central difference method based on the ratio of the temperature difference between two adjacent temperature nodes to their corresponding spatial distance. For temperature nodes located at the top and bottom boundaries of the furnace tube, forward and backward difference methods are used to determine the spatial temperature gradient, respectively. The direction of the spatial temperature gradient is determined according to the axial direction from the bottom to the top of the furnace tube.

[0113] In the temperature field prediction process, a Kalman filter algorithm is used to correct the predicted temperature state, and the prediction error covariance matrix is ​​calculated recursively. The dimension of the prediction error covariance matrix is... In this model, the diagonal elements represent the prediction error variance of each temperature node, and the off-diagonal elements represent the covariance between prediction errors of different temperature nodes. The process noise covariance and measurement noise covariance in the Kalman filter are determined based on the model prediction error and temperature sensor measurement error over multiple historical operating cycles under normal process conditions, and remain unchanged during model operation. When the trace of the prediction error covariance matrix is ​​less than 0.1, the prediction error covariance is considered to meet the preset convergence condition. In this embodiment, the process noise covariance matrix Q of the Kalman filter... W Set to 0.01×I n , where I n R is an n-order identity matrix; the measurement noise covariance matrix is ​​R. W Set to 0.04×I n When the trace of the prediction error covariance matrix is ​​less than 0.1, the prediction error covariance is deemed to meet the preset convergence condition. The threshold of 0.1 corresponds to an average prediction standard deviation of approximately 0.1 for each temperature node (after dimensionless scaling of the state variables), which is equivalent to an average prediction uncertainty of approximately 0.5°C on the actual temperature scale. The above Q... W and R W The numerical values ​​correspond to the state variables after linearization and dimensionlessness using the rated operating range of each control channel. If different standardization methods are used, Q... W and R W The absolute value needs to be converted according to the same standardized scale.

[0114] When the prediction error covariance meets the preset convergence condition, a temperature tracking error is established based on the difference between the predicted temperature and the target temperature at each temperature node within the next 10 sampling periods. The model predictive control objective function is then constructed, using the temperature tracking error and the control input increment as optimization objects. The objective function includes the weighted sum of squares of the temperature tracking error within the next 10 sampling periods, as well as the weighted sum of squares of the heater power increment and the gas flow rate increment. The weights corresponding to the temperature tracking error are used to limit the deviation of the predicted temperature from the target temperature, while the weights corresponding to the control input increment are used to suppress excessively rapid changes in heater power and gas flow rate.

[0115] The specific form of the objective function is,

[0116] J=Σ_{k=1}^{10}||y pred (k)-y ref (k)||²Q T +Σ_{k=0}^{9}||Δu(k)||²R T ;

[0117] Q T=diag(10,10,...,10) is the temperature tracking error weight matrix (an n×n diagonal matrix), R T To control the incremental weight matrix, the diagonal elements corresponding to the heater power increment are set to 1, and the diagonal elements corresponding to the gas flow increment are set to 0.5.

[0118] The optimization variables for model predictive control are the power increment of the heaters in each temperature zone and the flow increment of each gas flow control channel over the next 10 sampling periods. The optimization process satisfies the following control constraints: the power of each temperature zone heater is not less than the minimum output power of the corresponding heater and not more than the rated maximum output power; the gas flow rate of each gas flow control channel is not less than the minimum allowable flow rate of the corresponding mass flow controller and not more than the rated maximum flow rate; the power and flow rate changes between adjacent control periods do not exceed the maximum allowable rate of change of the corresponding actuator. The above power upper and lower limits, flow upper and lower limits, and rate of change limits are preset based on the equipment parameters of the heaters and mass flow controllers in the semiconductor vertical furnace.

[0119] Under control constraints, a quadratic programming algorithm is used to solve the objective function, obtaining the optimal control increment sequence corresponding to the next 10 sampling periods. The optimal control increment sequence is then added to the heater power and gas flow rate corresponding to each sampling period to obtain the updated control input sequence. The updated control input sequence is then substituted into the discrete state-space model for recursive calculation to update the initial temperature field evolution trajectory, obtaining the predicted temperature values ​​corresponding to each temperature node within the next 5 seconds. Based on the spatial location of each predicted temperature value within the furnace tube, a predicted temperature distribution is generated.

[0120] It should be noted that the predicted temperature distribution is used to represent the predicted temperature state of each temperature zone or spatial node of the semiconductor vertical furnace within a preset time range in the future, and is used for comparison with the preset target temperature distribution in step S5. When the prediction error covariance does not meet the preset convergence condition, the current state of the temperature field is updated according to the actual temperature data collected at the current moment, and the initial temperature field evolution trajectory and prediction error covariance are recalculated until the prediction error covariance meets the preset convergence condition or reaches the preset number of iterations.

[0121] like Figure 3As shown in the figure, this graph compares the temperature distribution of each temperature zone in the semiconductor vertical furnace at the current time, the predicted future time, and the target time. The comparison clearly shows that the current temperature curve (blue solid line, marked with dots) exhibits a linear upward trend within the temperature zone from zone 0 to zone 4, with the temperature increasing from approximately 799°C to approximately 820°C, generally about 5-6°C lower than the target temperature curve (green dashed line, marked with triangles). The predicted temperature curve (red solid line, marked with squares) generated by the model predictive control algorithm lies between the current and target temperatures. At zone 0, it is approximately 802°C, an increase of about 3°C ​​from the current temperature; at zone 4, it is approximately 823°C, an increase of about 3°C ​​from the current temperature, and gradually approaches the target temperature of 825°C. The difference between the predicted and target temperatures converges from approximately 3°C at zone 0 to approximately 2°C at zone 4. After the state transition matrix is ​​corrected based on real-time heat exchange efficiency and the initial temperature field evolution trajectory is generated using a model predictive control algorithm, the current temperature distribution, which was originally far from the target temperature, is effectively guided to converge towards the target temperature curve. The deviation between the predicted temperature and the target temperature in each temperature zone is significantly reduced in the prediction time domain. This result intuitively proves that the method described in this scheme for predicting the future temperature field based on real-time heat exchange efficiency and real-time gas flow data can obtain the temperature change trend and solve the optimal control increment before the temperature deviation further expands. This provides a solid data foundation for subsequently solving the control quantity in reverse based on the difference between the predicted temperature distribution and the target temperature distribution, and for coordinating the adjustment of the heater power and mass flow controller in each temperature zone.

[0122] In step S5, if the difference between the predicted temperature distribution and the preset target temperature distribution is greater than a preset temperature deviation threshold, then the control quantity used for coordinated adjustment of the heater power and mass flow controller of each temperature zone of the semiconductor vertical furnace is solved in reverse based on the difference, and coordinated adjustment is performed based on the control quantity, including:

[0123] The predicted temperature distribution is compared with the target temperature distribution node by node to obtain the node temperature difference vector, and the compensation heat flux density is determined based on the node temperature difference vector.

[0124] The compensated heat flux density is subjected to reverse heat transfer to determine the heater power adjustment amount and the flow rate adjustment amount of the mass flow controller for each temperature zone.

[0125] The power of the heaters in each temperature zone and the mass flow controller are adjusted in a coordinated manner based on the power adjustment amount and the flow rate adjustment amount.

[0126] In one implementation, if the difference between the predicted temperature distribution and the preset target temperature distribution is greater than a preset temperature deviation threshold, the method further includes:

[0127] If the difference between the predicted temperature distribution and the preset target temperature distribution is not greater than the preset temperature deviation threshold, then the current control parameters of the heater power and mass flow controller in each temperature zone remain unchanged, and the next control cycle begins.

[0128] Specifically, the target temperature value corresponding to each spatial node in the target temperature distribution is subtracted from the predicted temperature value of the corresponding spatial node in the predicted temperature distribution to obtain the node temperature difference vector. For the i-th spatial node, the difference between the target temperature value and the predicted temperature value corresponding to that spatial node is determined as the signed temperature difference component.

[0129] When the signed temperature difference component is positive, it indicates that the predicted temperature of the corresponding space node is lower than the target temperature, and the input heat needs to be increased; when the signed temperature difference component is negative, it indicates that the predicted temperature of the corresponding space node is higher than the target temperature, and the input heat needs to be reduced.

[0130] The absolute values ​​of each signed temperature difference component are compared with a temperature deviation threshold. A temperature compensation process is triggered when the absolute value of at least one signed temperature difference component in the node temperature difference vector exceeds the temperature deviation threshold. The absolute value of the signed temperature difference component is used to determine whether temperature compensation is triggered, and its sign (positive or negative) determines the direction of compensation: increasing or decreasing heat.

[0131] For example, when the target temperature of a certain temperature zone is 800 degrees Celsius, the predicted temperature is 796 degrees Celsius, and the preset temperature deviation threshold is 2 degrees Celsius, the temperature difference component corresponding to the temperature zone is 4 degrees Celsius, which is greater than the temperature deviation threshold. Therefore, the compensation heat flux density calculation is triggered, and the corresponding heater power adjustment and gas flow rate adjustment are further determined. The above values ​​are only used to illustrate this embodiment and do not constitute a limitation. In one implementation, a node temperature difference field is constructed based on the signed temperature difference component corresponding to each spatial node and the spatial coordinates of each spatial node inside the furnace tube. For spatial nodes located inside the furnace tube, the spatial gradient of the node temperature difference field is determined by the central difference method based on the ratio between the difference of the signed temperature difference components of the adjacent spatial nodes before and after the spatial node and the corresponding node spacing. For spatial nodes located at the top and bottom boundaries of the furnace tube, the spatial gradient is determined by the backward difference method and the forward difference method, respectively. The positive direction of the spatial gradient is determined according to the axial direction from the bottom of the furnace tube to the top of the furnace tube. In this embodiment, the preset temperature deviation threshold is set to 2 degrees Celsius. The threshold is determined based on the process temperature tolerance range of the semiconductor vertical furnace during the isothermal stage of the processed wafer: for conventional silicon oxidation processes, the process temperature tolerance during the isothermal stage is usually ±2°C to ±5°C. The lower limit of tolerance of 2°C is taken as the deviation threshold to achieve early prediction and timely compensation of temperature deviation.

[0132] Based on the spatial gradient of the nodal temperature difference field and the thermal conductivity of the furnace tube material, the thermal conductivity compensation component corresponding to each spatial node is determined using the Fourier thermal conductivity relation. The transient heat compensation component corresponding to each spatial node is determined based on the signed temperature difference component corresponding to each spatial node, the equivalent heat capacity of the corresponding spatial region, the effective heat transfer area, and the duration of the current control cycle. The Fourier thermal conductivity relation is q. cond =-k×∇(ΔT), where q cond denoted as the thermal conductivity compensation component (in W / m²), k = 1.4 W / (m·K) is the thermal conductivity of the quartz furnace tube material, and ∇(ΔT) is the spatial gradient of the nodal temperature difference field.

[0133] The equivalent heat capacity is determined based on the mass and specific heat capacity of the furnace tube material and process gas within the corresponding spatial region. The thermal conductivity compensation component is used to compensate for uneven temperature distribution between adjacent spatial nodes, while the transient heat compensation component is used to compensate for the heat required for the overall deviation of the corresponding spatial node from the target temperature. The thermal conductivity compensation component and the transient heat compensation component are superimposed to obtain the compensated heat flux density corresponding to the spatial node. For the k-th spatial region, its equivalent heat capacity is...

[0134] C k =(ρ wall ×V wall,k ×c p,wall +ρ gas ×V gas,k ×c p,gas ), where ρ wall =2200kg / m³ is the density of quartz, the material of the furnace tube. p,wall =750J / (kg·K) is the specific heat capacity of quartz, ρ gas c is the density of the process gas at the current temperature and pressure. p,gas V represents the specific heat capacity of the process gas. wall,k and V gas,k These represent the volumes occupied by the furnace tube material and process gas in the k-th spatial region, respectively. The thermophysical parameters of the quartz material mentioned above are typical engineering values ​​around 800°C. In practical applications, these values ​​can be adjusted according to the specific furnace tube material and operating temperature range.

[0135] It should be noted that the compensation heat flux density is used to characterize the increase or decrease in heat flux required per unit time through a unit area to reduce the temperature deviation of the corresponding spatial node, rather than directly representing the total compensation heat required to eliminate the corresponding temperature deviation.

[0136] The compensated heat flux density corresponding to each spatial node is multiplied by the effective heat transfer area of ​​the corresponding spatial region, and then integrated over time within the current control period to obtain the compensated heat demand for each spatial region within the current control period. The compensated heat demands for each spatial region are combined according to the fixed arrangement order of the spatial nodes to obtain the compensated heat demand vector, and the compensated heat demand vector is used as the target boundary condition for the inverse heat transfer model.

[0137] Each row in the thermal response matrix corresponds to a spatial node, and each column corresponds to a temperature zone heater or a gas flow control channel. Each element in the matrix represents the degree of influence of a unit change in control quantity on the heat change of the corresponding spatial node. The thermal response matrix is ​​pre-identified based on the control quantity changes of each control channel and the temperature response of the corresponding spatial node over multiple historical operating cycles, and is fixed and saved before system operation. The identification of the thermal response matrix uses the least squares method. For no less than 30 historical operating cycles, the unit step response of each control channel is recorded (heater power increases by 10% of rated power each time, and gas flow increases by 10% of rated flow each time), and the steady-state increment of temperature change at each spatial node is measured. The ratio of the steady-state increment to the corresponding control quantity change is used as the element at the corresponding position in the thermal response matrix. The thermal response matrix has n rows (m... h +m f A column matrix, where the element H in the i-th row and j-th column is... ij This represents the steady-state temperature change of the i-th spatial node when the j-th control channel undergoes a unit change in control quantity.

[0138] The compensation heat demand vector and thermal response matrix are input into the inverse heat transfer model. Under constraints of the upper and lower limits of heater power in each temperature zone, the upper and lower limits of flow rate in each gas flow control channel, the maximum power change, and the maximum flow rate change, the conjugate gradient method is used for iterative solution to obtain the target heat source distribution matrix. req Where H is the thermal response matrix, Δu is the control variable change vector to be solved, and Q req To compensate for the heat demand vector, the iterative convergence condition is that the 2-norm of the residual vector is less than 1 × 10⁻³ multiplied by Q. req The 2-norm is used, and the maximum number of iterations is set to 200. If the iterations fail to converge after reaching the maximum number of iterations, the result of the current iteration is taken as an approximate solution.

[0139] The target heat source distribution matrix is ​​used to represent the target regulation state corresponding to each temperature zone heater and each gas flow control channel to meet the compensation heat demand of each spatial node. When there are multiple solution results satisfying the conditions in the reverse heat transfer model, the solution result with the smallest weighted square sum of heater power regulation and gas flow regulation is selected under the condition that the compensation heat demand residual is not greater than the preset residual threshold, in order to reduce the amplitude of control quantity variation and determine the compensation distribution relationship between power regulation and flow regulation. The preset residual threshold is consistent with the convergence condition of the conjugate gradient method, that is, the 2-norm of the residual vector is less than 1×10⁻³ multiplied by Q. req The 2-norm of the equation. In the weighted sum of squares, the weight corresponding to the heater power regulation is set to 1, and the weight corresponding to the gas flow regulation is set to 0.2, so that gas flow regulation is used preferentially over heater power regulation under the same compensation capability.

[0140] Before PCA dimensionality reduction, the columns of the target heat source distribution matrix are centered. The number of principal components after dimensionality reduction is determined based on the principle that the cumulative variance contribution rate reaches 95%. When the original column dimension of the target heat source distribution matrix (i.e., the number of control channels m) is... h +m f If the cumulative variance contribution rate of the first two principal components is found to be above 95% by PCA analysis when the value is not greater than 10, then the first two principal components are retained as driving adjustment signals.

[0141] A pre-constructed power-flow dynamic coupling matrix is ​​invoked, and the power regulation component and flow regulation component are input into the power-flow dynamic coupling matrix. The power-flow dynamic coupling matrix is ​​pre-identified based on the dynamic relationship between heater power changes, gas flow rate changes, and corresponding temperature responses over multiple historical operating cycles. The identification method for the power-flow dynamic coupling matrix is ​​as follows: during multiple historical operating cycles, independent pseudo-random binary sequence excitations are applied to each heater channel and each gas flow channel, the response data of each control channel is recorded, and the steady-state gain and dynamic coupling coefficient between each control channel are identified using a multivariate output error least squares algorithm. The power-flow dynamic coupling matrix is ​​(m h +m f ) order square matrix.

[0142] In the power-flow dynamic coupling matrix, the diagonal elements represent the direct impact of heater power or gas flow rate changes on the corresponding controlled object, while the off-diagonal elements represent the cross-influence between heater power and gas flow rate changes. Singular value decomposition is performed on the power-flow dynamic coupling matrix to obtain mutually decoupled control components. These control components are then converted into control duty cycles for each temperature zone heater and valve opening signals for each mass flow controller. For (m... h +m fSingular value decomposition of the power-flow dynamic coupling matrix M of order 1: M = U·Σ·V T U and V T Let V be the left and right singular vector matrices, respectively, and let Σ be a diagonal matrix arranged in descending order of singular values. The decoupled control components are obtained by right-multiplying the drive adjustment signal by V. T Get: u decoupled =V T ·u driven , where u driven The driving adjustment signal is obtained after PCA dimensionality reduction. The decoupled control component u... decoupled The first m h Each component is proportionally mapped to the control duty cycle (0% to 100%) of the heater in each temperature zone, and then the m f Each component is proportionally mapped to the valve opening signal (0% to 100%) of each mass flow controller. Wherein, u decoupled For each decoupled control component that corresponds one-to-one with a control channel, its dimension is (m h +m f ).

[0143] When determining the control duty cycle and valve opening signal, they are compared with the minimum allowable value, maximum allowable value, and maximum allowable change within a unit control cycle of the corresponding actuator. When the calculated control signal exceeds the allowable range of the corresponding actuator, it is limited to the corresponding boundary value; when the change of the calculated control signal between adjacent control cycles exceeds the maximum allowable change of the corresponding actuator, it is truncated according to the maximum allowable change.

[0144] When any control channel reaches its upper or lower allowable limit, that control channel is fixed at the corresponding boundary value, and the remaining compensation heat demand is recalculated using other control channels that have not reached their boundary values. Based on the processed control duty cycle and valve opening signals, the heater power and mass flow controller of each temperature zone are synchronously adjusted to achieve coordinated control of the temperature compensation demand of multiple nodes. The remaining compensation heat demand is obtained from the original compensation heat demand vector Q. req The compensation contribution provided by the control channels that have reached the boundary and are fixed at their boundary values ​​is determined by subtracting the compensation contribution from the boundary value. Specifically, let B be the set of indices of the control channels that have reached the boundary, then the remaining compensation heat demand is Q. residual =Q req -Σ_{j∈B}h j· Δu bound,j , where h j Let Δu be the j-th column of the thermal response matrix H. bound,j This represents the change in control quantity when the j-th control channel is fixed at the boundary value. The conjugate gradient method is used to simplify the dimension-reduced equation H. reduced· Δuremaining =Q residual Resolve, where H reduced The result is obtained by removing columns from index set B from H. If control channels still reach boundary values ​​after resolving, the above boundary fixing and resolving process is repeated, up to a maximum of 3 times; if the convergence condition is still not met after 3 times, the current solution result is used as the final control output.

[0145] When the absolute value of each signed temperature difference component in the node temperature difference vector is not greater than the temperature deviation threshold, the current control parameters of the heater power and mass flow controller in each temperature zone remain unchanged, and the next control cycle begins.

[0146] It should be noted that keeping the current control parameters unchanged does not mean stopping the control of the semiconductor vertical furnace. The heaters in each temperature zone and the mass flow controller will still operate according to the current control parameters, only without generating additional power and flow regulation in the current control cycle.

[0147] It should be noted that the node temperature difference vector is used to characterize the degree of deviation of the predicted temperature of each spatial node from the target temperature; the compensation heat flux density is used to characterize the increase or decrease of heat flux required per unit time through a unit area to reduce the corresponding temperature deviation; the power regulation is used to adjust the output power of the heaters in each temperature zone; and the flow rate regulation is used to adjust the flow rate of the process gas entering the furnace tube.

[0148] In summary, this invention provides a multi-process software control method for semiconductor vertical furnaces, which can solve the technical problem of insufficient temperature control accuracy in vertical furnaces due to the difficulty in identifying the dynamic coupling relationship between temperature, gas flow rate and cavity pressure.

[0149] The second embodiment of the present invention provides a multi-process software control system for a semiconductor vertical furnace, including: a memory, a processor, and a computer program stored in the memory and executable on the processor. When the processor executes the program, it implements the steps of the method described above.

[0150] It should be noted that the semiconductor vertical furnace multi-process software control system provided in this embodiment of the invention is used to execute all process steps of the semiconductor vertical furnace multi-process software control method in the above embodiment. The working principles and beneficial effects of the two are one-to-one, so they will not be described again.

[0151] The specific embodiments described above further illustrate the purpose, technical solution, and beneficial effects of the present invention. It should be understood that the above descriptions are merely specific embodiments of the present invention and are not intended to limit the scope of protection of the present invention. In particular, it should be noted that any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of the present invention should be included within the scope of protection of the present invention for those skilled in the art.

Claims

1. A multi-process software control method for a semiconductor vertical furnace, characterized in that, Executed by a computer, including: Historical temperature data, historical gas flow data, and historical cavity pressure data of each temperature zone of the semiconductor vertical furnace during its historical operating cycle are obtained, and long-term and short-term time series analysis is performed to extract correlation time series features to characterize the dynamic correlation of multiple physical quantities. When it is determined that disturbance propagation has occurred based on the associated time series characteristics, cross-physical quantity coupling analysis is performed on the associated time series characteristics to determine the disturbance propagation path; Based on the disturbance propagation path, the heat transfer efficiency is calculated using the real-time gas flow rate data and real-time cavity pressure data at the current moment to obtain the real-time heat exchange efficiency. Based on the real-time heat exchange efficiency and the real-time gas flow rate data, the temperature field at future times is predicted, and the predicted temperature distribution is obtained. If the difference between the predicted temperature distribution and the preset target temperature distribution is greater than the preset temperature deviation threshold, then the control quantity used to coordinate the adjustment of the heater power and mass flow controller of each temperature zone of the semiconductor vertical furnace is solved in reverse based on the difference, and coordinated adjustment is performed based on the control quantity.

2. The multi-process software control method for semiconductor vertical furnaces according to claim 1, characterized in that, The extraction of correlation time-series features used to characterize the dynamic correlation of multiple physical quantities includes: The historical temperature data, historical gas flow data, and historical cavity pressure data are standardized, and the standardized data are segmented according to a preset sliding time window to obtain a multidimensional time series data sequence. The multidimensional time-series data sequence is input into a preset long short-term memory network to generate the associated time-series features.

3. The multi-process software control method for semiconductor vertical furnaces according to claim 1, characterized in that, When it is determined that disturbance propagation has occurred based on the associated time-series characteristics, the step of performing cross-physical quantity coupling analysis on the associated time-series characteristics to determine the disturbance propagation path includes: A cross-physical quantity coupling analysis was performed on the associated time-series features to obtain the partial correlation matrix; Extract elements whose absolute values ​​of partial correlation coefficients are greater than a preset correlation threshold from the partial correlation matrix, and construct a topological connection graph of physical quantity interaction; The out-degree of each physical quantity node in the topology connection graph is counted, and the physical quantity node with the largest out-degree is determined as the disturbance source node; The disturbance propagation path is determined based on the time lag between the disturbance source node and other physical quantity nodes, as well as the topology connection diagram.

4. The multi-process software control method for semiconductor vertical furnaces according to claim 3, characterized in that, Extracting elements from the partial correlation matrix whose absolute value of the partial correlation coefficient is greater than a preset correlation threshold further includes: If the absolute value of the partial correlation coefficient of an element in the partial correlation matrix is ​​not greater than a preset correlation threshold, then the partial correlation coefficient is set to zero to obtain a sparse partial correlation matrix. Based on the non-zero elements in the sparse partial correlation matrix, the topological connection relationship between each physical quantity node is determined, and the topological connection graph is generated based on the topological connection relationship.

5. The multi-process software control method for semiconductor vertical furnaces according to claim 1, characterized in that, The step of calculating the heat transfer efficiency based on the real-time gas flow rate data and real-time cavity pressure data at the current moment according to the disturbance propagation path, to obtain the real-time heat exchange efficiency, includes: Based on the disturbance propagation path, the flow fluctuation amplitude is extracted from the real-time gas flow data, and the pressure gradient change is extracted from the real-time cavity pressure data; Convection heat transfer characteristics are constructed based on the flow rate fluctuation amplitude and the pressure gradient change, and the heat transfer delay in the gas flow process is determined based on the convection heat transfer characteristics. A thermal resistance mapping diagram is generated based on the heat transfer delay. If at least one region in the thermal resistance mapping map satisfies a preset thermal resistance threshold, then the dynamic heat flux density corresponding to the region is analyzed. The real-time heat exchange efficiency is determined based on the ratio between the dynamic heat flux density and the total input heat.

6. The multi-process software control method for semiconductor vertical furnaces according to claim 5, characterized in that, If at least one region in the thermal resistance mapping map satisfies a preset thermal resistance threshold, then the dynamic heat flux density corresponding to the analyzed region further includes: If the preset thermal resistance threshold is not met in the thermal resistance mapping, the real-time heat exchange efficiency corresponding to the thermal resistance mapping is determined based on the average value of the dynamic heat flux density corresponding to each region.

7. The multi-process software control method for semiconductor vertical furnaces according to claim 1, characterized in that, The step of predicting the temperature field at future times based on the real-time heat exchange efficiency and the real-time gas flow rate data, to obtain the predicted temperature distribution, includes: Based on the real-time heat exchange efficiency and the real-time gas flow rate data, determine the state transition matrix and control input sequence corresponding to the temperature field; Based on the state transition matrix and the control input sequence, a model predictive control algorithm is used to generate the initial temperature field evolution trajectory within a future preset time window; The spatial temperature gradient and prediction error covariance are determined based on the initial temperature field evolution trajectory, and the optimal control increment is solved based on the spatial temperature gradient when the prediction error covariance meets the preset convergence condition. The initial temperature field evolution trajectory is updated based on the optimal control increment to obtain the predicted temperature distribution.

8. The multi-process software control method for semiconductor vertical furnaces according to claim 1, characterized in that, If the difference between the predicted temperature distribution and the preset target temperature distribution is greater than a preset temperature deviation threshold, then the control quantity used for coordinated adjustment of the heater power and mass flow controller of each temperature zone of the semiconductor vertical furnace is solved in reverse based on the difference, and coordinated adjustment is performed based on the control quantity, including: The predicted temperature distribution is compared with the target temperature distribution node by node to obtain the node temperature difference vector, and the compensation heat flux density is determined based on the node temperature difference vector. The compensated heat flux density is subjected to reverse heat transfer to determine the heater power adjustment amount and the flow rate adjustment amount of the mass flow controller for each temperature zone. The power of the heaters in each temperature zone and the mass flow controller are adjusted in a coordinated manner based on the power adjustment amount and the flow rate adjustment amount.

9. The multi-process software control method for semiconductor vertical furnaces according to claim 8, characterized in that, If the difference between the predicted temperature distribution and the preset target temperature distribution is greater than a preset temperature deviation threshold, the method further includes: If the difference between the predicted temperature distribution and the preset target temperature distribution is not greater than the preset temperature deviation threshold, then the current control parameters of the heater power and mass flow controller in each temperature zone remain unchanged, and the next control cycle begins.

10. A multi-process software control system for a semiconductor vertical furnace, comprising a memory, a processor, and a computer program stored in the memory and executable on the processor, characterized in that, When the processor executes the program, it implements the multi-process software control method for semiconductor vertical furnaces as described in any one of claims 1 to 9.

Images