High-frequency wafer test probe card multi-temperature-zone environment self-adaptive adjustment method and control system

Abstract

The application discloses a high-frequency wafer test probe card multi-temperature zone environment adaptive adjustment method, and belongs to the technical field of semiconductor wafer testing. Firstly, a three-dimensional entity space domain is divided into finite element units, a global heat conduction matrix and a global heat capacity matrix are established, and control input and environmental disturbance are explicitly separated, and a nominal three-dimensional thermal model is obtained through time discretization; then, a degradation sensitive channel sensitive to temperature response is identified, a heater efficiency degradation factor and a contact thermal resistance drift coefficient are implanted as latent parameters to form an evolutionary three-dimensional thermal model, and the latent parameters are updated online by using ensemble Kalman filtering; a safety control strategy network generates a candidate action vector, which is verified based on an adaptive safety barrier based on a control Lyapunov function, a model prediction error drives an adaptive tightening of a safety criterion, and a rejected action is corrected in a gradient direction until the criterion is met.

Description

Technical Field

[0001] This invention relates to the field of semiconductor wafer testing technology, and more specifically, to a method and control system for adaptive adjustment of multi-temperature zone environment of high-frequency wafer testing probe cards. Background Technology

[0002] High-frequency wafer test probe cards are key components in semiconductor wafer-level testing, used to accurately transmit electrical signals from the tester to the pads of individual chips on the wafer. To meet the performance characterization requirements of chips under different temperature conditions, the chuck carrying the wafer is typically divided into multiple independent and controlled temperature zones. Each temperature zone integrates heating and cooling elements. By actively adjusting the heat flow injection in each temperature zone, a specified temperature spatial distribution is established and maintained on the wafer surface. This thermal system involves various heterogeneous materials, including the probe card ceramic substrate, the chuck body, and the multi-temperature zone heating and cooling elements. These components form a complex three-dimensional heat conduction channel through mounting surfaces, thermally conductive interface materials, and probe contact points. The evolution of its temperature field is described by unsteady-state heat conduction partial differential equations and is simultaneously affected by multiple factors such as ambient airflow, thermal cross-linking of adjacent temperature zones, and the self-heating of the components. During long-term continuous testing, the resistance wire of the heating element gradually ages, leading to a decrease in heating efficiency. The probe wears down during repeated contact, causing changes in contact force, resulting in a slow drift in the contact thermal resistance between the probe and the pad. These latent degradation processes continuously alter the dynamic response characteristics and static gain of the entire thermal system, causing the control quality of control strategies based on factory settings to gradually decline over time.

[0003] To address the inherent strong coupling, unsteady-state, and time-varying degradation characteristics of the aforementioned thermal systems, various multi-temperature zone temperature control methods have been proposed in the industry. One mainstream approach is based on proportional-integral-derivative control combined with static or dynamic decoupling networks. This approach weakens the thermal cross-linking effects between temperature zones by designing independent controllers for each temperature zone and introducing feedforward compensation. Building upon this, to further improve dynamic response performance and disturbance rejection capabilities, some solutions employ model predictive control based on lumped-parameter thermal models. This simplifies the chuck into several isothermal blocks, identifies the low-order state-space model offline or online, and generates power commands by solving a constrained finite-time optimization problem in each control cycle. Furthermore, to endow the control system with adaptability to slow degradation, existing technologies have introduced online parameter estimation methods. These methods utilize algorithms such as recursive least squares and extended Kalman filtering to adjust a few characteristic parameters in the simplified model in real time, thereby updating the internal model of the controller and attempting to make the temperature control strategy follow the changes in the characteristics of the physical system.

[0004] However, existing multi-temperature zone adaptive control methods are still not perfect in terms of degradation perception accuracy, model fidelity, and safety assurance mechanisms. First, online estimation methods based on lumped parameter simplified models, because they abstract the entire chuck into a few uniform temperature zones and ignore the spatially continuous heat flow paths represented by multi-material interfaces and hundreds of nodes, lack sufficient spatial resolution for degradation effects occurring in local areas such as probe contact areas and thermal interface layers. This makes it difficult to effectively identify and locate key degradation sources, and the parameter estimation results may deviate significantly from the actual degradation state of the physical entity. Second, current safety constraint strategies mostly adopt fixed thresholds or fixed margins, such as temperature upper limits and power saturation limits. These static constraints do not establish a quantitative correlation with the real-time fluctuations in model prediction accuracy. When the model prediction error increases due to unmodeled degradation, the constraints may be too conservative, limiting the regulation performance of the temperature control system, or too lenient, failing to effectively prevent potentially dangerous actions generated by inaccurate models. Furthermore, when optimized control actions are rejected due to violations of safety constraints, existing methods typically only limit or discard the actions, failing to utilize the physical gradient information of the system's thermal dynamics to automatically generate alternative actions that balance safety constraints and temperature control performance. This results in limitations in maintaining the system's temperature uniformity under degradation conditions. These shortcomings mean that the reliability of existing methods in maintaining high-precision and uniform temperature control of the wafer surface throughout the probe card's entire lifecycle still has room for improvement.

[0005] In view of this, we propose a multi-temperature zone environmental adaptive adjustment method and control system for high-frequency wafer test probe cards. Summary of the Invention

[0006] The purpose of this invention is to solve the problems in the prior art, such as insufficient spatial resolution of degraded perception, disconnect between safety constraints and model confidence, and lack of physical guidance and correction for dangerous actions.

[0007] To achieve the above objectives, the present invention provides a method for adaptive adjustment of multi-temperature zone environment of a high-frequency wafer test probe card, comprising the following:

[0008] Step S1: Construct a nominal three-dimensional thermal model: Divide the three-dimensional solid space domain into multiple elements, establish the global thermal conductivity matrix and global thermal capacity matrix of each element, explicitly separate the control input vector and environmental disturbance vector from them and discretize them over time to obtain a nominal three-dimensional thermal model including multiple heat conduction channels.

[0009] Step S2, Implanting latent parameters and estimating online: Identify the degradation sensitive channels in the nominal three-dimensional thermal model, implant heater efficiency degradation factor and contact thermal resistance drift coefficient as latent parameters into the degradation sensitive channels, form an evolutionary three-dimensional thermal model containing latent parameters, and use ensemble Kalman filtering to update the latent parameter set online based on the measured temperature values ​​of multiple temperature zones;

[0010] Step S3: Train a safety control strategy network to map the observable temperature state output by the evolutionary three-dimensional thermal model into candidate action vectors. Verify the candidate action vectors using an adaptive safety barrier. Only when the candidate action vector satisfies the adaptive stability criterion based on the control Lyapunov function will it be used as a control input for execution. Otherwise, the candidate action vector will be corrected along the negative direction of the gradient of the control Lyapunov function until the adaptive stability criterion is satisfied.

[0011] Compared with the prior art, the beneficial effects of the present invention are as follows:

[0012] This invention constructs a nominal three-dimensional thermal model through finite element analysis, fully preserving multi-material interfaces, complex geometric boundaries, and spatially continuous heat conduction paths, providing a high-fidelity physical basis for degradation sensing. Based on this, degradation-sensitive channels are identified, and heater efficiency degradation factors and contact thermal resistance drift coefficients are embedded into these channels as latent parameters, forming a latent parameter vector. This vector ultimately constructs an evolutionary three-dimensional thermal model, and ensemble Kalman filtering is used to update the latent parameter set online. Since these latent parameters are directly embedded in the element thermal conductivity integral and input matrix scaling stages, each parameter update synchronously corrects the system matrix and input matrix of the thermal model, allowing the model's thermal response characteristics to follow the actual degradation trajectory of the physical system in real time. This overcomes the limitation of lumped parameter models in being unable to locate local degradation, further ensuring high spatial resolution of degradation sensing and physical consistency of parameter estimation.

[0013] Furthermore, the aforementioned latent parameter vector, combined with the global node temperature vector, serves as an augmented state. Using ensemble Kalman filtering, the estimation uncertainty is directly characterized by the statistical dispersion of the posterior set of latent parameters, eliminating the need for linear approximation of the high-dimensional system. This makes it highly suitable for nonlinear estimation scenarios where latent parameters are strongly coupled with the global node temperature vector and the thermal system matrix changes with the degradation state. The process involves iterative recursion within each control cycle, predicting, observing, and updating the latent parameter set based on measured temperatures across multiple temperature zones. Moreover, the latent parameters are directly embedded in the integral calculation of the unit thermal conductivity matrix and the scaling of the input matrix. Therefore, each update of the ensemble Kalman filter re-corrects the system matrix and input matrix of the evolving three-dimensional thermal model based on the posterior latent parameters, further enabling the model's thermal response to follow the actual degradation trajectory of heater efficiency decay and contact thermal resistance drift in the physical system online.

[0014] When a candidate action is rejected by a safety barrier, this invention calculates the gradient of the control Lyapunov function for the candidate action based on an evolutionary three-dimensional thermal model using a chain rule. The candidate action is then corrected by stepping along the negative gradient direction, with the correction step size adaptively determined by the maximum prediction deviation mapped from the current model's prediction error. Specifically, this gradient guidance provides a physically based optimization direction for the rejected action, ensuring that alternative control quantities that balance safety constraints and temperature control performance can be automatically generated under any deterioration condition.

[0015] In addition to the objectives, features, and advantages described above, the present invention has other objectives, features, and advantages. The invention will now be described in further detail with reference to the figures. Attached Figure Description

[0016] Figure 1 This is a flowchart illustrating the overall process of the present invention. Detailed Implementation

[0017] The technical solutions in 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.

[0018] This embodiment uses a four-temperature zone chuck for an 8-inch wafer probe station as an example. The probe chuck's ceramic substrate is made of aluminum nitride, and the chuck itself is made of copper alloy. Each temperature zone is independently equipped with a thin-film heating element and a thermoelectric cooler, with rated power of 150W, 120W, 150W, and 120W respectively. The temperature control objective is to establish a uniform temperature field on the wafer surface that can be set within the range of 25°C to 150°C, with a temperature uniformity requirement of ±0.5°C. During long-term continuous testing, the heating element resistance wire gradually ages due to high-temperature oxidation. The contact force between the probe tip and the wafer pad decreases with probe wear, and the contact thermal resistance drifts accordingly. This causes the control parameters calibrated at the factory to become inapplicable after several months. (Refer to...) Figure 1 The high-frequency wafer test probe card multi-temperature zone environment adaptive adjustment method shown includes the following steps:

[0019] Step S1: Construct a nominal three-dimensional thermal model: Divide the three-dimensional solid space domain into multiple elements, establish the global thermal conductivity matrix and global thermal capacity matrix of each element, explicitly separate the control input vector and environmental disturbance vector from them and discretize them over time to obtain a nominal three-dimensional thermal model including multiple heat conduction channels.

[0020] Step S2, Implanting latent parameters and estimating online: Identify the degradation sensitive channels in the nominal three-dimensional thermal model, implant the heater efficiency degradation factor and contact thermal resistance drift coefficient as latent parameters into the degradation sensitive channels, form an evolutionary three-dimensional thermal model containing latent parameters, and use ensemble Kalman filtering to update the latent parameter set online based on the measured temperature values ​​of multiple temperature zones;

[0021] Step S3: Train a safety control strategy network to map the observable temperature state output by the evolutionary three-dimensional thermal model into candidate action vectors. Use an adaptive safety barrier to verify the candidate action vectors. Only when the candidate action vector satisfies the adaptive stability criterion based on the control Lyapunov function will it be used as a control input for execution. Otherwise, the candidate action vector will be corrected along the negative direction of the gradient of the control Lyapunov function until the adaptive stability criterion is satisfied.

[0022] Steps S1, S2, and S3 are not executed in isolation: the global thermal conductivity matrix assembled in step S1 provides a quantitative benchmark for the thermal conductivity value of each heat transfer path in step S2. Based on this, step S2 calculates the partial derivatives of the steady-state temperature response of the wafer surface temperature measurement nodes under thermal conductivity perturbation for each path to screen for degradation-sensitive channels; in step S2, after each posterior update of the latent parameters by the ensemble Kalman filter, the system matrix and input matrix of the evolving three-dimensional thermal model are immediately recalculated with the updated latent parameters, thereby changing the thermal response characteristics of the model itself, so that the temperature prediction value of the next cycle is directly driven by the estimation result of this round; the adaptive stability criterion in step S3 is rigorous The severity is determined by the predicted residual sequence generated by the three-dimensional thermal model in step S2 combined with the mean of the current latent parameters through sliding window statistics. The statistical value of the residual sequence is converted into the level of the safety threshold in real time. When the candidate action vector is rejected by the safety barrier, step S3 uses the three-dimensional thermal model corresponding to the mean of the latent parameters in step S2 to calculate the gradient of the control Lyapunov function on the candidate action vector through the chain rule. The gradient direction and the correction step are determined by the maximum prediction deviation mapped by the current degenerate state of the model and the prediction error sequence, respectively. The corrected action is re-input into the three-dimensional thermal model to deduce the energy change and is compared with the current adaptive stability criterion.

[0023] Furthermore, step S2's online estimation of latent parameters enables the model to continuously track the degradation trajectory of the entity system; step S3's adaptive safety barrier quantifies the model's current prediction confidence into the stringency of the safety criterion, automatically tightening the protection threshold when the model is inaccurate and releasing more optimization space when the model is reliable; the action correction mechanism based on controlling the gradient of the Lyapunov function utilizes the model's current degradation state information to provide physically based alternative directions for rejected actions. Based on the above synergistic operation, degradation perception provides a model basis for risk assessment, risk assessment provides a quantitative basis for action correction, and action correction ensures safety under the constraints of the degraded model, ultimately achieving adaptive and precise maintenance of wafer temperature uniformity throughout the probe card system's entire lifecycle. The working principles of each of the above steps are detailed below:

[0024] In this embodiment, the probe card ceramic substrate, multi-temperature zone heating and cooling elements, and the carrier chuck in the temperature-controlled chuck occupy a three-dimensional solid space domain after assembly. Because this space domain contains contact interfaces between various heterogeneous materials, embedded assembly boundaries between the heating and cooling elements of each temperature zone and the chuck body, and convective heat transfer boundaries between the chuck's outer surface and the cleanroom air, the spatial distribution of these boundary conditions is extremely uneven and geometrically complex. Therefore, the temperature field within the space domain is spatially continuous, and heat conduction follows physical laws described by partial differential equations, which cannot be directly solved numerically. Taking a 4-temperature zone chuck of an 8-inch wafer probe station as an example: at the interface between the third and fourth temperature zones, three materials converge in the same tiny area: the chuck body is a copper alloy, the fourth temperature zone heating element is an aluminum nitride ceramic package, and a very thin layer of thermally conductive silicone grease fills the space between the heating element and the chuck body. The thermal conductivity of these three materials differs significantly. When the heating element in the fourth temperature zone is working, heat originates from the heating element, passes through the thermal grease layer, enters the copper alloy chuck body, and then diffuses along the chuck plane across the temperature zone boundary into the third temperature zone. The partial differential equations governing this process are mathematically complete, meaning that the thermal conductivity in the equations is a function of spatial position. However, due to the irregular geometry of the three materials, the non-simple plane of the material interfaces, and the real-time changes in boundary conditions with heating power and ambient airflow, the partial differential equations constitute a complex boundary value problem with spatially distributed parameters, and therefore lack an analytical solution.

[0025] To solve the above problem, step S1 divides the three-dimensional solid space domain into a finite number of elements, and each element contains multiple finite element nodes for approximating the continuous temperature field distribution inside the element through shape function interpolation, thereby defining the spatial distribution of the continuous temperature field with a finite number of discrete finite element node temperature values. Specifically:

[0026] All boundary surfaces of the 3D solid space domain are read, and boundary finite element nodes are placed on the boundary surfaces according to preset element size parameters. When placing boundary finite element nodes, the coordinates of the corresponding boundary finite element nodes fall on the boundary surfaces to ensure that the final discrete mesh can accurately reproduce the original geometry. After the boundary finite element nodes are placed, new finite element nodes are generated layer by layer from the boundary until the entire 3D solid space domain is filled. All finite element nodes are connected into independent tetrahedral elements based on the Delaunay triangulation criterion. Each element after triangulation is a sub-region with a simple geometry. Elements are interconnected through shared finite element nodes, collectively forming the original geometry of the components in the 3D solid space domain. Each element defines material thermal properties consistent with the physical component to which it belongs. For example, elements located in the probe card ceramic substrate region are assigned the thermal conductivity and specific heat capacity of the ceramic material; elements located in the multi-temperature zone heating and cooling element region are assigned the thermal properties of the corresponding functional material; and elements located in the bearing chuck region are assigned the thermal properties of metal alloys or composite ceramics.

[0027] The continuous temperature field corresponding to each element (referring to the temperature as a continuous function of spatial coordinates within the spatial subdomain occupied by the corresponding element, which is defined at infinitely many spatial points, and the temperature value changes continuously with any tiny change in spatial position, without jumps or discontinuities, thus having infinite degrees of freedom, and cannot be directly stored and solved in a computer) is approximately expressed by the temperature values ​​of each finite element node of the element through shape function interpolation. The shape function satisfies that it takes a value of 1 at its own finite element node and a value of 0 at the other finite element nodes.

[0028] Since the three-dimensional solid space domain is divided into a finite number of elements, the heat flow that originally propagated freely in any direction in the continuous medium is constrained to be transferred only along the element boundaries. The only heat exchange channel between adjacent elements is their shared finite element nodes. Therefore, heat is transferred between adjacent elements through these shared finite element nodes, and a heat conduction channel is established based on this. Specifically, the evolution of the temperature field inside the three-dimensional solid space domain satisfies the three-dimensional unsteady-state heat conduction partial differential equation. The Galerkin weighted residual method is used to obtain an equivalent weak integral form. The shape function interpolation of the temperature field is approximated into the weak form, and integral calculations are performed for each element. For any element, the heat conduction relationship between its internal finite element nodes is described by the element thermal conductivity matrix. The element values ​​in the element thermal conductivity matrix are determined by the integral of the dot product of the element's material thermal conductivity and the shape function gradient over the element's subdomain. The specific element values ​​comprehensively reflect the thermal conductivity of the element material, the geometric positional relationship between the finite element nodes, and the influence of the element volume. Meanwhile, the heat storage coupling relationship between each finite element node within the unit due to temperature changes is described by the unit heat capacity matrix. The element values ​​in the aforementioned unit heat capacity matrix are determined by the integral of the product of the unit's material density, specific heat capacity, and the shape functions of the two finite element nodes over the unit's subdomain. The values ​​comprehensively reflect the unit material's ability to store heat and the degree of spatial coupling of temperature changes between finite element nodes.

[0029] To establish a complete heat conduction path and thermal capacity coupling relationship across element boundaries, the element thermal conductivity matrix and element thermal capacity matrix of all elements are assembled into a global thermal conductivity matrix and a global thermal capacity matrix according to the finite element node numbering, respectively. The assembly rules for the global thermal conductivity matrix and the global thermal capacity matrix are the same, specifically:

[0030] For each element, all finite element node pairs within it are traversed. The corresponding element values ​​in the element's thermal conductivity matrix are accumulated to the corresponding row and column positions in the global thermal conductivity matrix, and the corresponding element values ​​in the element's thermal capacity matrix are also accumulated to the corresponding row and column positions in the global thermal capacity matrix. When two adjacent elements share a finite element node, each element contributes to the matrix elements related to the shared finite element node. These contributions are accumulated to the same position in both the global thermal conductivity matrix and the global thermal capacity matrix. In the global thermal conductivity matrix, the specific accumulation process is a natural merging of the heat flow contributions at the shared finite element node. This means that when the temperature of the shared node changes, all adjacent elements will transfer heat to or absorb heat from the node through their respective heat conduction channels. All heat flows converge at the shared node, jointly determining the net heat flow of that node.

[0031] After assembly, each non-zero off-diagonal element in the global thermal conductivity matrix corresponds to a heat conduction channel connecting two finite element nodes, and each non-zero off-diagonal element in the global thermal capacity matrix corresponds to the heat storage coupling relationship between a pair of finite element nodes. The set of heat conduction channels formed between all adjacent finite element nodes constitutes all the heat conduction channels in the nominal three-dimensional thermal model.

[0032] Based on this, the result of multiplying the global thermal conductivity matrix by the global node temperature vector describes the total heat flow transmitted in parallel through all heat conduction channels under the current temperature field distribution; the result of multiplying the global thermal capacity matrix by the time derivative of the global node temperature vector describes the total rate at which each finite element node stores or releases heat due to the material's thermal capacity effect under the current temperature change rate; the above two types of heat flow contributions, together with the heat load vector including internal heat sources, boundary heat flow, and environmental convection and radiation contributions, satisfy the instantaneous power balance relationship required by the first law of thermodynamics. From this, a set of ordinary differential equations describing the unsteady heat conduction process in the entire three-dimensional solid space domain is naturally derived, namely: the product of the global thermal capacity matrix and the time derivative of the global node temperature vector plus the product of the global thermal conductivity matrix and the global node temperature vector equals the heat load vector.

[0033] The balance between the net heat flow and the external heat load at any time in the three-dimensional solid space domain determines the dynamic trajectory of the temperature vector of the whole domain node as it evolves over time. The role of the above set of ordinary differential equations is to transform the heat conduction problem of the continuous body composed of the probe card ceramic substrate, multi-temperature zone heating and cooling elements and the carrier chuck into numerically solvable finite degree of freedom ordinary differential equations.

[0034] In the aforementioned set of ordinary differential equations, the heat load vector is a composite term representing multiple physical contributions, including internal heat sources, boundary heat flow, and environmental convection and radiation. The heat flow actively applied by heating and cooling elements in each temperature zone is mixed with external conditions such as temperature fluctuations in the cleanroom environment, making it impossible to directly distinguish which parts of the temperature field evolution are actively driven by the controller and which are passively caused by the external environment. Given the fundamental differences between control input and environmental disturbance in their impact on temperature field evolution—control input being the decision variable to be optimized, and environmental disturbance being an uncontrollable external condition—this example further explicitly separates the control input and environmental disturbance. This allows the subsequent safety control strategy network to directly use the control input vector as the optimization object, while incorporating the environmental disturbance vector as a known external condition into the state evolution equation. The explicit separation is as follows:

[0035] Based on the geometric assembly positions of the heating and cooling elements in each temperature zone within the bearing chuck, the set of finite element nodes covered by each temperature zone is determined. Assuming that the heat flow is uniformly distributed on its mating surface, the proportion of the mating surface area represented by each node is used as the heat flow distribution coefficient to construct a control input matrix. Multiplying the control input vector with this control input matrix yields the heat load increment obtained by each finite element node due to active heating or cooling. Secondly, based on all the outer surface finite element nodes marked as convection and radiation boundary conditions, and combining the convective heat transfer coefficient, surface emissivity, and the boundary area represented by each node, a disturbance input matrix is ​​constructed. Multiplying the environmental disturbance vector with this disturbance input matrix yields the heat load increment obtained by each boundary node due to changes in external thermal conditions.

[0036] The convective heat transfer coefficient and surface emissivity of each node are as follows:

[0037] Convection heat transfer coefficient: By consulting the criteria and correlations for natural and forced convection of air in thermal design manuals or heat transfer textbooks, and combining the characteristic dimensions and surface orientation of the probe chuck, as well as the airflow velocity range in the cleanroom, the theoretical reference range of the convection heat transfer coefficient can be estimated.

[0038] Surface emissivity: Initial values ​​are determined based on the emissivity data tables for corresponding materials and surface conditions in the standard radiative heat transfer data handbook. First, the substrate material is determined according to the material (copper alloy or aluminum nitride, etc.). Then, the corresponding emissivity reference value range is found in the above data table according to the surface treatment state (such as mechanical polishing, natural oxidation, sandblasting, coating with high emissivity coating, etc.). After comprehensive judgment, the initial emissivity values ​​for each outer surface area are given.

[0039] After the above explicit separation is completed, the contributions of the thermal load vector related to the control input vector and the environmental disturbance vector are extracted independently, and the ordinary differential equation system is then transformed into the standard state-space form: the time derivative of the global nodal temperature vector is equal to the thermal system matrix multiplied by the global nodal temperature vector, plus the control input matrix multiplied by the control input vector, plus the disturbance input matrix multiplied by the environmental disturbance vector.

[0040] Considering that the actual controller and control algorithm are both built into the digital processing chip, and their instruction execution and data sampling strictly follow a fixed control cycle, while the above continuous-time state-space equation describes the ideal dynamics within an infinitesimally small time interval, it cannot provide a recursive basis for discrete-time decision-making. Therefore, the above equation is discretized in time with the control cycle as the step size, transforming it into a recursive relationship of the state during adjacent control cycles, resulting in the nominal three-dimensional thermal model;

[0041] The specific control cycle mentioned above needs to be further determined based on the dynamic characteristics and hardware constraints of the probe station temperature control system, taking into account the requirements for temperature signal sampling, actuator refresh frequency, and model computation.

[0042] At the signal sampling level, the fastest mode of temperature change on the chuck surface is determined by the thermal diffusion process near the contact surface of the heating element, and its characteristic time constant is on the order of several seconds. The control cycle should be significantly smaller than this time constant to ensure complete capture of dynamic temperature field change information.

[0043] Actuator refresh: too short a control cycle will shorten the lifespan of power devices and increase electromagnetic interference, while too long a cycle will cause lag in temperature deviation response. A balance must be struck between actuator durability and control response speed.

[0044] The model calculations require sequentially completing the state recursion of the evolutionary three-dimensional thermal model, the prediction and update of the latent parameter set by the ensemble Kalman filter, and the reasoning and verification of the safety control strategy network within each control cycle. Sufficient time margin must be allocated for the above calculation tasks within the control cycle.

[0045] Considering the above constraints, this embodiment sets the control cycle to one hundred milliseconds. This value is much smaller than the thermal time constant of the chuck, matches the sampling refresh cycle of mainstream industrial-grade temperature measurement modules, and the calculation for a single cycle can be completed in real time on the aforementioned embedded industrial control computer platform. If the actual probe station chuck size or material changes significantly, the control cycle can be adaptively adjusted within the range of fifty to two hundred milliseconds according to the proportion of the thermal time constant.

[0046] The time discretization specifically adopts the backward Euler scheme (i.e. implicit Euler) to perform time integration on the ordinary differential equation system: the derivative of the global node temperature vector with respect to time is approximated by the temperature difference between the end time of the current control cycle and the end time of the previous cycle, and the rate of change corresponding to this difference is determined by the temperature distribution at the current time rather than the temperature distribution at the previous time.

[0047] The nominal three-dimensional thermal model is as follows:

[0048] ;

[0049] in:

[0050] : No. The global node temperature vector at the end of each control cycle, i.e., all nodes in the nominal three-dimensional thermal model. The complete spatial temperature distribution formed by the temperature values ​​corresponding to the finite element nodes at the end of the control cycle.

[0051] : No. The global node temperature vector at the end of each control cycle.

[0052] : No. The control input vector applied in each control cycle has a dimension equal to the total number of temperature zones. Each component represents the target heat flow injection amount of the heating and cooling elements in the corresponding temperature zone during the control cycle.

[0053] : No. The environmental disturbance vector for each control cycle includes external variables such as ambient temperature, and is used to introduce thermally influential factors in the external environment that are not controlled by the controller into the nominal three-dimensional thermal model.

[0054] This is a state transition function used to transfer the state of the first state. Global node temperature vector under each control cycle Control input vector Environmental disturbance vector Mapped to the first The global node temperature vector at the end of each control cycle is determined by the finite element discretization and time integration algorithm of the nominal three-dimensional thermal model.

[0055] Process noise is a random vector used to characterize model uncertainty. It is used to quantify the unavoidable deviation between the nominal 3D thermal model and the physical entity, specifically the process noise. The statistical properties of the symmetric matrix are: ... The model follows a Gaussian distribution. Its mean of zero physically means that the model is statistically unbiased; that is, the model's predictions will neither consistently overestimate nor consistently underestimate the actual temperature. Its covariance matrix... The value is obtained offline by comparing the residual statistical characteristics between the model prediction value and the measured value after conducting a step response test on the physical probe card-chuck system, which characterizes the amplitude and spatial distribution of the model uncertainty.

[0056] Regarding the computational scale and real-time performance of the aforementioned nominal 3D thermal model in actual deployment, this embodiment describes the typical meshing parameters of an 8-inch wafer probe station four-temperature zone chuck (approximately 310 mm in diameter and 45 mm in thickness). During finite element mesh generation, the global element size is set to 2.5 mm, with local refinement to 0.8 mm near the probe contact area and the heating element contact surface, ultimately generating approximately 4.2 × 10⁻⁶ tetrahedral elements. 4 There are approximately 8.7 × 10⁻⁶ finite element nodes. 3 One, global thermal conductivity matrix ( The heat conduction pathways corresponding to non-zero, non-diagonal elements in the matrix are approximately 2.1 × 10⁻⁶. 5 strip.

[0057] At the aforementioned grid size, the core computation for a single state recursion is a sparse linear system of equations. Solving for (implicit time integration scheme) =100ms as the control cycle), when using the preconditional conjugate gradient method for solving, the Intel Xeon Gold 6248 processor (single core) takes approximately 1825ms for a single iteration, while OpenMP multi-threaded parallelism (46 cores) can compress it to 6-9ms. In this embodiment, the actual controller uses an industrial-grade embedded industrial PC (Intel Core i7-8700T, 6 cores, 12 threads, 2.4GHz, 16GB DDR4, real-time patched kernel Ubuntu 22.04). The complete control cycle includes step S2. The parallel computation of state recursion and updates for each set member, forward inference of the safety control strategy network in step S3, and verification and correction of the adaptive safety barrier: the total time is approximately 4565ms, less than a 100ms control cycle, with a computational margin of approximately 35%-55%. The state recursion between set members is completely independent. OpenMP is used to evenly distribute the 100 members across 6 physical cores for parallel execution. Combined with the automatic vectorization of sparse matrix operations using the Intel MKL math library, the aforementioned real-time performance is achieved. The mesh size described above is a typical configuration for an 8-inch probe station, balancing spatial resolution of the temperature field (sufficient to resolve heat flow distribution details on the heating element contact surface at approximately 10mm) and real-time computing power. If expanded to a 12-inch probe station (approximately 450mm in diameter), the number of nodes and units is estimated to be approximately 1.8-2.2 times that of the 8-inch platform based on volume ratio. Real-time requirements can be met by increasing the control cycle to 150-200ms or upgrading to a GPU-accelerated platform (such as NVIDIA Jetson AGX Orin, with a CUDA speedup approximately 35 times that of CPU). Those skilled in the art can adjust the global unit size within the range of 2.0-4.0 mm according to the specific probe stage size and temperature control accuracy requirements, so as to achieve an appropriate balance between calculation accuracy and real-time performance.

[0058] Step S2 is specifically used to endow the nominal 3D thermal model with the ability to perceive and follow the implicit degradation of the physical system online, and to identify degradation-sensitive channels in the nominal 3D thermal model that reflect heat transfer links that have a significant impact on the thermal response. Specifically:

[0059] For all heat conduction channels in the global thermal conductivity matrix, the impact of a small change in the thermal conductivity value of each heat conduction channel on the steady-state temperature response of each temperature measurement node on the wafer surface is evaluated one by one; a corresponding sensitivity threshold is preset; the impact degree is compared with the sensitivity threshold, and the heat conduction channel whose impact degree exceeds the sensitivity threshold is defined as a degradation sensitive channel; the sensitivity threshold is as follows: all heat conduction channels are arranged in descending order of their respective impact degree, and the impact degree of each channel is accumulated one by one until the accumulated value reaches 90% of the total impact degree of all channels. The channel that has been accumulated is determined to be a degradation sensitive channel, and the impact degree value at the corresponding cumulative contribution inflection point is taken as the sensitivity threshold.

[0060] The aforementioned level of influence is used to quantify the magnitude of the disturbance to the steady-state temperature distribution of each temperature measurement node on the wafer surface when the thermal conductivity value of any heat conduction channel in the global thermal conductivity matrix changes slightly. Specifically, the nominal three-dimensional thermal model's ordinary differential equations are placed under steady-state conditions, i.e., the time derivative of the global node temperature vector is zero. At this point, the ordinary differential equations degenerate into a system of linear algebraic equations, meaning that under the current thermal load, the heat flow into each node through all heat conduction channels is exactly balanced with the external thermal load received by each node. Solving this system of linear algebraic equations yields the steady-state global node temperature vector, which describes the final temperature distribution that each finite element node in the entire three-dimensional solid space domain eventually tends to, under the condition that the current control input and external disturbance remain unchanged. From this steady-state global node temperature vector, by retaining only the observation matrix corresponding to the temperature measurement nodes on the wafer surface, the steady-state temperature values ​​of each temperature measurement node on the wafer surface are extracted as the benchmark for subsequent sensitivity analysis.

[0061] A non-zero, non-diagonal element connecting node p and node q in the global thermal conductivity matrix is ​​selected; this element represents a specific heat conduction path. When the thermal conductivity of this element changes slightly, the structure of the global thermal conductivity matrix changes accordingly, leading to a change in the steady-state global node temperature vector. The physical meaning of this change is that the altered conductivity of the heat conduction path disrupts the original balance between heat flow and external load at each node, requiring the system to redistribute the temperature of each node to establish a new equilibrium.

[0062] The sensitivity of the steady-state temperature of each temperature sensing node on the wafer surface to the thermal conductivity of the heat conduction channel is determined by two factors. The first factor is the temperature of node q in its current steady state: the higher the temperature of node q, the greater the base amount of heat that can be transferred across the heat conduction channel, and the more significant the impact of its thermal conductivity change on the overall temperature field. The second factor is the system response characteristic described by the p-th column of the inverse of the global thermal conductivity matrix: this column vector reflects the steady-state response distribution of each node's temperature when a unit heat flux is injected at node p; restricting this to the temperature sensing nodes on the wafer surface yields the influence weight of node p on the temperatures of each temperature sensing node on the wafer surface. Finally, multiplying the steady-state temperature of node q by this influence weight vector yields the sensitivity vector of the steady-state temperature of each temperature sensing node on the wafer surface to the thermal conductivity of the heat conduction channel.

[0063] Each component of the sensitivity vector represents the change in steady-state temperature of the corresponding temperature sensing node on the wafer surface as a function of the thermal conductivity of the heat conduction channel. To synthesize the vector-based sensitivity into a comparable scalar index, the Euclidean norm of the sensitivity vector is taken, which is the square root of the sum of the squares of the components. This scalar value reflects the overall magnitude of the steady-state temperature response of all temperature sensing nodes on the wafer surface when the heat conduction channel undergoes a unit change in thermal conductivity; that is, the degree of influence of the heat conduction channel. The greater the degree of influence, the more severe the disturbance to the uniformity of the temperature field on the wafer surface caused by the thermal conduction drift of the heat conduction channel, and the more necessary it is to include it in the online tracking range of latent parameters.

[0064] It should be noted that the above steady-state sensitivity analysis calculates the impact of changes in thermal conductivity of each heat conduction channel on the wafer surface temperature under the condition that the system reaches thermal equilibrium. This is used to spatially locate heat transfer paths that have a global dominant effect on the temperature field distribution. However, the probe card is not always in a steady state during actual testing, and fluctuations in the airflow organization of the cleanroom will continuously change the convective heat transfer conditions on the outer surface of the chuck. This results in the temperature field being in continuous change under the above transient conditions, and the magnitude of the effect of different heat conduction channels in the dynamic response process may differ from the steady-state analysis results. To address this, this embodiment further incorporates a transient activity weighting strategy in the determination of degradation-sensitive channels, ensuring that the screening results simultaneously consider both steady-state dominance and transient participation. During the step response test, specifically in the offline calibration phase of step S2, while applying a step heating power to each temperature zone, the temperature response curves of all thermocouples throughout the transient process are recorded. The root mean square value of the temperature change rate is extracted from each curve as the activity index of that temperature-measuring node during the transient process. The more drastic the temperature change, the more concentrated the heat flow changes along the heat transfer path where the node is located during the dynamic process, and the greater the impact of the accuracy of its thermal conductivity parameters on transient temperature tracking. The transient activity index of each temperature-measuring node is then reversibly distributed to adjacent heat conduction channels according to its spatial position in the finite element mesh, giving channels in areas of drastic temperature change a higher transient weight. Finally, the composite sensitivity of the degradation-sensitive channel is selected, and the influence obtained from the steady-state sensitivity analysis is weighted and synthesized with the aforementioned transient activity weight according to a preset ratio: the steady-state component accounts for 70%, and the transient component accounts for 30%.

[0065] Furthermore, considering that the number of heat conduction channels corresponding to non-zero, off-diagonal elements in the global thermal conductivity matrix is ​​typically tens of thousands to hundreds of thousands, the computational complexity of independently evaluating the impact on each channel would increase linearly with the number of channels, which is unsustainable in practical engineering. This embodiment employs a batch sensitivity analysis method based on matrix perturbation theory: decomposing the impact of minute perturbations in the thermal conductivity of each heat conduction channel into the product of three independent, calculable physical factors:

[0066] One is the temperature difference between the two nodes at the two ends of the channel under the current steady-state temperature distribution, which reflects the base amount of heat flux currently carried by the channel;

[0067] Secondly, the influence weight of the upstream node connected to the channel on the temperature measurement node on the wafer surface reflects the spatial transfer efficiency of heat flow disturbance at the node to the wafer surface.

[0068] Thirdly, there is the direction factor, which reflects the influence of the channel's position in the thermal conductivity matrix on the global temperature field.

[0069] The above three factors can be obtained by solving them all once under the nominal steady-state temperature distribution. Then, the combined calculation of the above factors is performed on all channels. The channels are completely decoupled and the influence of all channels can be evaluated at one time through parallel computing. The total amount of computation is about three to five times that of a single steady-state solution, rather than tens of thousands of times that of solving each channel individually.

[0070] It is worth noting that identifying the degradation-sensitive channel only locates the heat transfer links that significantly affect the thermal response. However, at this point, the thermal conductivity parameters of the degradation-sensitive channel remain fixed nominal values. Therefore, the nominal three-dimensional thermal model does not yet possess the ability to sense and track the degradation state of the physical system. To endow it with the ability to quantitatively track the aforementioned degradation online, ensuring that its thermal response remains synchronized with the physical entity throughout its entire lifecycle, rather than adhering to the nominal thermal characteristics at the time of manufacture, a heater efficiency degradation factor and contact thermal resistance drift coefficient are also embedded as latent parameters in the degradation-sensitive channel. All latent parameters are then stacked in a fixed order to form a latent parameter vector. Specifically:

[0071] The heater efficiency degradation factor is connected in series with the drive terminals of the heating and cooling elements in each temperature zone. Its implementation depends on the scaling input matrix. More details:

[0072] The scaling input matrix is ​​specifically a linear transformation matrix constructed when defining the control input vector to map the control commands of the heating and cooling elements in each temperature zone from the temperature zone control space to the physical space of the finite element nodes. Its construction method is as follows:

[0073] Based on the geometric assembly positions of the heating and cooling elements in each temperature zone within the chuck, the set of finite element nodes covered by the contact surfaces of each heating and cooling element with the chuck body is determined. It is assumed that heat flow is uniformly distributed on the corresponding contact surfaces. The proportion of the contact surface area represented by each node to the total contact surface area is used as the heat flow distribution coefficient. This coefficient is filled into the row of the corresponding node and the column of the corresponding temperature zone in the input matrix. Zero values ​​are filled into the corresponding rows of uncovered nodes. This scaled input matrix is ​​used when the control input vector is multiplied by it. The product result is the incremental vector of heat load obtained by each finite element node due to the active heating or cooling of each temperature zone, thus realizing the spatial mapping from temperature zone control commands to the corresponding finite element node heat load.

[0074] The heater efficiency degradation factor is specifically implemented by embedding a scaling input matrix: When performing the integral calculation of the unit heat load vector, for each unit covered by the temperature zone, the original heat flow injection based on the nominal efficiency is replaced with the effective heat flow injection multiplied by the heater efficiency degradation factor to participate in the integration. After the replacement of all temperature zones is completed and the entire domain is assembled (corresponding to the assembly of the global thermal conductivity matrix and the global thermal capacity matrix), each non-zero element of the column vector corresponding to the temperature zone in the scaling input matrix is ​​scaled proportionally by the heater efficiency degradation factor of that temperature zone. This makes the nodal heat load increment generated by the same control input vector decrease with efficiency decay, thereby truly reflecting the impact of heater efficiency degradation on the temperature field driving capability.

[0075] The contact thermal resistance drift coefficient is distributed at the contact interface between the probe tip and the wafer pad. It is implanted simultaneously in the integral calculation of the unit thermal conductivity matrix: for each unit containing the corresponding contact interface, when calculating its unit thermal conductivity matrix, the interface thermal conductivity corresponding to the reciprocal of the nominal contact thermal resistance is replaced by the actual interface thermal conductivity given by the sum of the nominal contact thermal resistance and the contact thermal resistance drift coefficient. After all contact interfaces are replaced and the entire domain is assembled, the thermal conductivity element values ​​of the nodes on both sides of the corresponding contact interface in the global thermal conductivity matrix are updated. The larger the drift, the smaller the thermal conductivity element value, indicating that the heat conduction channel has a weaker conduction capacity. This causes the corresponding elements in the thermal system matrix obtained by multiplying the global thermal conductivity matrix and the inverse of the global thermal capacity matrix to also change, thus truly reflecting the interface heat transfer degradation caused by probe wear and contact force changes.

[0076] To ensure the integrity of this solution, the specific range and physical constraints of the contact thermal resistance drift coefficient in this embodiment are as follows: Physically, the contact thermal resistance cannot be negative; the sum of the nominal contact thermal resistance and the drift coefficient must always be greater than zero. Based on this physical constraint, the lower bound of the drift coefficient is the negative value of the nominal contact thermal resistance, meaning the drift coefficient must not be less than the opposite of the nominal contact thermal resistance. Otherwise, the actual contact thermal resistance will become negative, and heat flow will be reversed on both sides of the interface, fundamentally violating the second law of thermodynamics. The upper bound of the drift coefficient is theoretically unrestricted; as the probe wears, the contact thermal resistance can gradually increase and approach infinity. However, in actual engineering, when the drift coefficient increases to several times the nominal contact thermal resistance, the contact performance of the probe has severely degraded, and the interface heat transfer capacity has decreased significantly. Continued use will cause the wafer surface temperature field to exceed the uniformity control tolerance. At this point, the probe is no longer worth using and must be replaced within the maintenance cycle.

[0077] The update process of ensemble Kalman filtering adjusts the latent parameter values ​​based on observational information and ensemble statistical covariance, and it does not contain prior encoding of physical constraints such as the non-negativity of contact thermal resistance. Under this unconstrained update mechanism, due to the random fluctuations of observational noise, local model mismatch, or finite sampling errors in ensemble dispersion, the drift coefficients in the ensemble members may be pushed into physically infeasible regions after the update: when the sum of nominal thermal resistance and drift coefficients tends to zero or becomes negative, the integral calculation of the element thermal conduction matrix of the corresponding contact interface in the evolving three-dimensional thermal model will show numerical anomalies with zero or negative denominators. At best, this will cause the heat flow transfer at the contact interface to be erroneously truncated or reversed; at worst, it will cause numerical divergence in the solution process of the entire sparse linear equation system, making the temperature prediction results lose physical meaning and jeopardizing the decision-making of subsequent safety control strategy networks.

[0078] For the reasons mentioned above, this embodiment can add a physical feasibility projection after the analysis step of the ensemble Kalman filter is completed and before the posterior set of latent parameters is used for model prediction in the next cycle. This projection iterates through each member of the posterior set in each control cycle, checking the update results of each contact thermal resistance drift coefficient. If the updated value of a drift coefficient causes the sum of the nominal contact thermal resistance and the drift coefficient to be less than or equal to zero, then the drift coefficient is forcibly projected to the physically permissible lower bound. If a drift coefficient exceeds the engineering upper bound by several times the nominal contact thermal resistance, although it does not violate thermodynamic feasibility, it exceeds the normal degradation range, and the projection operation similarly compresses it to the set engineering upper bound. The projected drift coefficients are placed on the boundary of the feasible region, while the remaining latent parameters within the same set remain unchanged. This maximizes the preservation of effective estimation information obtained from the observation data by the ensemble Kalman filter while maintaining numerical stability and physical interpretability.

[0079] After incorporating the heater efficiency degradation factor and the contact thermal resistance drift coefficient, the nominal three-dimensional thermal model is corrected to a latent parameter dependent form: the heater efficiency degradation factor changes the effective heat flux injection by scaling the element values ​​of the corresponding temperature zone column in the input matrix; the contact thermal resistance drift coefficient affects the local heat transfer capacity by modifying the thermal conductivity element values ​​between the finite element nodes on both sides of the contact interface in the thermal system matrix. The combined effect of both is ultimately reflected in the changes in the matrix element values ​​of the input matrix and the thermal system matrix as the latent parameter values ​​change. Since the actual controller and control algorithm operate in the digital system with a fixed control cycle, the continuous-time state-space equation cannot be directly used for online recursive calculation. Therefore, time discretization is performed with the control cycle to obtain an evolutionary three-dimensional thermal model containing latent parameters, as follows:

[0080] ;

[0081] in:

[0082] For the first The global node temperature vector at the end of each control cycle; For the first The global node temperature vector at the end of each control cycle; For the first The control input vector applied in each control cycle has a dimension equal to the total number of temperature zones. ; For the first The environmental disturbance vector for each control cycle includes external variables such as ambient temperature; : State transition function; This refers to process noise; the specific definitions of the above parameters are the same as those of the parameters corresponding to the nominal three-dimensional thermal model in step S1.

[0083] : No. The latent parameter vector for each control cycle contains the heater efficiency degradation factor for all temperature zones. Contact thermal resistance drift coefficient of all probe contact areas ;

[0084] The latent parameter vector in step S2 is as follows:

[0085] First, initial latent parameters are obtained based on the manufacturer's specifications or material handbooks for each component. The initial value of the heater efficiency degradation factor is set to 1, corresponding to the state where the heater has not experienced any efficiency degradation at the time of manufacture; the initial value of the contact thermal resistance drift coefficient is set to 0, corresponding to the state where the contact thermal resistance has not experienced any drift at the time of manufacture. These initial latent parameters are only theoretical nominal values. Due to manufacturing tolerances and assembly differences, the actual initial state of each physical probe card-chuck system will deviate from the above theoretical nominal values. If these are directly used as the starting point for online recursive estimation, the latent parameters will converge slowly in the initial stage or even deviate from the actual degradation trajectory. Therefore, it is necessary to perform offline calibration of the initial latent parameters through step response testing to obtain a statistical distribution that more closely approximates the actual initial state of the physical system. The specific operation of the step response test is as follows: A heating power step of known amplitude is applied to each temperature zone sequentially, and the transient temperature response curves corresponding to each temperature zone are recorded simultaneously. The measured curves are then fitted to the predicted response of the nominal three-dimensional thermal model under the same step excitation using a system identification method. An optimization algorithm is used to adjust the values ​​of the latent parameters to minimize the deviation between the predicted and measured values. After calibration, the initial prior mean and initial prior covariance of the latent parameter vector are obtained. The specific working principle of the system identification method is as follows:

[0086] Record No. When a step excitation is applied to each temperature zone individually, at the sampling time... ( , (Total number of sampling points) The measured temperature value collected by the embedded thermocouple is... The temperature prediction output from the simulation of the nominal three-dimensional thermal model containing latent parameters at the same step input and the same node location is... ,in Let `t` be the latent parameter vector for the current iteration step. Define a scalar objective function. Sum of squared prediction biases for all temperature zones at all sampling times ,in This represents the total number of temperature zones. The goal of system identification is to solve for... Minimize the optimal latent parameter vector ;

[0087] Solving the above optimization problem can be done using nonlinear least squares algorithms such as the Gauss-Newton method or the Levenberg-Marquardt method. Taking the Gauss-Newton method as an example, its iterative update formula is: ,in For the first The latent parameter vector of the next iteration The prediction bias for all temperature zones at all sampling times The residual vector formed by stacking in a fixed order. Let be the Jacobian matrix of the residual vector with respect to the latent parameter vector, where each row corresponds to a residual component and each column corresponds to a latent parameter component. The matrix elements are... That is, the first The residual component for the first... The partial derivatives of each latent parameter are used to characterize the sensitivity of the model's predicted temperature to changes in each latent parameter.

[0088] After iterative convergence, the initial prior mean of the latent parameter vector is taken as... Initial prior covariance matrix Determined by the approximate Hessian matrix of the objective function at the convergence point: ,in This is an unbiased estimate of the residual variance. Let be the dimension of the latent parameter vector. The resulting initial prior mean is... and initial prior covariance ;

[0089] Then from the initial prior mean The mean and initial prior covariance are given. Independent random sampling in a Gaussian distribution with covariance Each sampling yields a result from... The construction of the latent parameter value vector includes an initial set of latent parameters with multiple parameters to be selected. Each member in the initial set of latent parameters is assigned a corresponding initial value of the global node temperature vector. Specifically, the initial value is taken as the steady-state temperature field at the end of the step response test, that is, a column vector formed by arranging the temperature values ​​of all finite element nodes at this moment in a fixed numbered order.

[0090] After running online, in the first... For each control cycle, obtain the corresponding multi-temperature zone temperature measurement value vector, and use it as the initial set of latent parameters. Each parameter to be selected in The corresponding multi-temperature zone temperature measurement vector These parameters are input into the evolutionary three-dimensional thermal model in step S2, and the model is calculated forward for one control cycle to obtain the predicted temperature values ​​corresponding to each parameter to be selected. A temperature prediction value; ultimately determined by... Each temperature prediction value constitutes a set of prior prediction states for the temperature field.

[0091] Because the latent parameter values ​​carried by each parameter to be selected in the initial set of latent parameters are different, the above The predicted temperature values ​​are dispersed in the state space. The degree of dispersion reflects the uncertainty of the three-dimensional thermal model's estimation of the temperature state before the measured values ​​arrive. The statistical structure of this uncertainty is the basis for calculating the ensemble Kalman gain. Therefore, the average predicted temperature is calculated based on all predicted temperature values ​​for all selected parameters.

[0092] An observation matrix is ​​pre-constructed based on the installation locations of physical sensors within the finite element nodes. Each row of this matrix corresponds to a sensor installation node, with the element in the column corresponding to the node's global number set to 1 and the elements in the remaining columns set to 0. The observation matrix maps the predicted temperature values ​​of each parameter to the sensor locations, yielding the observed predicted values ​​and their ensemble average. Furthermore, ensemble statistics are used to obtain the state-observation cross-covariance matrix and the observation autocovariance matrix. The result is then calculated from these two matrices to obtain the... The Kalman gain over one control cycle ,in The state-observation cross-covariance matrix, To observe the autocovariance matrix, For the ensemble average predicted temperature, To measure noise The covariance matrix;

[0093] Call out the first Vector of measured temperature values ​​in multiple temperature zones under each control cycle Calculate the observed predicted value for each parameter to be selected. Vector of measured temperature values ​​in multiple temperature zones The deviation between these parameters is taken as the innovation. Since the innovation reflects the degree of deviation between the prediction of the current temperature state by each set member and the measured value of the physical entity, but the innovation itself only exists in the sensor observation space and cannot be directly used to correct the high-dimensional global node temperature vector and latent parameter vector, the ensemble Kalman gain is used to back-map the innovation into the state correction amount of each parameter to be selected, thus obtaining the updated global node temperature vector. ,in For a mean of zero and a covariance of The perturbation of random sampling in a Gaussian distribution is used to maintain the statistical dispersion of the set.

[0094] Meanwhile, the latent parameter vector As part of the augmented state, it participates in the above update together with the global node temperature vector, in which the latent parameter value of each parameter to be selected is taken. Naturally, a posteriori correction is obtained, thus yielding the first... The updated posterior set of latent parameters under each control cycle That is, the set of latent parameters corresponding to step S2;

[0095] To ensure the completeness of this solution, the following explanation is provided regarding the selection and engineering processing of several key parameters for ensemble Kalman filtering in practical deployment:

[0096] Number of members in the set The determination of the covariance matrix requires a balance between estimation accuracy and computational burden: if the number of members is too small, the covariance matrix calculated from a finite sample will have large sampling noise, and the matrix condition number will deteriorate or even approach singularity; if the number of members is too large, the computational load will increase linearly, which may exceed the constraints of the real-time control cycle. In this embodiment, the number of set members is set to one hundred, which is about two to five times the dimension of the latent parameter vector (twenty to fifty), which is sufficient to ensure reliable estimation of the covariance matrix. At the same time, on the aforementioned embedded industrial control computer platform, the forward recursion and update calculation of all members can be completed within a hundred millisecond control cycle.

[0097] The construction of the initial prior covariance of the latent parameters adopts a conservative amplification strategy: the dynamic range excited by the step response test is mainly based on the dominant thermal mode of the system, and its frequency band is mainly concentrated in the low-frequency range. The observability of the high-frequency dynamic related part of the latent parameters is limited, and the initial covariance matrix obtained by a single step test may be too small. In this embodiment, based on the covariance estimate obtained by approximating the Hessian matrix, different amplification magnitudes are set for the latent parameters with different physical properties: the initial standard deviation of the heater efficiency degradation factor is amplified by 150% to 200% of the estimated value, and the initial standard deviation of the contact thermal resistance drift coefficient is amplified by 200% to 300% of the estimated value, so that the prior distribution has a wider confidence interval.

[0098] Process noise covariance matrix and measurement noise covariance matrix The results were determined through offline estimation and static calibration, respectively. The residual sequence between the measured temperature values ​​and the model prediction values ​​during the step response test is used for estimation: During the step response test, the measured temperatures of all thermocouples and the predicted temperatures of the nominal three-dimensional thermal model under the same excitation are recorded simultaneously. The difference between the two is used to obtain the residual vector at each sampling time. The residual sequence is divided into several segments according to time intervals, and the statistical average of the covariance of each segment is taken as the residual vector. The initial estimate. While the probe station is in a constant-temperature stable state, the output readings of each thermocouple channel are continuously acquired. The sampling lasts from tens of seconds to hundreds of control cycles. The standard deviation of each channel reading in the time dimension is calculated, and its square is taken as the measurement noise variance. The variances of each channel are arranged along the diagonal to form the standard deviation of the measurement noise variance. matrix.

[0099] The random perturbation injection method in the ensemble Kalman filter update step is to add random perturbations with the same distribution as the measurement noise to the observations. The disturbance has a mean of zero and a covariance of Independent sampling is performed within a Gaussian distribution, affecting the calculation of information for each member. Its function is to maintain the statistical dispersion of the set members after updates; without perturbation injection, each member converges to the measured value after each update, the set dispersion gradually shrinks, the covariance matrix tends to degenerate, the filter gain decreases accordingly, and the corrective effect of subsequent measured values ​​on the state is gradually lost. The physical meaning of observation perturbation is that the simulated actual measurement value itself contains random errors; each set member faces a noisy real measurement value, thus maintaining the set covariance at a reasonable level.

[0100] It is worth further explaining that the posterior set of latent parameters obtained by ensemble Kalman filtering update... It is not merely a mathematical moderating variable within the model, but rather has a verifiable and traceable physical correspondence with the actual degradation state of the physical hardware. Specifically:

[0101] Heater efficiency degradation factor Between 0 and 1 Corresponding to the The heating element in the temperature zone is operating at its factory-rated efficiency. Each decrease of 0.01 indicates a 1% drop in the electrothermal conversion efficiency of the heating element in that temperature range relative to its nominal value. In practical engineering, the heater efficiency degradation factor... A quantitative correlation can be directly established with the measured DC resistance value of the heating element. Specifically, for thin-film heating elements, their resistance value... The output power gradually increases with the increase of cumulative energizing time and thermal cycling count. Under constant voltage drive, the output power... It is inversely proportional to the resistance value, therefore the efficiency degradation factor Specifically During regular maintenance, the actual resistance values ​​of the heating elements in each temperature zone can be obtained through four-wire resistance measurement and compared with the efficiency degradation factor estimate output by the ensemble Kalman filter. If the deviation between the two exceeds a preset threshold, it indicates that the model estimation may be affected by unmodeled factors, requiring model verification or recalibration.

[0102] Contact thermal resistance drift coefficient Indicates the first The increase in actual contact thermal resistance relative to the nominal value at the probe contact interface. In actual probe card maintenance, the contact thermal resistance drift coefficient... The wear and indentation depth of the probe tip are directly related to the cumulative number of probe drops, the amount of tip wear, and the cleanliness of the probe. Regularly measuring the wear area and indentation depth of the probe tip using an optical microscope or scanning electron microscope allows us to apply the inverse relationship between contact thermal resistance and contact area. The change in wear area is converted into an increase in contact thermal resistance, which is then compared with the ensemble Kalman filter estimate. Perform cross-validation.

[0103] Mean of the posterior set of latent parameters The optimal estimate of the current degradation state is directly used in step S3 to update the system matrix and input matrix of the evolving 3D thermal model, ensuring that the model's thermal response characteristics remain consistent with the entity system. The covariance matrix of the posterior set... The uncertainty of the current estimate is quantified by the square root of the estimated variance of each latent parameter corresponding to the diagonal elements of the covariance matrix. In engineering applications, when the estimated standard deviation of a degradation factor exceeds 20% of its current mean, it indicates that the measured temperature value is insufficient to identify the degradation factor. This may be due to strong coupling between the degradation factor and other parameters or insufficient sensor placement to provide sufficient observational information. In this case, an early warning of insufficient parameter identification can be issued to the operator, suggesting increasing the density of temperature sensors in the area or performing a specialized step excitation test to improve identification accuracy.

[0104] Simultaneously, the estimated values ​​of each degradation factor output from the Kalman filter can be directly mapped to specific hardware maintenance actions. The specific mapping rule is: when the heater efficiency degradation factor for any temperature zone... When the contact thermal resistance coefficient remains below 0.85 and shows a monotonically decreasing trend, it is determined that the heating element in that temperature zone has reached the end of its lifespan, and the system outputs a maintenance recommendation to replace the heating element in that temperature zone; when the contact thermal resistance coefficient of any contact area drifts... When the increase exceeds 50% of the nominal contact thermal resistance, it is determined that the contact performance of the probe has been severely degraded, and maintenance recommendations for cleaning or replacing the probe are output. When the degradation factors of multiple temperature zones or contact areas show abnormal changes at the same time, it is further determined that there may be a common failure mode (such as insufficient coolant flow or abnormal ambient temperature), and corresponding systemic troubleshooting recommendations are output.

[0105] It should be noted that the location of degradation-sensitive channels and the inversion of latent parameters follow a causal chain of measurable physical effects, spatial location of channels, and numerical identification of parameters. To facilitate understanding of the complete mapping relationship between finite element nodal temperature measurements and numerical estimation of heater efficiency degradation factors, the following explanation is presented from both the forward and inverse problem perspectives.

[0106] At the core of the problem: how does the degradation effect manifest in measurable temperature? The heater efficiency degradation factor is specifically the ratio of the actual output heat power of the heating element to the nominal drive power. When the heater resistance wire undergoes microstructural changes due to high-temperature oxidation and thermal cycling fatigue, its resistivity increases and its heating efficiency decreases. At this time, even if the PWM duty cycle given by the control command remains unchanged, the actual heat flux density injected into the chuck body at the contact surface also decreases. The decrease in heat flux density, as a change in boundary conditions, propagates throughout the three-dimensional solid space domain via the finite element heat conduction equation (the discrete form of the partial differential equation described by the global thermal conductivity matrix and the global thermal capacity matrix) established in step S1, ultimately manifesting as a systematic reduction in the temperature response amplitude of each temperature measurement node. Similarly, the physical meaning of the contact thermal resistance drift coefficient is the increment of the actual contact thermal resistance at the interface between the probe tip and the wafer pad relative to the nominal value. As the number of probe drops increases, wear, oxide film thickening, and particulate matter adhesion occur on the probe tip surface, reducing the actual contact area and increasing the resistance to heat flow transfer at the interface. The change in interfacial thermal resistance directly alters the thermal conductivity element values ​​of the corresponding node pairs on both sides of the contact interface in the element's thermal conductivity matrix, thereby changing the steady-state distribution and dynamic response characteristics of the entire temperature field. Both of these degradation effects, described by the finite element model as heat conduction physical processes, ultimately leave measurable response characteristics at the observable finite element node temperatures.

[0107] Inverse problem level: The degradation factor is deduced from the temperature measurement. The ensemble Kalman filtering used in step S2 is specifically a sequential data assimilation method, which is used to transform the latent parameter vector... (Including all heater efficiency degradation factors and contact thermal resistance drift coefficients to be estimated) and global nodal temperature vector The augmented state vector is constructed by combining the measured temperature values. The augmented state is sequentially updated. Specifically, in each control cycle, the ensemble Kalman filter performs two steps: prediction and analysis. Prediction uses the evolutionary three-dimensional thermal model under the current latent parameter values ​​to advance the augmented state of the previous cycle to the current time, obtaining the ensemble distribution of temperature prediction values. Analysis, based on the deviation (news) between the measured and predicted temperature values, uses the Kalman gain calculated by ensemble statistics to back-project the news to the high-dimensional augmented state space, while simultaneously correcting the global node temperature vector and latent parameter vector.

[0108] The aforementioned back projection specifically involves using the covariance information implicit in the distribution of ensemble members in the state space to decompose the one-dimensional bias signal in the observation space into corrections for each latent parameter component. Taking the heater efficiency degradation factor as an example, when the measured temperature in a certain temperature zone is systematically lower than the model prediction, the component in the ensemble Kalman gain corresponding to that temperature zone will convert the negative bias in the innovation into a negative correction of the heater efficiency degradation factor for that temperature zone (i.e., a decrease in the efficiency factor value). This reduces the effective heat flux injection in that temperature zone in the model prediction of the next cycle, causing the model temperature response to approach the measured value. Similarly, when the measured temperature near a probe contact area shows a trend that is inconsistent with adjacent areas (higher or lower), the gain matrix will map this local temperature anomaly into a corresponding correction of the contact thermal resistance drift coefficient of that contact area. Since there is a definite correspondence between the implantation location of the latent parameters (the heater efficiency degradation factor acts on the corresponding temperature zone column of the input matrix, and the contact thermal resistance drift coefficient acts on the corresponding contact interface node pair of the thermal conductivity matrix) and the spatial location of the measured temperature value, each set of latent parameter values ​​estimated by the ensemble Kalman filter has a clear spatial positioning meaning: the efficiency degradation factor belongs to the heating element of a specific temperature zone, and the contact thermal resistance drift coefficient belongs to the contact interface of a specific probe.

[0109] Step S3, based on step S2, the first Based on the evolutionary three-dimensional thermal model and latent parameter set under each control cycle, a corresponding safety control strategy network with embedded adaptive safety barriers is established.

[0110] The safety control strategy network is specifically an actuator network composed of deep neural networks, wherein:

[0111] The actuator network receives observable temperature states as input, passes through a multi-layer fully connected neural network for forward propagation, and outputs a deterministic candidate action vector. The dimension of this vector equals the total number of temperature zones, and each component represents the target power command for the heating and cooling elements in the corresponding temperature zone. Nonlinear activation functions are used between the hidden layers, enabling the actuator network to fit the complex mapping relationship between high-dimensional temperature states and optimal power allocation. Finally, the actuator network outputs the candidate action vector.

[0112] Because the direct output of the evolutionary 3D thermal model is a high-dimensional global node temperature vector, which contains the temperature values ​​of all finite element nodes in the entire 3D solid space domain, but in a real probe station controller, temperature feedback at a specific location can only be obtained through a limited number of embedded thermocouples. To ensure that the temperature information conditions relied upon by the safety control strategy network are completely consistent during the virtual training and physical deployment phases, the observation matrix in step S2 is used again to map the global node temperature vector output by the evolutionary 3D thermal model to a low-dimensional observable temperature state corresponding one-to-one with the physical sensors;

[0113] The safety control strategy network maps observable temperature states to specific power commands for heating and cooling elements in each temperature zone. Specifically, this network is an actuator network composed of a deep neural network, consisting of an input layer, multiple hidden layers, and an output layer. Information is transmitted between these layers via full connections. Its function is to map the observable temperature states to specific power commands for heating and cooling elements in each temperature zone, thereby completing end-to-end decision-making from limited sensor temperature information to multi-temperature zone control actions.

[0114] The input layer receives the first Observable temperature state for each control cycle As input. The input layer contains There are 10 neurons, each neuron receiving a corresponding... One component is responsible only for receiving data and does not perform any mathematical transformations.

[0115] Hidden layer by The first fully connected layer is connected in series. The number of neurons in each hidden layer is denoted as . The output vector of each hidden layer The output vector from the previous layer is obtained through a linear transformation and then through a nonlinear activation function. The first hidden layer receives the observable temperature state transmitted from the input layer. As input, its output is ,in This is the weight matrix of the first hidden layer, with dimension 1. Each row corresponds to the connection weights between a neuron in the first hidden layer and each neuron in the input layer. Let be the bias vector of the first hidden layer, with dimension . , where is the offset of each neuron in the first hidden layer; It is a nonlinear activation function, which introduces nonlinear transformation capability into the network, enabling the network to fit the complex mapping relationship between temperature state and power command.

[0116] For the Hidden layers ( ), its input is the first Output of hidden layer The output is ,in For the first The weight matrix of the hidden layer has a dimension of . Each row corresponds to the first A neuron in the hidden layer and the first Connection weights of neurons in each layer; For the first The bias vector of the hidden layer has a dimension of . .

[0117] The output layer is a fully connected layer with a certain number of neurons. This equals the total number of temperature zones. The output layer receives the output of the last hidden layer. As input, it undergoes a linear transformation to produce candidate action vectors. ,in Here is the weight matrix of the output layer, with dimension 1. Each row corresponds to a temperature zone and the connection weights of neurons in the last hidden layer. Let be the bias vector of the output layer, with dimension . Each component corresponds to an output offset for a specific temperature range. Candidate action vectors. Each component ( ) represents the first The target power command for each temperature zone heating and cooling element in the current control cycle. To ensure that its numerical range conforms to the upper and lower limits of the physical power of each temperature zone heating and cooling element, a range clipping operation can be added after the output layer to limit each component to within... Within the range.

[0118] In summary, the security control strategy network starts from the observable temperature state. to candidate action vectors The complete forward propagation process can be represented as ,in This is the set of all parameters to be optimized for the security control policy network. Indicated by network parameters The defined mapping function, through layer-by-layer linear transformation and nonlinear activation, gradually transforms the high-dimensional observable temperature state into the target power command for each temperature zone.

[0119] The training process for the aforementioned actuator network is as follows:

[0120] The training of the safety control policy network is an offline reinforcement learning process performed within the evolutionary 3D thermal model. The goal of training is to iteratively update the network parameters. This ensures that the candidate action vectors output by the network maximize the expected cumulative discounted reward starting from the current state, while always satisfying the constraints of the adaptive safety barrier. Training employs a deterministic policy gradient-based approach, executing the following steps in each virtual control cycle.

[0121] Action generation and security verification:

[0122] In the current virtual control cycle, the evolving three-dimensional thermal model will use the global node temperature vector. Observation matrix Mapped to observable temperature state This serves as the input to the security control policy network. The security control policy network outputs candidate action vectors through forward propagation. The candidate action vector must first be submitted to the adaptive safety barrier for verification. It will only be released for execution if it meets the current adaptive stability criterion; otherwise, it will be discarded, and the safety control policy network will need to correct it along the negative direction of the gradient of the control Lyapunov function and re-output it until it passes the verification.

[0123] Action execution and reward acquisition:

[0124] The validated candidate action vectors are used as the control input for the current cycle. , and the mean of the current posterior set of latent parameters Environmental disturbance vector By inputting the state transition function of the three-dimensional thermal evolution model, the global nodal temperature vector at the next time step can be obtained. The environment from The maximum temperature difference on the wafer surface, the total energy consumption across multiple temperature zones, and the temperature recovery and stabilization time are extracted and substituted into the composite reward function to calculate the scalar reward signal. .at the same time, The observable temperature state at the next moment is obtained after mapping the observation matrix. The sample of this interaction Store in the experience replay buffer.

[0125] Policy network updates:

[0126] The security control policy network employs a deterministic policy gradient algorithm for parameter updates. The value estimate required for the update is not provided by an independent judge network, but rather directly estimated from the actual reward sequence of samples in the experience replay buffer. For a batch of samples randomly sampled from the buffer, the value of each action is... The value is calculated by accumulating the immediate reward of the sample with the discounted reward of subsequent states. The optimization objective of the security control policy network is to maximize the mean value of all actions within the sampling batch. .in This is the batch size, i.e., the number of samples used in a single parameter update; For the first The action value corresponding to each sample represents the value in the state. Next action The expected cumulative discount reward obtained afterward.

[0127] Network parameters Update the parameters along the gradient direction of the objective function, specifically as follows: ,in The learning rate controls the step size of each parameter update. For the objective function For network parameters The gradient indicates the direction in which the parameters should be adjusted to increase the objective function value. The gradient of this policy is calculated as follows: for each state in the batch sampling... Calculate the gradient of the action value with respect to the network's output action. This gradient indicates in which direction the action should be fine-tuned to increase value; then, using the chain rule, this gradient is backpropagated to the parameters of each layer of the network to obtain the gradient of the objective function with respect to each parameter. The mean of the gradients of each sample is the policy gradient.

[0128] ;

[0129] in Indicates the value of an action Action The gradient in The value at that point reflects how small changes in the action will affect its value assessment in the vicinity of the current policy output. Output network parameters for the policy network The gradient of the objective function reflects how changes in network parameters will affect the output action; multiplying the two and averaging them gives the gradient estimate of the objective function with respect to the network parameters.

[0130] Policy iteration convergence:

[0131] The three steps described above (action generation and security verification, action execution and reward acquisition, and policy network update) are executed cyclically within the evolutionary 3D thermal model, with network parameters updated after each virtual control cycle. As training progresses, the candidate action vectors output by the security control policy network gradually converge from random exploration to the optimal policy that maximizes cumulative reward under adaptive security barrier constraints. The parameters of the security control policy network after training are shown below. It is fixed, compressed by distillation, and then deployed to the actual controller for online temperature control decisions;

[0132] Obtain the global node temperature vector output by the evolutionary three-dimensional thermal model, and map the global node temperature vector to the observable temperature state through the observation matrix in step S2, which serves as the control input vector of the safety control strategy network; the safety control strategy network outputs the control input vector for the current period based on the observable temperature state.

[0133] When the safety control strategy network outputs the control input vector for the current period, it first outputs the candidate action vector and then uses an adaptive safety barrier to test the candidate action vector for an adaptive stability criterion. If the test is passed, the candidate action vector is used as the control input vector for the current period.

[0134] If the test fails, the candidate action vector is discarded, and the security control policy network generates a new candidate action vector and submits it for test again until it passes the test of the adaptive security barrier.

[0135] The security control policy network with embedded adaptive security barriers in step S3 is as follows:

[0136] Set a sliding window of fixed length, and perform two statistical operations on the prediction residual vectors of multiple consecutive control cycles within the sliding window: the first is the root mean square value of all residual vectors within the sliding window, and the second is the root mean square value of the difference between the residuals of adjacent cycles within the sliding window; the two statistical operations are weighted and synthesized according to preset weights to obtain the model prediction error sequence for the current cycle.

[0137] Regarding the prediction residual vector, this embodiment uses the vector formed by the difference between the predicted and measured temperatures at the corresponding positions of all embedded thermocouples within each control cycle as the prediction residual vector. The dimension of this vector is equal to the total number of physical thermocouples, rather than the dimension of the global node temperature vector. This is because the actual controller can only obtain temperature feedback information from the sensor deployment locations, and the number of sensors deployed is far less than the total number of finite element nodes (sixteen thermocouples are deployed in this embodiment). Therefore, the prediction residual vector is a low-dimensional vector, and its computational and storage requirements are far less than those of the high-dimensional node temperature residual scheme. At the same time, it is consistent with the observation space on which the observation update of the ensemble Kalman filter in step S2 is based.

[0138] The low error interval threshold is determined by the residual statistical characteristics of the step response test. After completing the step response test in step S2, the temperature prediction value of the nominal three-dimensional thermal model under the same step excitation is subtracted from the measured value at the corresponding time, resulting in a residual vector sequence for all sampling times. The absolute value of the residual for each channel at each time is taken, and its mean and standard deviation are calculated in the time direction. The sum of the mean and twice the standard deviation is used as the upper limit threshold of the low error interval. The principle for this value is: under normal model prediction accuracy, approximately 95% of the prediction residuals should fall within this threshold range. Prediction residuals exceeding this range indicate a significant deviation between the model prediction and the actual measurement, requiring the triggering of the safety barrier tightening criterion. This threshold can be further fine-tuned in the early stages of online operation based on the actual statistical characteristics of the measured residuals.

[0139] The maximum prediction deviation is obtained by subtracting the low-error interval threshold from the current model prediction error sequence. Physically, it represents the magnitude by which the current model prediction error exceeds the normal confidence interval. This magnitude characterizes the degree to which the model's reliability in predicting temperature decreases at the current moment—the greater the error exceeds the threshold, the less reliable the model prediction, and the greater the potential risk of the control action assessed based on this prediction. The theoretical basis for using the excess error magnitude as the maximum prediction deviation is as follows: when the deviation between the model prediction and the measured value does not exceed the normal confidence interval, the prediction reliability is sufficiently high, and it is only necessary to ensure that the control Lyapunov function value decreases; when the deviation exceeds the confidence interval, the unreliable component in the prediction is positively correlated with the excess magnitude. The safety criterion requires that the energy decay generated by the control action at least covers the prediction error that may be caused by this unreliable part, i.e., the excess magnitude is used as the minimum requirement for safety margin.

[0140] The candidate action vector output by the actuator network is used as the control input. Together with the global node temperature vector output by the nominal three-dimensional thermal model in the current control cycle of step S2, the mean of the latent parameter set in step S2, and environmental disturbances, these are input into the evolving three-dimensional thermal model in step S2. The model is then calculated one control cycle backward to obtain the predicted value of the global node temperature vector in the next control cycle after executing the corresponding candidate action vector. The control Lyapunov function values ​​at the current and next time moments are calculated, and the difference between the two is obtained to obtain the change in the generalized energy of the system before and after executing the candidate action vector, which is the energy change.

[0141] ;

[0142] in:

[0143] To control the value of the Lyapunov function, The target temperature vector represents the desired temperature value at the corresponding finite element node. Specifically, in a uniform temperature field, all nodes take the same set value; in a non-uniform temperature gradient distribution, each temperature zone has an independent target, the set value of the temperature zone is evenly distributed among the nodes within the temperature zone, and the nodes at the boundary are weighted by inverse distance interpolation based on the distance from the contact surface of the heating element in the adjacent temperature zone. The set temperature difference between adjacent temperature zones needs to be verified based on the maximum allowable thermal stress of the chuck material. If it exceeds the thermal shock tolerance limit, the set temperature is iteratively and smoothly adjusted until the thermal stress constraint is met. Given a symmetric positive definite weight matrix, the solution is obtained by finding the matrix of the nominal thermal system. Related Lyapunov algebraic equations get, This is a preset positive definite matrix used to determine... The specific value;

[0144] Regarding nominal values ​​under degradation conditions Solving Whether stability can still be guaranteed is illustrated in the following example: when the system matrix is ​​perturbed. At that time, as long as The norm is less than that of and Determined stability margin, matrix inequalities That is, it is established. This embodiment selects... For a positive definite matrix with large diagonal elements, such that It has a large gain to reserve a stability margin that can cover system matrix perturbations where heater efficiency degrades to 70% of the nominal value and contact thermal resistance drifts to twice the nominal value. If the degradation exceeds this range, the current posterior mean is used. Refactoring And resolve the Lyapunov equations online to update .

[0145] To generate the predicted value of the global node temperature vector at the next time step after executing the candidate action vector in the evolutionary three-dimensional thermal model. After substituting the control Lyapunov function, the predicted generalized energy of the system at the next moment after executing the candidate action vector is calculated. ;

[0146] Obtain the low-error interval; compare the model's predicted error sequence with the low-error interval, and analyze whether the candidate action vector passes the verification based on the comparison results:

[0147] The validation passes when the model's prediction error sequence is in the low error range and the energy change is negative.

[0148] When the model prediction error sequence exceeds the low error range, if only the energy change is required to be negative, the dangerous action may be misjudged as safe due to the inaccuracy of the evolutionary three-dimensional thermal model itself. Therefore, the verification standard is automatically tightened to require that the magnitude of the energy change must exceed the maximum prediction deviation mapped by the current model prediction error sequence.

[0149] When a candidate action vector is determined by the adaptive safety barrier to not meet the adaptive stability criterion, the adaptive safety barrier simultaneously calculates the gradient of the current control Lyapunov function value with respect to each component of the candidate action vector. The control Lyapunov function is a function of the predicted global node temperature vector value at the next time step, and the predicted global node temperature vector value at the next time step is a function of the candidate action vector. Therefore, by using the chain rule, the gradient of the control Lyapunov function with respect to the global node temperature vector is multiplied by the Jacobian matrix of the predicted global node temperature vector value with respect to the candidate action vector to obtain the gradient of the control Lyapunov function with respect to the candidate action vector.

[0150] Furthermore, the gradient of the control Lyapunov function with respect to the global nodal temperature vector is directly and analytically given by the definition of the control Lyapunov function, specifically as follows: , To control the gradient of the Lyapunov function with respect to the global node temperature vector, according to the definition of controlling the Lyapunov function, the gradient of the global node temperature vector can be analytically expressed as: ; Let be the Jacobian matrix of the predicted global node temperature vectors for the next time step against the candidate action vectors, and its elements are... Indicates the first The predicted temperature of the finite element node affects the ... The sensitivity of power commands in each temperature zone is determined by the state transition function of the evolving three-dimensional thermal model. The local linearization at a given point can be obtained by calculating using the finite difference perturbation method or automatic differentiation techniques.

[0151] In more detail, this embodiment uses the finite difference perturbation method to calculate the Jacobian matrix. : Using the current candidate action vector Based on this, a small perturbation is applied to each component (corresponding to the power command for each temperature zone). The perturbed action vector is then sequentially input into the evolutionary three-dimensional thermal model to perform state recursion, yielding the corresponding global node temperature prediction vectors. The difference between the perturbed temperature prediction vector and the baseline temperature prediction vector is divided by the perturbation amount to obtain the corresponding column of the Jacobian matrix. The perturbation amount is set as small as possible while ensuring that the temperature response change is higher than the numerical rounding error; typically, it is taken as one-thousandth to one-hundredth of the power range for each temperature zone. This calculation is performed... Sub-state recursion ( This represents the total number of temperature zones in this embodiment. Each iteration solves an independent system of sparse linear equations, and the iterations are completely decoupled from each other. This can be executed in parallel on the aforementioned embedded industrial control computer using OpenMP.

[0152] Specifically, the gradient calculated above... For one dimensional vector, This equals the total number of temperature zones. Its first... The component represents the first... The rate of change of the system's generalized energy prediction when the power command changes slightly in each temperature zone. The positive direction of the gradient indicates the direction in which the system's generalized energy increases the fastest, while the negative direction indicates the direction in which the system's generalized energy decays the fastest. Therefore, the safety control policy network, starting with the rejected candidate action vectors, corrects along the negative direction of the gradient to make the system's generalized energy decay most rapidly.

[0153] The safety control policy network corrects the candidate action vectors along the negative direction of the gradient, generating new candidate action vectors. ,in Let be the Euclidean norm of the gradient; To correct the stride, it corresponds to the lower limit of energy decay required by the stringent criterion in the adaptive stability criterion, i.e., the current model prediction error sequence value. Exceeding the low error range threshold The maximum prediction bias mapped by the part (when (when), and the above model predicts the error sequence value The low error range threshold is calculated and output in real time for each control cycle. The adaptive safety barrier is determined based on the residual statistical characteristics of the step response test. When the value is large, the maximum prediction deviation is correspondingly large, and the correction step size increases accordingly, allowing for a faster traversal of the ambiguity zone with high model uncertainty; when near When the maximum prediction deviation is smaller, the correction step size is reduced accordingly, allowing for a more precise approach to the safety boundary.

[0154] The above-mentioned correction step The maximum deviation in the prediction of temperature dimension needs to be calculated. The correction factor for converting to power dimensions is introduced in this embodiment by introducing a conversion coefficient. To achieve this dimensional transformation, the correction step size is taken as... Conversion factor The physical meaning of this is the amount of power adjustment required to compensate for a unit temperature deviation, which is obtained through pre-calibration offline. After the above conversion, the correction step size... The dimensions of the vector are consistent with those of the action vector.

[0155] The revised new candidate action vector The data is resubmitted to the adaptive safety barrier for verification. Using this as control input, the corresponding energy change is derived through an evolutionary three-dimensional thermal model and compared with the adaptive stability criterion at the corresponding time point. If the decay of the energy change meets the criterion requirements, then… If the test succeeds, the released code is executed as the control input for the current cycle; if it still fails, then... Starting anew, the gradient calculation and correction process is repeated, continuing the incremental correction along the new negative gradient direction until the energy change satisfies the current adaptive stability criterion. Since each correction proceeds along the steepest descent direction of the current system's generalized energy, and the correction step size adaptively adjusts with the model's prediction error sequence value, a candidate action vector that meets the safety requirements will be found after a finite number of iterations.

[0156] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely preferred examples and are not intended to limit the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for adaptive adjustment of multi-temperature zone environment of a high-frequency wafer test probe card, characterized in that, Including the following: Step S1: Construct a nominal three-dimensional thermal model: Divide the three-dimensional solid space domain into multiple elements, establish the global thermal conductivity matrix and global thermal capacity matrix of each element, explicitly separate the control input vector and environmental disturbance vector from them and discretize them over time to obtain a nominal three-dimensional thermal model including multiple heat conduction channels. Step S2, Implanting latent parameters and estimating online: Identify the degradation sensitive channels in the nominal three-dimensional thermal model, implant heater efficiency degradation factor and contact thermal resistance drift coefficient as latent parameters into the degradation sensitive channels, form an evolutionary three-dimensional thermal model containing latent parameters, and use ensemble Kalman filtering to update the latent parameter set online based on the measured temperature values ​​of multiple temperature zones; Step S3: Train a safety control strategy network to map the observable temperature state output by the evolutionary three-dimensional thermal model into candidate action vectors. Verify the candidate action vectors using an adaptive safety barrier. Only when the candidate action vector satisfies the adaptive stability criterion based on the control Lyapunov function will it be used as a control input for execution. Otherwise, the candidate action vector will be corrected along the negative direction of the gradient of the control Lyapunov function until the adaptive stability criterion is satisfied.

2. The multi-temperature zone environmental adaptive adjustment method for high-frequency wafer test probe cards according to claim 1, characterized in that: The construction of the nominal three-dimensional thermal model includes: Read all boundary surfaces of the three-dimensional solid space domain and set boundary finite element nodes on the boundary surfaces, with the coordinates of the boundary finite element nodes accurately falling on the boundary surfaces; Starting from the boundary, the internal finite element nodes are generated layer by layer inward, and all finite element nodes are connected into tetrahedral elements based on the Delaunay triangulation criterion. Each element is assigned the same material thermal conductivity, density and specific heat capacity as the physical components of the region where the element is located. The Galerkin weighted residual method is used to calculate the element thermal conductivity matrix and element thermal capacity matrix for each element. The element thermal conductivity matrix and element thermal capacity matrix of all elements are accumulated and assembled into a global thermal conductivity matrix and a global thermal capacity matrix according to the global number of finite element nodes. Each non-zero off-diagonal element in the global thermal conductivity matrix corresponds to a heat conduction channel connecting two finite element nodes. Based on the contact surface position of heating and cooling elements in each temperature zone inside the bearing chuck, a control input matrix is ​​constructed using the area ratio of the contact surface represented by each node as the heat flow distribution coefficient. A disturbance input matrix is ​​constructed based on the convective heat transfer coefficient, surface emissivity, and boundary area of ​​the finite element nodes on the outer surface marked as convective and radiative boundary conditions. Thus, the heat load vector is separated into a linear combination of the control input vector and the environmental disturbance vector. The assembled continuous-time state-space equations are discretized in time with the control period as the step size to obtain the nominal three-dimensional thermal model.

3. The multi-temperature zone environmental adaptive adjustment method for high-frequency wafer test probe cards according to claim 2, characterized in that: The degradation-sensitive channel: For all heat conduction channels in the global thermal conductivity matrix, evaluate the degree of influence of a small change in the thermal conductivity value of each heat conduction channel on the steady-state temperature response of each temperature measurement node on the wafer surface; Preset a corresponding sensitivity threshold; compare the degree of influence with the sensitivity threshold, and define the heat conduction channel whose degree of influence exceeds the sensitivity threshold as a degradation sensitive channel.

4. The multi-temperature zone environmental adaptive adjustment method for high-frequency wafer test probe cards according to claim 2, characterized in that: The implantation latency parameters include: The heater efficiency degradation factor is applied to the construction process of the control input matrix. When the unit heat load vector is integrally calculated, the heat flow injection amount of each temperature zone heating and cooling element is multiplied by the heater efficiency degradation factor of the temperature zone and then integrated. After global assembly, the non-zero elements of the corresponding temperature zone column vector in the input matrix are scaled proportionally by the corresponding degradation factor. The contact thermal resistance drift coefficient is applied to the integral calculation process of the unit thermal conductivity matrix. For the unit containing the contact interface between the probe tip and the wafer pad, the interface thermal conductivity is replaced by the reciprocal of the nominal contact thermal resistance and the reciprocal of the sum of the nominal contact thermal resistance and the contact thermal resistance drift coefficient before being integrated. After global assembly, the thermal conductivity element values ​​of the corresponding node pairs on both sides of the contact interface in the global thermal conductivity matrix are updated. The thermal system matrix and input matrix after the latent parameters are implanted are discretized over time using a control period to obtain an evolutionary three-dimensional thermal model containing latent parameters.

5. The multi-temperature zone environmental adaptive adjustment method for high-frequency wafer test probe cards according to claim 1, characterized in that: The method of using ensemble Kalman filtering to update the latent parameter set online based on measured temperature values ​​from multiple temperature zones includes: Before online operation, heating power steps are applied to each temperature zone one by one, transient temperature response curves are collected, and the initial prior mean and initial prior covariance of the latent parameter vector are obtained through system identification and fitting. Then, an initial set of latent parameters containing multiple members is generated by sampling from the Gaussian distribution determined by the initial prior mean and initial prior covariance. In each control cycle, the latent parameter values ​​of each member in the latent parameter posterior set obtained in the previous cycle and the current measured temperature values ​​of the multi-temperature zones are input into the evolutionary three-dimensional thermal model for forward calculation to obtain a set of temperature prediction values. Using the observation matrix pre-constructed based on the installation positions of physical sensors in the finite element nodes, each temperature prediction value is mapped to the sensor position to obtain the observation prediction value. The state-observation cross-covariance matrix and the observation autocovariance matrix are calculated through ensemble statistics, and then the ensemble Kalman gain is calculated. Based on the deviation between the measured and predicted temperatures in the multi-temperature zones and the Kalman gain of the set, the global node temperature vector and the latent parameter values ​​of each member are updated to obtain the posterior set of latent parameters for the current period.

6. The multi-temperature zone environmental adaptive adjustment method for high-frequency wafer test probe cards according to claim 1, characterized in that: The adaptive stability criterion is dynamically adjusted based on the model prediction error sequence: A fixed-length sliding window is set, and the root mean square value of the prediction residual vector and the root mean square value of the difference between adjacent residuals are calculated for multiple consecutive control cycles within the window. The prediction error sequence of the model is then synthesized by weighting the residuals according to a preset weight. When the model prediction error sequence is in the low error range, it is required to reduce the Lyapunov function value after executing the candidate action vector; When the model prediction error sequence exceeds the low error range, it is required that the decrease in the Lyapunov function value after executing the candidate action vector exceed the maximum prediction deviation mapped by the current model prediction error sequence.

7. The multi-temperature zone environmental adaptive adjustment method for high-frequency wafer test probe cards according to claim 6, characterized in that: When a candidate action vector fails validation, the adaptive safety barrier corrects it in the following way: Using the chain rule, the gradient of the control Lyapunov function with respect to the global node temperature vector is multiplied by the Jacobian matrix of the global node temperature vector prediction value with respect to the candidate action vector to obtain the gradient of the control Lyapunov function with respect to the candidate action vector. Starting with the rejected candidate action vector, correction is performed along the negative direction of the corresponding gradient. The correction step size is determined by the maximum prediction deviation mapped from the current model's prediction error sequence. The revised candidate action vector is resubmitted for verification. If it still fails, the gradient calculation and revision are repeated until a candidate action vector that satisfies the current adaptive stability criterion is found.

8. The multi-temperature zone environmental adaptive adjustment method for a high-frequency wafer test probe card according to claim 1, characterized in that, The security control strategy network is trained offline through reinforcement learning within the evolutionary three-dimensional thermal model. In each virtual control cycle, the security control policy network generates candidate action vectors, submits them to the adaptive security barrier for verification, and executes them as control inputs after passing the verification. The vectors also acquire a composite reward signal consisting of the maximum temperature difference on the wafer surface, the total energy consumption of multiple temperature zones, and the temperature recovery stabilization time. A deterministic policy gradient algorithm is adopted to estimate the action value using the actual reward sequence of samples in the experience replay buffer, and update the network parameters along the gradient direction that maximizes the action value. After training converges, the network parameters are fixed and deployed to a real controller for online temperature control decision-making.

9. The multi-temperature zone environmental adaptive adjustment method for a high-frequency wafer test probe card according to claim 1, characterized in that: The observable temperature state is obtained in the following manner: An observation matrix is ​​pre-constructed based on the installation position of the embedded thermocouples in the finite element nodes. The global node temperature vector output by the evolutionary three-dimensional thermal model is multiplied with the observation matrix to obtain the low-dimensional observable temperature state corresponding one-to-one with the physical sensor.

10. A multi-temperature zone environmental adaptive adjustment and control system for a high-frequency wafer test probe card, characterized in that, include: The model building unit is used to divide the three-dimensional solid space domain occupied by the probe card ceramic substrate, multi-temperature zone heating and cooling elements and carrier chuck into finite element units, establish the global thermal conductivity matrix and global thermal capacity matrix, and separate the control input and environmental disturbance. After time discretization, the nominal three-dimensional thermal model is obtained. The latent parameter estimation unit is used to identify the temperature-sensitive degradation channels in the nominal three-dimensional thermal model, implant heater efficiency degradation factor and contact thermal resistance drift coefficient to construct an evolutionary three-dimensional thermal model with latent parameters, and use ensemble Kalman filtering to update the latent parameter set online based on the measured temperature values ​​of multiple temperature zones. The safety control unit is used to map the observable temperature state into candidate action vectors through a deep neural network, and to verify and correct the candidate action vectors using an adaptive safety barrier, and outputs a control input vector that satisfies the adaptive stability criterion based on the control Lyapunov function.

Images