A method for predicting the lifetime of a photovoltaic inverter

Abstract

The application discloses a kind of photovoltaic inverter life prediction methods.The method includes: obtaining the operating waveform data of photovoltaic inverter and device structure parameters;Based on the transient electrical event contained in operating waveform data, and in combination with device structure parameters, the physical state parameters that characterize the internal characteristics of power device are determined by inversion, and the corrected temperature history is generated;According to the corrected temperature history and physical state parameters, the damage state that characterizes cumulative effect with path dependence is updated;Based on the damage state with path dependence, the critical stability margin of power module is evaluated, and the remaining life of module is determined.The present application associates millisecond-level transient with power redistribution, thermal path degradation and other physical state changes in the device interior through physical closed-loop inverse problem, and uses incremental damage kernel with memory and system-level critical gain criterion, realizes the accurate evaluation of traceability, early warning from transient event to module life.

Description

TECHNICAL FIELD

[0001] The application belongs to the technical field of photovoltaic inverters, and particularly relates to a photovoltaic inverter life prediction method. BACKGROUND

[0002] As the core energy conversion equipment connecting photovoltaic arrays and power grids, the reliability and life of photovoltaic inverters in long-term operation are directly related to the investment return rate, operation and maintenance cost of photovoltaic power stations, and the stability and safety of power grids. Photovoltaic power generation systems are developing towards large capacity, high integration and long service life, which puts forward more stringent requirements for the reliability of inverters. Therefore, developing a life prediction method that is accurate, online and has early warning capability is of great theoretical value and engineering significance for realizing the condition-based maintenance of photovoltaic power stations, optimizing asset management and reducing the levelized cost of electricity (LCOE).

[0003] At present, the life prediction research of photovoltaic inverter power modules mainly focuses on several aspects. One is the offline evaluation method based on mission profile. This method uses long-term statistical data of solar radiation, environmental temperature and other data in a specific region, combined with the efficiency model of the inverter, to estimate the annual temperature cycle load of the power module, and then substitutes it into the classic physical failure model (such as Coffin-Manson model or its improved form) to estimate the life. Another is the statistical analysis method based on field operation data. Through mining a large number of historical operation data and failure data of inverters, machine learning or artificial intelligence algorithms are used to establish a correlation model between operating conditions and device failure probability, so as to realize the probabilistic prediction of remaining life. In addition, there are also researches using offline condition monitoring technology, such as periodic thermal imaging analysis or electrical parameter testing of power modules, to evaluate their health status.

[0004] However, the existing technology still faces several deep-seated technical problems in dealing with the increasingly complex operating conditions and early fault warning needs. These problems mainly come from the observation blind area of the internal microphysical state evolution of the device and the simplified treatment of the multi-physical field coupling degradation mechanism, resulting in insufficient accuracy and timeliness of the prediction results.

[0005] Firstly, the existing methods generally ignore the rich degradation information contained in the millisecond-level electrical-thermal transient events. For example, under the conditions of grid-connected disturbance, reactive power support switching or MPPT micro-disturbance, the electrical transient response of the device (such as on-state voltage drop, switch ringing, etc.) will change slightly due to early damage such as internal bond line lift-off. Traditional analysis methods based on low sampling rate or statistical mean cannot capture and interpret these "transient fingerprints", thus missing the perception of early physical defects, and even more cannot include the power redistribution phenomenon between parallel chips caused thereby as a key accelerator of aging into the life assessment model. Secondly, the classic damage accumulation model essentially regards thermal cycling as an independent "amplitude-frequency" event, lacking the ability to remember the loading history and path. This makes it difficult for them to accurately depict the "thermal ratchet" one-way cumulative damage effect produced by the combined action of temperature rise / drop rate, high / low temperature residence time in the stress singular zone at the edge of the chip. Finally, the existing evaluation system lacks a critical criterion from the micro device degradation to the macro system instability. The electro-thermal positive feedback of the parallel bridge arm is the key factor leading to the catastrophic failure of the module, but the traditional method cannot establish a deterministic physical correlation between the local damage state (such as slow change of on-state resistance) and the system dynamic bifurcation point (i.e. thermal runaway critical point), leading to inaccurate prediction of the critical safety margin of the module and lagging warning. SUMMARY

[0006] The present application provides a photovoltaic inverter life prediction method.

[0007] Technical scheme: A photovoltaic inverter life prediction method, comprising:

[0008] Obtain the operating waveform data and device structure parameters of the photovoltaic inverter, and the operating waveform data contains transient electrical events;

[0009] Based on the transient electrical events, and combined with the device structure parameters, the physical state parameters characterizing the internal characteristics of the power device are determined by inversion, and the corrected temperature history is generated;

[0010] According to the corrected temperature history and physical state parameters, update the damage state with path dependence characterizing the cumulative effect;

[0011] Based on the damage state with path dependence, evaluate the critical stability margin of the power module, and determine the remaining life of the module.

[0012] Beneficial effects: The application realizes online and quantitative sensing of early degradation precursors such as bond wire lift-off, and captures the key aging accelerator-chip interpower imbalance caused thereby, filling the observation blind area of traditional methods; makes damage accumulation no longer a simple superposition of isolated thermal cycles, but a more physically realistic process with directionality and historical dependence, improving the accuracy of life prediction under complex and variable working conditions; provides a critical point prediction method with clear physical meaning from device degradation to system bifurcation, improves the judgment of life termination from an empirical threshold to the level of system stability analysis, and enhances the reliability and advance of early warning. BRIEF DESCRIPTION OF DRAWINGS

[0013] Figure 1 The technical scheme of the application is a whole flowchart.

[0014] Figure 2 The flowchart of the embodiment of the application for determining physical state parameters through inversion.

[0015] Figure 3 The flowchart of the embodiment of the application for updating the damage state with path dependence representing cumulative effects.

[0016] Figure 4 The flowchart of the embodiment of the application for evaluating the critical stability margin of the power module and determining the remaining life of the module. DETAILED DESCRIPTION

[0017] In order to enable those skilled in the art to better understand the present application, the technical solutions in the embodiments of the present application will be described clearly and completely below with reference to the accompanying drawings of the embodiments of the present application. Obviously, the described embodiments are only a part of the embodiments of the present application, rather than all the embodiments. Based on the embodiments in the present application, all other embodiments obtained by those skilled in the art without creative efforts should fall within the protection scope of the present application.

[0018] It should be noted that the terms first, second, etc. in the specification and the above drawings are used to distinguish similar objects, and do not necessarily indicate a specific order or sequence. It should be understood that the data thus used can be interchanged under appropriate circumstances, so that the embodiments of the application described herein can be implemented in an order other than that illustrated or described herein. In addition, the terms include and have and any variations thereof are intended to cover non-exclusive inclusion, for example, a process, method, system, product or device including a series of steps or units does not necessarily limit to the clearly listed steps or units, but can include other steps or units not clearly listed or inherent to the process, method, product or device.

[0019] Embodiment one

[0020] As Figure 1 shown in the drawings, according to an embodiment of the present application, a photovoltaic inverter life prediction method is provided, and specifically, the method comprises the following steps:

[0021] In step S101, the operating waveform data and device structure parameters of the photovoltaic inverter are obtained.

[0022] In this embodiment, the operating waveform data can be high-frequency data collected from the controller or monitoring unit of the photovoltaic inverter, and exemplary includes the output current waveform of the inverter, the DC bus voltage, the gate drive voltage signal, and the switching timing, etc. The device structure parameters can be prior information provided by the power module design or factory, such as the size of the power semiconductor chip, the topology layout of the parallel chip, the bonding wire layout, the thickness of each layer of solder, and the material thermal physical parameters of the lead and substrate, etc.

[0023] Further, after obtaining the original data, data preparation and model initialization can be performed first. Preferably, the process can include uniform alignment of the time reference of the collected operating waveform data and elimination of abnormal data segments. Further, through a short window power loss estimation algorithm, combined with a multi-node thermal resistance and thermal capacity network model, the junction temperature of the power device is preliminarily calculated to obtain a high-resolution temperature trajectory preliminary value, which is denoted as Tj_highres here. Further, the time stamp sequence of transient electrical events is identified and extracted from the operating data, such as power steps, reactive power adjustments, perturbations in the maximum power point tracking (MPPT) process, etc. These event windows will serve as trigger signals for subsequent analysis.

[0024] As a preferred scheme, the multi-node thermal network model can be self-corrected online once after initialization using the event window data identified above. Specifically, by comparing the deviation of the model response during the event from the actual temperature trajectory Tj_highres, the thermal resistance and thermal capacity parameters in the model are fine-tuned in a small range, so as to obtain a basic thermal network model that is closer to the current actual working condition and can be called by subsequent steps, which is denoted as callable thermal network model_base here. Through the self-correction mechanism of this scheme, the accuracy of the baseline of all subsequent calculations can be improved.

[0025] In step S102, based on the transient electrical events contained in the operating waveform data and combined with the device structure parameters, the physical state parameters characterizing the internal characteristics of the power device are inversely determined, and a corrected temperature history is generated.

[0026] The transient electrical event is an electrical-thermal transient process identified in step S101, caused by a working condition change such as a power step, and lasting for milliseconds. The physical state parameter is used to quantify the microscopic characteristic changes in the power module caused by aging and degradation. In this embodiment, the physical state parameter preferably includes two parts: a power redistribution matrix M(t) for characterizing the uneven heating power between parallel chips, and a thermal network parameter perturbation delta_theta(t) for characterizing the degradation of the device thermal path. In a specific implementation, by constructing an inverse problem, the two internal state parameters can be deduced from the externally measurable electrical transient response. Further, after obtaining M(t) and delta_theta(t), they are injected into the callable thermal network model_base to correct the initial estimate of the temperature trajectory obtained in step S101, thereby generating a more accurate corrected temperature history considering the actual unevenness and degradation effects, denoted as Tj_corr.

[0027] In step S103, the damage state with path dependence characterizing the cumulative effect is updated according to the corrected temperature history and the physical state parameter.

[0028] In this embodiment, the traditional life model (such as the Coffin-Manson model) regards each thermal cycle as an independent event, ignoring the influence of the historical path. To solve this problem, this embodiment introduces an internal state variable denoted as xi_edge, which is used to characterize the one-way cumulative damage with memory caused by thermal ratcheting effect and other stress singularities at the chip edge. Specifically, Tj_corr obtained in step S102 is divided into a series of discrete thermal cycles; further, for each thermal cycle, an incremental function with path dependence is used to perform a recursive update on xi_edge. As a preferred, the input of the incremental function not only includes the temperature characteristics of the current thermal cycle (such as uniform temperature, swing, temperature rise and fall slope, dwell time), but also couples the local stress term determined by M(t) and delta_theta(t) inverted in step S102, so that the cumulative process of damage is closely related to the real physical state change.

[0029] In step S104, based on the damage state with path dependence, the critical stability margin of the power module is evaluated, and the remaining life of the module is determined.

[0030] In this embodiment, as xi_edge accumulates, the macroscopic electrical parameters (such as on-resistance R_on) and thermal parameters (such as thermal resistance R_th) of the power devices drift, exacerbating the electro-thermal positive feedback between parallel chips, which can lead to thermal runaway in severe cases. This step aims to assess this risk from the perspective of system stability. Specifically, the damage state represented by xi_edge is mapped to time-varying terms of R_on and R_th, and a closed-loop model describing the electro-thermal coupling effect of the parallel chip system is constructed. Further, by calculating the closed-loop small-signal gain matrix K_loop and its spectral radius rho of this model, the critical stability margin, denoted as margin, can be obtained, where margin = 1 - rho. A margin approaching 0 means that the system is approaching the imbalance critical point. Further, by comprehensively considering the cumulative damage of each chip (derived from xi_edge) and the critical stability margin of the entire module, the remaining lifetime RL_module of the module is determined.

[0031] Example 2

[0032] like Figure 2 As shown, this embodiment, based on Embodiment 1, elaborates on the process of determining the physical state parameters through inversion in step S102. In a specific embodiment, this process may include:

[0033] Step S201: For any transient electrical event, extract a set of transient electrical fingerprints from its corresponding operating waveform data. The transient electrical fingerprints are used to quantify the electrical transient response during the event. Further, the transient electrical fingerprints are used as inputs for subsequent inversion operations to determine physical state parameters.

[0034] Transient electrical fingerprints are a set of features used to characterize subtle changes in electrical behavior during transient processes. In a preferred embodiment, this fingerprint may include, but is not limited to, the following four items:

[0035] On-state voltage drop short window offset delta_Vce_on: This represents the offset of the average on-state voltage drop of the device within a short time window during a transient event relative to the steady-state value before the event. It is more sensitive to early degradation caused by changes in current path resistance, such as bond wire lift-off.

[0036] Rising edge overshoot amplitude V_ov: This represents the voltage overshoot peak value that occurs at the rising edge of the collector-emitter voltage waveform during device turn-on. It is related to the stray inductance and current change rate of the commutation circuit, and the stray inductance may change due to bonding wire deformation or desoldering.

[0037] Turn-off ringing decay constant zeta: indicates the speed of voltage or current waveform ringing decay during device turn-off process, its variation can reflect the variation of device parasitic capacitance and inductance, and indirectly relates to the aging of physical structure.

[0038] Transient thermal resistance short window estimation delta_Zth: an estimation of transient thermal resistance by calculating the ratio of instantaneous power and instantaneous junction temperature variation within a transient event window, and calculating its deviation from the nominal value, reflecting the degradation of thermal path (such as solder layer).

[0039] In some embodiments, in order to facilitate subsequent unified processing, the extracted fingerprints can be normalized and weighted fused. Exemplarily, a normalized fingerprint weight is obtained by a linear combination w_norm = a0 + k1 x delta_Vce_on + k2 x V_ov + k3 x zeta + k4 x delta_Zth; wherein w_norm is the normalized weight; a0 is the bias coefficient; k1, k2, k3, k4 are the weighting coefficients of each fingerprint, which values can be calibrated by experiments according to the sensitivity of different device types and degradation modes.

[0040] Step S202, construct and solve an inverse problem with transient electrical fingerprints as input, the solution as physical state parameters, including at least: power redistribution matrix for characterizing the uneven heating power between parallel chips; or thermal network parameter perturbation for characterizing the degradation of thermal path. In some embodiments, in addition to the transient electrical fingerprints, some other data are also included.

[0041] In this embodiment, the inverse problem refers to the process of inversely inferring the internal state parameters that cannot be directly measured from the external observations (i.e. transient electrical fingerprints and temperature response) of the system. This corresponds to the positive problem of calculating the external response from the internal parameters.

[0042] As a preferred implementation, before directly solving the inverse problem, the method further includes an initialization step, i.e. generating power redistribution matrix initial value and thermal network parameter perturbation initial value directly from transient electrical fingerprints through linear prior mapping; further, taking the power redistribution matrix initial value and the thermal network parameter perturbation initial value as the starting point of iterative solution of the inverse problem.

[0043] Specifically, the approximate linear relationship between the transient fingerprint and the internal state parameters can be established in advance through experiments or simulations. For example, m_i(t) = a0_i + a1_i x delta_Vce_on + a2_i x V_ov + a3_i x zeta + a4_i x delta_Zth; delta_theta_Rj = b0 + b1 x delta_Zth + b2 x zeta; where m_i(t) is the initial value of the i-th chip in the power redistribution matrix M_init(t); delta_theta_Rj is the initial value of the component related to the junction shell thermal resistance in the thermal network parameter perturbation delta_theta_init(t); a0_i to a4_i and b0 to b2 are mapping coefficient matrices, which are obtained through offline calibration. This prior mapping is used to generate an iteration starting point with high quality, which can improve the convergence speed and robustness of the inverse problem solution and avoid falling into a local optimal solution.

[0044] According to a further improvement of the present application, the solution of the inverse problem is realized by minimizing a joint cost function on this basis. The joint cost function at least includes: a temperature difference term representing the difference between the simulated temperature and the high-resolution temperature trajectory; and a total variation regularization term for constraining the physical reasonableness of the power redistribution matrix or the thermal network parameter perturbation.

[0045] Specifically, the joint cost function J can be represented as: J = wT x ||Tj_highres - Tsim(M, delta_theta)||2 2 + lambda1 x TV(M) + lambda2 x TV(delta_theta); where J is the total cost to be minimized; wT is the weight of the temperature difference term, which can be in the range of 0.5 to 0.9; Tj_highres is the high-resolution temperature trajectory observation value obtained in Embodiment One; Tsim(M, delta_theta) is the temperature trajectory simulated by calling the thermal network model_base using the power redistribution matrix M and the thermal network parameter perturbation delta_theta of the current iteration; ||...||2 2 represents the square of the L2 norm, i.e., the sum of squared errors; lambda1 and lambda2 are regularization term coefficients, which can be in the range of 1e-4 to 1e-2, for example, to balance the fitting accuracy and the smoothness of the solution; TV(M) and TV(delta_theta) are total variation regularization terms, which are used to punish the drastic changes of the solution, making it more physically reasonable, i.e., the degradation process is usually gradual rather than sudden.

[0046] Optionally, the minimization problem can be solved by an iterative algorithm of alternating minimization. For example, fix delta_theta, update M by one step of gradient descent, then fix the updated M, update delta_theta by one step of gradient descent, and so on until convergence. Further, after each iteration update, the solution also needs to be projected back into its physical constraint set, for example, the sum of elements of M must be 0 (total power conservation), the elements of delta_theta must be non-negative, etc.

[0047] Step S203, further, after solving the inverse problem, evaluate the sensitivity of the problem solution to generate an observability index; when the observability index is lower than a preset threshold, generate an adaptive perturbation plan, which is used to actively adjust the operating condition or data sampling strategy in the subsequent evaluation period to improve the observability.

[0048] Since not all transient events can provide sufficient information for the inverse problem. If the temperature change caused by an event is not sensitive to the internal state parameters M and delta_theta, the confidence of the inversion result is low. Specifically, the observability index J_obs can be obtained by calculating the norm of the Jacobian matrix of the simulated temperature with respect to the state parameters: ; Where J_obs is the observability index; ||...||_F represents the Frobenius norm of the matrix; Tsim / M and Tsim / delta_theta are the sensitivity (Jacobian) matrices of the simulated temperature with respect to M and delta_theta respectively; gamma1 and gamma2 are weighting coefficients, for example, both can be taken as 0.5.

[0049] As a preferred embodiment, when the calculated J_obs is lower than a preset threshold J_min (for example, J_min = 0.1 x M_scale, where M_scale is the typical scale of M), it is considered that the data quality of the current event window is insufficient, and the system will generate an adaptive perturbation plan. The plan can impose a small, grid-compliant active disturbance on the operating condition through the controller of the inverter in the next suitable evaluation period (for example, short-time adjustment within 1% of the rated power), to excite a more sensitive system response to the internal state parameters, so as to obtain higher quality data in the subsequent evaluation, forming an adaptive closed loop of risk-observation-modeling.

[0050] Embodiment three

[0051] As Figure 3As shown, the embodiment is based on the embodiments one and two, and elaborates the process of updating the damage state with path dependence and cumulative effect in step S103 in detail. The process can specifically include:

[0052] In step S301, the corrected temperature history is divided into a series of discrete thermal cycles; further, for each thermal cycle, an intrinsic state variable for representing the cumulative effect of the chip edge singular zone is updated recursively.

[0053] In the embodiment, the division of thermal cycles can be performed according to the method of zero crossing and residence threshold. Specifically, the potential cycle inflection point is identified by detecting the sign change (zero crossing point) of the time derivative (i.e. the temperature rising and falling rate) of the corrected temperature history Tj_corr, and a threshold of minimum temperature difference and minimum residence time is set to filter out the pseudo cycles caused by high frequency noise, so as to obtain a series of effective thermal cycles with clear starting point, peak point and end point. The intrinsic state variable, denoted as xi_edge here, is a scalar or vector, and its physical meaning is to represent the cumulative state of irreversible micro plastic strain or micro crack in the stress concentration area such as the chip edge and the solder layer due to the mismatch of material properties (such as the difference in thermal expansion coefficient between the chip and the substrate).

[0054] The recursive update of xi_edge reflects the memory of damage, that is, the damage increment caused by the current cycle is based on the state at the end of the last cycle.

[0055] In step S302, the recursive update of the intrinsic state variable is performed according to a path-dependent increment function. In a preferred embodiment, the inputs of the increment function include the temperature rising and falling slope and the residence time extracted from the thermal cycle, and further include a local stress term depending on the physical state parameters.

[0056] In this embodiment, the recursive update of xi_edge for any thermal cycle with index k can be expressed as xi_edge(k+1) = xi_edge(k) + Delta_xi_edge(k); where xi_edge(k) is the state at the beginning of the kth cycle, xi_edge(k+1) is the state at the end of the cycle, and Delta_xi_edge(k) is the state increment generated by the current cycle. The state increment Delta_xi_edge(k) is calculated by an increment function g1, which can have the specific form: Delta_xi_edge(k) = p0 + p1 x T_mean(k) + p2 x Delta_T(k) + p3 x |dT_dt(k)| - S_relax(k) + p5 x sigma_local(k); where T_mean(k), Delta_T(k), and dT_dt(k) are the mean junction temperature, temperature swing, and absolute value of temperature ramping rate extracted from the kth thermal cycle, respectively; S_relax(k) is a relaxation term, which can be expressed as S_relax(k) = beta_relax x xi_edge(k), for example, to simulate the physical effects of stress relaxation or partial damage recovery that can occur during continuous low temperature or long dwell, with beta_relax being the relaxation coefficient; p0 to p5 are model coefficients calibrated by experiments; and sigma_local(k) is the local stress term. The local stress term sigma_local(k) couples the effects of the physical state parameters M(t) and delta_theta(t) inverted in Embodiment Two into the damage model, and has the calculation formula: sigma_local(k) = s0 + s1 x Delta_T(k) + s2 x T_mean(k) + s3 x |dT_dt(k)| + s4 x t_dwell(k) + s5 x phi(M, delta_theta); where s0 to s5 are coefficients; t_dwell(k) is the dwell time of the kth thermal cycle in the high and low temperature intervals; and phi(M, delta_theta) is the core coupling function, which represents the local additional heat flow or temperature gradient introduced by the power redistribution M and the thermal parameter perturbation delta_theta, resulting in local stress concentration. By introducing this term, when the chip has wire debonding (affecting M) or solder layer voids (affecting delta_theta), even if the macro temperature cycle characteristics remain unchanged, the damage increment generated will also increase accordingly, thus more truly reflecting the physical degradation process.

[0057] Further, the path-dependent incremental function differentiates the heating and cooling phases of the thermal cycle. In the present embodiment, this is achieved by using different gains for the heating and cooling slope terms in the incremental function. Specifically, in the calculation of Delta_xi_edge(k), the coefficient p3 is adjusted according to the sign of dT_dt(k): when dT_dt(k) > 0 (heating phase), the heating gain p3_up = p3 x (1 + r_up) is used; when dT_dt(k) < 0 (cooling phase), the cooling gain p3_down = p3 x (1 + r_down) is used. Here, r_up and r_down are asymmetric gain coefficients, for example, r_up can take 0.2 to 0.6, and r_down can take 0.1 to 0.4. It is found through analysis that the plastic deformation of the material during the heating process and the damage accumulation rate caused by creep or reverse yield during the cooling process are often different. By introducing asymmetric gain, the present method can more accurately simulate the thermal ratchet effect, i.e. a net plastic deformation accumulation after each cycle.

[0058] Step S303, from the transient electrical fingerprint, identify and generate a lift-off event list; and when performing the recursive update of the intrinsic state variable, if it is determined that the current thermal cycle is in the event neighborhood in the lift-off event list, adjust the damage increment calculated by the path-dependent incremental function by a magnification weight to simulate the accelerated degradation effect caused by the event.

[0059] Specifically, in step S201 of Embodiment Two, a threshold criterion based on the transient electrical fingerprint can be set. When the deviation of delta_Vce_on exceeds a certain threshold, it is determined that a lift-off preliminary event has occurred, and the time stamp and an initial weight w_lo0 derived from w_norm are recorded in the lift-off event list.

[0060] In the embodiment, when performing the recursive update in step S302, Delta_xi_edge(k) in the neighborhood of the lift-off event timestamp (e.g., several thermal cycles before and after the event) is amplified with a weighting. Exemplarily, the adjustment is performed in a gated manner: Delta_xi_edge_gated(k) = Delta_xi_edge(k) x w_event(k); wherein the amplification weight w_event(k) can be calculated according to w_lo0, for example, w_event(k) = alpha_e + k_e x w_lo0(k), alpha_e and k_e being coefficients. Through the technical solution of the embodiment, the problem that the occurrence of lift-off deteriorates the current distribution and causes local hot spots and accelerates the degradation process of the region in subsequent thermal cycles is solved.

[0061] In step S304, the state increment calculated by the incremental function is mapped to the incremental damage of the thermal cycle.

[0062] In the embodiment, the mapping process can include applying a preset damage threshold to the state increment, and only when the state increment exceeds the threshold, it is determined to contribute to the damage; and the part exceeding the damage threshold is processed by a nonlinear mapping function to determine the incremental damage.

[0063] Specifically, the calculation of the incremental damage DeltaD_edge(k) of the kth cycle is divided into two steps. The damage threshold S_th is applied to the state increment Delta_xi_edge_gated(k) to obtain the effective state increment Delta_xi_pos(k) = max(Delta_xi_edge_gated(k) - S_th, 0); the physical meaning of S_th is the fatigue limit of the material, i.e., a cycle with a small amplitude will not cause cumulative damage. Further, the effective state increment is mapped to damage by a nonlinear function g2: DeltaD_edge(k) = q0 + q1 x Delta_xi_pos(k) + q2 x (Delta_xi_pos(k)) gamma ; wherein q0, q1, q2 are mapping coefficients, and gamma is a nonlinear cumulative effect index, usually taking a value greater than 1 (e.g., 1.2 to 2.5), indicating that the cumulative rate of damage will accelerate with the degradation of the state.

[0064] In step S305, the damage increment determined by the path-dependent incremental function is calculated, and the sensitivity of at least one of the average junction temperature, the temperature swing, the temperature rise and fall slope, or the residence time of the thermal cycle is calculated; further, based on the calculation result of the sensitivity, a singular point sensitivity report is generated.

[0065] Specifically, by calculating the partial derivatives of the incremental damage DeltaD_edge against T_mean, Delta_T, |dT_dt|, t_dwell, the respective sensitivities can be obtained, for example, S_Tmean = DeltaD_edge / Delta_T_mean. DeltaD / T_mean. By summarizing these sensitivity values and their ranking, a singularity sensitivity report can be generated. This report can reveal which working condition (e.g., long-term high-temperature load or frequent power fluctuation) is the main factor leading to the aging of the current device under its specific operating profile, providing direct data support for formulating targeted operation and maintenance strategies or optimizing the operating curve.

[0066] As an optional solution, in order to ensure the long-term numerical stability of the algorithm, an upper limit clipping can be applied to the state variable xi_edge and the incremental damage DeltaD_edge after each recursive update, in order to prevent numerical explosion caused by abnormal working conditions or model divergence. Preferably, the values of the numerous model parameters involved in the present embodiment (such as the p series, s series, q series coefficients, gamma, S_th, etc.) can be determined through a two-stage calibration process: using specific bench thermal cycle test data to determine the threshold-related and asymmetry-related parameters; further, using a larger range of accelerated aging test data to optimize the coefficients in the linear region to minimize the error between the simulated damage curve and the measured degradation index (such as crack length, incremental on-state resistance, etc.).

[0067] Embodiment Four

[0068] As shown in Figure 4 the present embodiment further elaborates on the process of evaluating the critical stability margin of the power module in step S104 and determining the remaining life of the module based on embodiments one, two and three. The process can specifically include:

[0069] Step S401, updating the time-varying equivalent electrical parameters and time-varying equivalent thermal parameters of the power device according to the damage state with path dependence; as a preferred, the time-varying equivalent electrical parameters at least include on-state resistance; as a preferred, the time-varying equivalent thermal parameters at least include thermal resistance.

[0070] Specifically, the intrinsic state variable xi_edge_gate(i) (where i is the chip index) obtained in Example Three is updated by a pre-calibrated mapping relationship to update the on-resistance R_on(i) and the junction housing thermal resistance R_th(i) of each chip: R_on(i) = R_on0(i) + c_r × xi_edge_gate(i) R_th(i) = R_th0(i) + c_h × xi_edge_gate(i); where R_on0(i) and R_th0(i) are initial (in a healthy state) parameter values of the chip; c_r and c_h are mapping coefficients obtained by regression of device aging experiment data. Through the technical solution of this embodiment, a bridge from microscopic cumulative damage to macroscopic device parameter change is established, which is the basis for subsequent system-level stability analysis.

[0071] In step S402, a model describing the electrical-thermal coupling effect of the parallel chip system is constructed based on the updated time-varying equivalent electrical parameter and the time-varying equivalent thermal parameter, and a closed-loop gain matrix is determined therefrom. Further, the spectral radius of the closed-loop gain matrix is calculated, and the spectral radius is used to determine the critical stability margin.

[0072] By way of example, in the parallel chip system, there is a potential positive feedback loop: a slight increase in the temperature T_chip of any chip → an increase in its R_on → if it is a parallel voltage-sharing structure, its shared current I_share may decrease, but if it is a current-sharing structure or under certain dynamics, it may also cause its loss P_loss to change → the change in P_loss in turn further affects T_chip. When the gain of this loop is greater than 1, a small disturbance will be continuously amplified, eventually leading to thermal runaway. In this embodiment, the gain of this feedback loop can be quantified. In a specific implementation, the closed-loop gain matrix K_loop can be obtained by the chain rule: K_loop = (dI_share / dT_chip) × (dT_chip / dP_loss) × (dP_loss / dR_on) × (dR_on / dI_share). For a system containing N parallel chips, K_loop is an N x N matrix, where each item K_loop(i,j) represents the gain of the current change of chip j on the current change of chip i. The calculation method of each partial derivative in the matrix is as follows:

[0073] dI_share / dT_chip: calculated by applying a small disturbance to the parallel current sharing model. For example, increase the junction temperature T_chip(j) of chip j by a small disturbance dT (such as 1K), recalculate the current I_share of all chips, and thus obtain dI_share(i) / dT_chip(j);

[0074] dT_chip / dP_loss: approximately equal to the transient thermal resistance matrix of the chip Z_th, i.e. dT_chip(i) / dP_loss(j) = Z_th(i,j);

[0075] dP_loss / dR_on: approximately a diagonal matrix with diagonal elements (i,i) = I_rms(i) 2 , I_rms(i) being the current root mean square flowing through chip i;

[0076] dR_on / dI_share: this term is the key electro-thermal coupling term, and is calculated by chain rule as dR_on(i) / dI_share(j) = (dR_on(i) / dT_chip(i)) x (dT_chip(i) / dI_share(j)).

[0077] Further, after obtaining the K_loop matrix, its spectral radius rho(K_loop) is calculated, which is the maximum of the modulus of its eigenvalues. The physical meaning of the spectral radius is the gain of the system in the least stable mode. The critical stability margin margin is defined as margin = 1 - rho(K_loop). When margin starts to decrease from 1 and approaches to 0, it indicates that the system is approaching to the critical instability point.

[0078] Step S403, after determining the critical stability margin, a hysteresis criterion with different entering threshold and exiting threshold is adopted to process the critical stability margin, to generate a robust critical instability point flag.

[0079] In this embodiment, to avoid false alarms caused by the instantaneous jitter of margin due to calculation noise or temporary operating condition fluctuations, a hysteresis criterion can be used. Specifically, a lower entering threshold margin_min_enter (e.g. 0.05) and a higher exiting threshold margin_min_exit (e.g. 0.08) are set. Only when margin falls below margin_min_enter, a critical warning flag is generated; and once in the warning state, only when margin rises above margin_min_exit, the flag is removed. As an alternative or supplementary solution, a 5-point median filter can also be applied to the calculated margin sequence to further suppress the interference of instantaneous noise.

[0080] Step S404, the critical stability margin is monitored, and when it narrows to a preset risk threshold, a sampling and window scheduling suggestion is generated; preferably, the sampling and window scheduling suggestion is fed back to adjust the data sampling density or window selection strategy for transient electrical events.

[0081] By the technical scheme of the embodiment, a final closed loop of risk-observation-modeling of the application is formed. When the critical early warning flag generated in step S403 is triggered, it indicates that the system health is declining, and at this time, the demand for more accurate observation of the internal state of the system is also increased. Therefore, the system automatically generates a scheduling suggestion, which is fed back to the data acquisition and processing module of Embodiment Two. As a preferred technical scheme, the suggestion can include: increasing the sampling frequency of waveform data when the next transient event occurs; or, extending the window length of event analysis; or, directly triggering the adaptive perturbation plan described in Embodiment Two to actively create a more observable incentive.

[0082] In step S405, based on the accumulation of the incremental damage, the remaining life of each chip is determined; and in combination with the weakest value in the remaining life and the critical stability margin, the remaining life of the module is comprehensively determined. In other words, in combination with the weakest value of each chip and the critical stability margin, the comprehensive remaining life is determined.

[0083] In a specific embodiment, by accumulating the incremental damage DeltaD_edge_i(k) of each chip calculated in Embodiment Three by time, the cumulative damage D_chip(i) = sum_k(DeltaD_edge_i(k)) of each chip (each chip, each chip, and each chip) is obtained. The remaining life of each chip can be approximately obtained by the linear damage accumulation rule (Miner's rule) RL_chip(i) = L_ref(i) - D_chip(i), where L_ref(i) is the reference life.

[0084] Further, the final remaining life RL_module of the module is determined by the life of the weakest chip and the stability margin of the system level, for example: RL_module = min_i(RL_chip(i)) × (1 + k_m × margin), where min_i(RL_chip(i)) represents the short board of the barrel, and k_m is a coefficient. The actual life of the module depends not only on the chip that fails first, but also on the stability of the entire system. A system with a small margin may fail prematurely due to a slight disturbance even if the damage to the weakest chip is not serious. Through the technical scheme of the embodiment, a more comprehensive and reliable module remaining life prediction result can be obtained than considering only the chip life or the system stability.

[0085] Embodiment Five

[0086] The embodiment provides a specific numerical calculation case of a photovoltaic inverter life prediction method, which is used for illustratively explaining the complete process and calculation details of the method in embodiments one to four.

[0087] In the embodiment, it is assumed that a power module of a photovoltaic inverter contains two parallel IGBT chips, denoted as chip 1 and chip 2. At t=0, the inverter experiences a power step event caused by the rapid dispersion of a cloud layer. The initial conditions are set as follows:

[0088] The ambient temperature is 25°C, and the cooling system is working normally; both chips are in a healthy state, and the initial intrinsic state variables xi_edge(0) of both chips are 0, and the cumulative damage D_chip(0) is also 0. The nominal on-resistance R_on0 is 2.5 mΩ, and the junction shell thermal resistance R_th0 is 0.1 K / W; a calibrated callable thermal network model_base is ready; a high-frequency sampling system records the complete running waveform data_raw before and after the power step event.

[0089] Further, data preparation and temperature preliminary estimation, including: according to the collected running waveform data_raw, a short window power loss estimation is performed, and the callable thermal network model_base is convolved to obtain the high-resolution temperature trajectory preliminary estimation Tj_highres of the two chips. For the sake of simplicity, the temperature change of chip 1 is focused on: the junction temperature is stabilized at 60°C before the event, and after the power step, a short oscillation is experienced, the peak value reaches 102°C, and finally it is stabilized at 95°C. The timestamp of the power step event is recorded.

[0090] Further, event-driven physical closed-loop inversion, including the following steps:

[0091] 1. Transient fingerprint extraction: in the time window of the power step event, a set of transient fingerprints: delta_Vce_on = +6 mV; V_ov = +4 V; zeta = -0.025; delta_Zth = +0.012 K / W are extracted from the electrical waveform of chip 1. The non-zero values of these fingerprints imply that chip 1 may have slight bond wire degradation and thermal path degradation.

[0092] 2. Initial value generation for inverse problem: Using the linear prior mapping described in Example 2, generate initial values for inverse problem iterations based on the fingerprints described above. Illustratively, assume the calibrated mapping is: m_1(t) = 0.008 x delta_Vce_on(mV) - 0.001 x V_ov(V); delta_theta_Rj1 = 0.8 x delta_Zth. Substitute the fingerprint values to get the initial value of the power redistribution matrix M_init for chip 1: m_1_init = 0.0086 - 0.0014 = 0.044; further, due to the total power conservation constraint sum_i(m_i) = 0, the initial value of the power redistribution matrix for chip 2 is m_2_init = -0.044; further, the initial value of the thermal network parameter perturbation delta_theta_init for chip 1 is delta_theta_Rj1_init = 0.8 x 0.012 = 0.0096 K / W.

[0093] 3. Iterative solution and convergence: Construct the joint cost function J of Example 2, set the weights wT= 0.7, λ1= 0.001, λ2= 0.0005; further, starting from the initial values described above, use the alternating minimization algorithm to iteratively solve; further, after, e.g., 30 iterations, the algorithm converges, and the final physical state parameters are obtained: the power redistribution matrix M(t) = [0.06, -0.06]. This indicates that the actual heat generation of chip 1 is 6% higher than the nominal value, while chip 2 is 6% lower; the thermal network parameter perturbation delta_theta(t) = [0.015 K / W, 0.002 K / W]. This indicates that the junction-to-case thermal resistance of chip 1 has actually increased by 0.015 K / W.

[0094] 4. Temperature history correction: Inject the converged M(t) and delta_theta(t) into the callable thermal network model_base to recalculate the junction temperature; obtain the corrected temperature history Tj_corr. For chip 1, the peak temperature is corrected to 105°C, which is 3°C higher than the initial estimate. This difference is caused by the power imbalance and thermal resistance degradation that was not captured by the initial model.

[0095] Further, the path-dependent damage accumulation includes the following steps:

[0096] 1. Thermal cycle statistics extraction: From the Tj_corr curve of die 1, a complete thermal cycle induced by this power step event is extracted. The statistics are: mean junction temperature T_mean = 82.5°C; temperature swing Delta_T = 45K (from 60°C to 105°C); temperature ramp-up rate dT_dt = +20 K / s; high and low temperature dwell time t_dwell is negligible in this fast event.

[0097] 2. State delta calculation: State delta Delta_xi_edge is calculated according to the delta function g1 of embodiment three. Exemplarily, assume the model parameters are calibrated as: pi = le-3, p2 = 1.5e-3, p3 = 2e-4, p5 = le-4, r_up = 0.3, si = 0.02, s5 = 0.1; calculate local stress term sigma_local. phi(M, delta_theta) is simplified as the term related to m_1 and delta_theta_Rj1, and phi ~ 1.5 is calculated, then sigma_local ~ si x Delta_T + s5 x phi = 0.0245 + 0.11.5 = 1.05; further, since it is a temperature ramp-up phase, the gain of p3 is p3_up = 2e-4 x (1 + 0.3) = 2.6e-4; further, calculate state delta: Delta_xi_edge(1) ~ pi x T_mean + p2 x Delta_T + p3_up x |dT_dt| + p5 x sigma_local = le-382.5 + 1.5e-345 + 2.6e-420 + le-41.05 ~ 0.0825 + 0.0675 + 0.0052 + 0.0001 = 0.1553.

[0098] 3. Delta damage mapping: Map the state delta to damage according to the mapping function g2 of embodiment three. Exemplarily, assume the parameters are S_th = 0.01, qi = 5e-5, q2 = le-4, gamma = 1.5; calculate effective state delta Delta_xi_pos = max(0.1553 - 0.01, 0) = 0.1453; further, calculate delta damage DeltaD_edge(1) = 5e-5 x 0.1453 + le-4 x (0.1453) 1.5 ~ 7.27e-6 + 5.53e-6 ~ 1.28e-5.

[0099] 4. State and damage update: update the state variable and accumulated damage of chip 1: xi_edge(1) = xi_edge(0) + Delta_xi_edge(1) = 0 + 0.1553 = 0.1553; D_chip(1) = D_chip(0) + DeltaD_edge(1) = 0 + 1.28e-5 = 1.28e-5. (The damage calculation process of chip 2 is similar, but because its power sharing is reduced, the damage increment will be much smaller).

[0100] Further, the critical stability margin evaluation includes the following steps:

[0101] 1. Device parameter update: update the electro-thermal parameters of chip 1 according to the updated xi_edge(1). Exemplarily, assuming the mapping coefficients c_r = 0.01 mΩ, c_h = 0.001 K / W; R_on(1) = 2.5 + 0.01 x 0.1553 ≈ 2.50155 mΩ; R_th(1) = 0.1 + 0.001 x 0.1553 ≈ 0.100155 K / W.

[0102] 2. Closed-loop gain and margin calculation: at the current new operating point (95°C steady state) and updated device parameters, calculate the partial derivatives of each item of the K_loop matrix by the small perturbation method of embodiment four. Further, through calculation, the closed-loop gain matrix (illustrative value) is obtained: K_loop = [[0.02, -0.01], [-0.015, 0.025]]; further, the eigenvalues of the matrix are calculated, and the spectral radius rho(K_loop) ≈ 0.041 is obtained; further, the critical stability margin margin = 1 - rho(K_loop) = 1 - 0.041 = 0.959.

[0103] Through the technical scheme of this embodiment, the following technical effects are achieved: chip 1 has power over-stress and thermal resistance degradation, and its accumulated damage increases by 1.28e-5 (or 0.00128%); at the steady state point after the event, the critical stability margin of the entire power module is 0.959, and the system is in a very stable state; each subsequent thermal cycle will repeat steps three and four, continuously accumulating damage D_chip and updating state xi_edge, and continuously monitoring the trend of margin change, thereby achieving dynamic, online and high-precision prediction of the service life of the photovoltaic inverter.

[0104] In some optional embodiments, in addition to using the gradient descent-based alternating minimization algorithm, other advanced optimization techniques can also be used for the construction and solution of the inverse problem in step S202.

[0105] As an optional solution, optimization based on a proxy model: a machine learning proxy model (Surrogate Model), such as Gaussian Process Regression or neural network, is trained in advance through a large number of simulations (e.g., tens of thousands of simulations of different M and delta_theta parameter combinations on the callable thermal network model_base) or a small number of experiments. The proxy model can quickly and approximately predict the transient electrical fingerprint and temperature response under given M and delta_theta. Subsequently, the inverse problem can be converted into global optimization (e.g., using particle swarm optimization or genetic algorithm) of this computationally inexpensive proxy model, so that the global optimal solution can be found more efficiently, especially suitable for scenarios where the cost function form is complex or there are multiple local extreme points.

[0106] As an optional solution, Bayesian inference method: the inverse problem is solved in the Bayesian framework. In this framework, M and delta_theta are regarded as random variables to be inferred, and the initial value obtained by linear mapping in Embodiment II can be used as the prior distribution thereof. Through sampling methods such as Markov Chain Monte Carlo (MCMC), the posterior probability distribution of M and delta_theta can be obtained. This solution not only gives a point estimate of the optimal value, but also quantifies the uncertainty of the solution (i.e., the confidence interval), so that the reliability of the evaluation result can be measured, providing more information for subsequent risk decision-making.

[0107] In addition, the transient electrical fingerprint set in step S201 can also be further expanded. Optionally, gate charge characteristics (such as changes in Miller plateau duration), short window energy integral changes of switching loss, or spectral features (such as shifts in the main resonance peak frequency) obtained by fast Fourier transform (FFT) analysis of the switching ringing waveform can be introduced as additional fingerprints, in order to capture the degradation information of the device from more dimensions, thereby improving the constraint strength and solving accuracy of the inverse problem.

[0108] In some optional embodiments, the damage increment function g1 and the damage mapping function g2 described in steps S302 and S304 can have different function forms.

[0109] As an alternative, instead of the damage increment function, the increment function g1 can also take a non-linear form based on physical mechanism, in addition to the linear weighted form. For example, the local stress term sigma_local can be coupled with the temperature swing Delta_T through a product term or an exponential term, to better simulate the accelerated effect of stress on fatigue life at high temperature. For example, the relevant term in Delta_xi_edge(k) can be expressed as p_coup x (sigma_local(k)) n × (Delta_T(k)) m where p_coup, n, m are calibration coefficients.

[0110] As an alternative, instead of the damage mapping function, the mapping function g2 can also take other functions that can describe S-shaped cumulative curves, such as Logistic Function or Gompertz Function, in addition to the polynomial form, which can be used to simulate the complete process of damage accumulation, which is slow at the beginning, accelerated in the middle, and again tends to saturation due to material hardening and other factors, and may be more suitable for describing the entire fatigue life cycle of some materials.

[0111] As a preferred embodiment, the damage accumulation model of the present application can also be extended to include other degradation modes in addition to thermal fatigue. Optionally, a parallel internal state variable xi_cosmic_ray is introduced to accumulate the probabilistic damage caused by single event burnout (SEB) or single event gate rupture (SEGR) due to cosmic ray or other high-energy particle bombardment. The final total damage will be the comprehensive result of multiple state variables such as xi_edge and xi_cosmic_ray, thus forming a multi-physical field, multi-degradation mode coupled life prediction model.

[0112] In some alternative embodiments, in addition to calculating the spectral radius of the closed-loop gain matrix, the method for evaluating the stability of the system in step S402 can also use other stability criteria in control theory.

[0113] As an alternative, stability analysis based on system state space model: the electro-thermal coupled system of parallel chips is described as a linear time-varying (LTV) or nonlinear state space model; further, the stability of the system can be judged by calculating the real part of the eigenvalues (for continuous systems) or the modulus (for discrete systems) of the system matrix after linearization of the system at the current operating point, and all eigenvalues need to satisfy the stability condition. Optionally, Lyapunov stability theory is used to directly determine the stability of the system by finding a Lyapunov function, which is particularly suitable for nonlinear systems.

[0114] Among the above-mentioned alternatives and variations, event-driven physical closed-loop inversion, path-dependent damage accumulation, and critical criteria based on system stability are universal and scalable. Their specific implementation methods can be flexibly adjusted according to the accuracy requirements of the application scenario, the limitations of computing resources, and the available prior knowledge.

[0115] This embodiment details how to determine the power redistribution matrix M(t) and the thermal network parameter perturbation deltaTheta(t) by constructing and solving an inverse problem with transient electrical fingerprints as the core input, including the following steps:

[0116] Step S241, cost function assembly and dimensional normalization

[0117] Construct and solve an inverse problem with transient electrical fingerprint as input, the solution of which serves as a physical state parameter and includes at least: a power redistribution matrix to characterize the uneven heating power among parallel chips; or thermal network parameter perturbation to characterize the degradation of the thermal path of the device.

[0118] In this embodiment, it is necessary to normalize the observations of different dimensions. Specifically, the high-resolution temperature trajectory Tj_highres obtained in Embodiment 1 and the transient fingerprint sequence extracted in step S201 of Embodiment 2 are read and normalized. Preferably, the temperature normalization adopts the formula: Tnorm = (Tj_highres - Tref) / Tscale, where Tref is the reference temperature, usually taken as the ambient temperature of 25℃; Tscale is the temperature scale factor, taken as a typical temperature swing of 50K. Preferably, the transient fingerprint normalization adopts the formula: Enorm = (E - Emean) / Estd, where E is the original fingerprint vector, containing four components: delta_Vce_on, V_ov, zeta, and delta_Zth; Emean is the mean of each component in the healthy state; and Estd is the standard deviation.

[0119] Furthermore, the joint cost function J is then assembled: J = wT × sum_t[(Tnorm(t) - Tsim_norm(M, deltaTheta, t)) 2 ] + wE × sum_i[(Enorm_i - Esim_norm_i(M,deltaTheta)) 2+ lambda1 x TV(M) + lambda2 x TV(deltaTheta); where wT is the temperature term weight, taking value 0.6; wE is the electrical transient term weight, taking value 0.4; Tsim_norm is the normalized temperature obtained by simulating the thermal network model under given M and deltaTheta; Esim_norm is the transient fingerprint predicted according to M and deltaTheta; TV(M) is the total variation of the power redistribution matrix, calculated as TV(M) = sum_t|M(t+1) - M(t)|; lambda1 takes 1e-3, and lambda2 takes 5e-4, for constraining the smoothness of the solution.

[0120] Through the technical scheme of the embodiment, the unified optimization problem of multi-source heterogeneous data is solved. Through normalization processing, the temperature error and the electrical error can be compared and weighted under the same scale, avoiding that a certain item dominates the optimization process. The introduction of the total variation regular term prevents the parameters from jumping sharply in time, which conforms to the physical law of gradual change of device degradation.

[0121] In some optional embodiments, the weights wT and wE can be dynamically adjusted according to the data quality. Exemplarily, when the transient event is weak, wE can be increased to 0.6, and wT can be correspondingly reduced to 0.4, to rely more on electrical information. Further, an adaptive weight mechanism can be introduced: wT(k) = sigma_E 2 / (sigma_T 2 + sigma_E 2 ), where sigma_T and sigma_E are the variances of the temperature and electrical prediction errors, respectively.

[0122] Step S242, constraint definition and initial value loading

[0123] In the embodiment, through a linear prior mapping, the power redistribution matrix initial value and the thermal network parameter perturbation initial value are directly generated from the transient electrical fingerprint; further, the power redistribution matrix initial value and the thermal network parameter perturbation initial value are taken as the starting point of the iterative solution of the inverse problem.

[0124] Further, the transient fingerprint sequence is read, and a pre-calibrated linear mapping is used to generate the initial values. For the power redistribution matrix, the initial value of the i-th chip is calculated as: M_i(0) = a0_i + a1_i×delta_Vce_on + a2_i×V_ov + a3_i×zeta + a4_i×delta_Zth; the mapping coefficients are determined through offline experiments: a0_i = 0; a1_i = 0.008 / mV, indicating that the power sharing increases by 0.8% for every 1mV increase in the on-voltage drop; a2_i = -0.002 / V, the overshoot voltage is negatively correlated with the power; a3_i = 0.15, the faster the ring decay, the greater the damping; a4_i = 0.5 / (K / W), an increase in thermal resistance directly leads to power redistribution.

[0125] Further, the initial value of the thermal network parameter perturbation is calculated as: deltaTheta_Rj(0) = b0 + b1×delta_Zth + b2×zeta; deltaTheta_Cj(0) = c0 + c1×delta_Zth; where deltaTheta_Rj is the change in junction shell thermal resistance, and deltaTheta_Cj is the change in thermal capacity; b0 = 0, b1 = 0.8, b2 = 0.2; c0 = 0, c1 = -0.1, the thermal capacity usually decreases slightly with aging.

[0126] Further, physical constraints are applied to ensure the reasonableness of the solution. As a preferred embodiment, the power conservation constraint sum_i[M_i(t)] = 0 is realized by M_i = M_i - mean(M); the boundary constraint: -0.2 ≤ M_i ≤ 0.2, limits the power deviation of a single chip to no more than 20%; the non-negative constraint: deltaTheta_Rj ≥ 0, deltaTheta_Cj ≥ -0.1*Cj_nom, the thermal resistance only increases and the thermal capacity decreases within a limited range.

[0127] The technical solution of the present embodiment provides a starting point close to the true solution through prior mapping, greatly reducing the number of iterations. The setting of the constraint conditions is based on the physical characteristics of the device, preventing unreasonable solutions from being obtained in the optimization process. Experiments show that the algorithm with prior initial values converges 3-5 times faster than random initialization.

[0128] Alternatively, when the historical data is sufficient, a machine learning method can be used to improve the linear mapping. For example, a shallow neural network is trained, with the transient fingerprint as input and the initial values of M and deltaTheta as output. The network structure can be selected as [4, 16, 8, N+2], where N is the number of chips. Further, a Bayesian framework can be introduced, regarding the mapping coefficients a, b, and c as random variables, and updating their posterior distribution through historical data.

[0129] Step S243, Jacobian and gradient calculation

[0130] In this embodiment, the inverse problem is solved by minimizing a joint cost function, which includes at least a temperature difference term representing the difference between the simulated temperature and the high-resolution temperature trajectory, and a total variation regularization term for constraining the physical reasonability of the power redistribution matrix or the perturbation of the thermal network parameters.

[0131] In this step of the embodiment, the finite difference method is used to calculate the Jacobian matrix. A perturbation epsilon_M = 0.001 is applied to the i-th element of M, and the change in the temperature response is calculated: dTsim_dM_i = [Tsim(M + epsilon_M×e_i, deltaTheta) - Tsim(M, deltaTheta)] / epsilon_M; where e_i is the i-th unit vector.

[0132] Further, the Jacobian of deltaTheta is calculated: dTsim_dTheta_j = [Tsim(M, deltaTheta+ epsilon_Theta×e_j)-Tsim(M,deltaTheta)] / epsilon_Theta; epsilon_Theta takes 0.0001×Theta_nom, and Theta_nom is the nominal value of the thermal parameter.

[0133] Further, the Jacobian of the electrical transient to the parameters is calculated through an analytical relationship. For example, the sensitivity of the on-state voltage drop to the power distribution: dVce_dM_i = Ron_nom × I_nom × (1 + alpha_T × deltaT_i); where Ron_nom is the nominal on-state resistance 2.5 mΩ; I_nom is the rated current; alpha_T is the temperature coefficient 0.003 / K; and deltaT_i is the temperature deviation of the i-th chip.

[0134] Further, the total gradient is calculated: grad_M = 2 x wT x J_T' x (Tsim - Tj) + 2 x wE x J_E' x (Esim - E) + lambda1 x sign(diff(M)); grad_Theta = 2 x wT x J_Theta' x (Tsim - Tj) + 2 x wE x J_E_Theta' x (Esim - E) + lambda2 x sign(diff(deltaTheta)); where J_T, J_E, J_Theta, J_E_Theta are the corresponding Jacobian matrices; diff() is the difference operator; sign() is the sign function, used to calculate the sub-gradient of the TV term.

[0135] According to further improvements of the present application, to improve numerical stability, the gradient is clipped: ||grad_M||_inf <= 0.05, ||grad_Theta||_inf <= 0.01; to prevent oscillation caused by excessively large single-step update.

[0136] Although the finite difference method is general, it has a large amount of calculation, and needs (N+M) times of forward simulation each time, where N is the dimension of M and M is the dimension of deltaTheta. Gradient clipping ensures the robustness of the algorithm, which can stably converge even in ill-conditioned conditions.

[0137] In some optional embodiments, automatic differentiation techniques can be used to accelerate Jacobian calculation. By constructing a computation graph and backpropagation, the gradient can be accurately and efficiently obtained. Further, when the model has a special structure, the gradient calculation complexity can be reduced to O(1) times of forward solution by using the adjoint method.

[0138] Step S244, alternating minimization and projection update

[0139] In the step of the present embodiment, an alternating minimization strategy is used for iterative solution. Fixing deltaTheta, update M: M_new = M_old - eta_M x grad_M; M_proj = Proj_M(M_new); the learning rate eta_M is initially set to 0.02, and an adaptive adjustment strategy is used. The implementation of the projection operator Proj_M includes two steps: perform the conservation projection, M_i = M_i - mean(M); further, perform the boundary projection, M_i = max(-0.2, min(0.2, M_i)).

[0140] Further, fixing the updated M, update deltaTheta: ; deltaTheta_proj = Proj_Theta(deltaTheta_new); eta_Theta is initially set to 0.005. The projection operator Proj_Theta performs non-negative projection: deltaTheta_Rj = max(0, deltaTheta_Rj); deltaTheta_Cj = max(-0.1 x Cj_nom, deltaTheta_Cj).

[0141] Further, the improvement of the cost function is evaluated after each iteration. If J_new > J_old, backtracking line search is triggered: eta = 0.5 x eta, the step update is recalculated until J_new < J_old or eta < eta_min. eta_min is set to le-5.

[0142] Alternating minimization decomposes the complex joint optimization problem into two relatively simple sub-problems, each of which is convex or approximately convex, easy to solve. The projection operation ensures that each iteration is within the feasible region, avoiding additional penalty terms. Experiments show that the algorithm usually converges within 20-50 iterations.

[0143] Optionally, momentum can be used to accelerate convergence. A momentum term is introduced: M_new = M_old - eta_Mgrad_M + beta_M(M_old - M_prev), where beta_M takes 0.8. Further, a quasi-Newton method like L-BFGS can be used to approximate the inverse of the Hessian matrix using historical gradient information, achieving super-linear convergence.

[0144] Step S245, convergence criterion and termination condition

[0145] In one specific embodiment, the historical sequence of the cost function is read, and the relative improvement degree is calculated: rel_improve = (J_prev - J_current) / max(J_prev, le-6).

[0146] Further, the iteration is terminated when any of the following conditions is met: the relative improvement degree rel_improve < le-4, indicating that the optimal solution has been approached; the number of iterations iter ≥ 100, preventing infinite loops; the gradient norm ||grad_M||2 + ||grad_Theta||2 < le-5, indicating that the first-order optimal condition is reached.

[0147] In some optional embodiments, to improve the robustness of the convergence criterion, a moving average is used: J_smooth(k) = 0.7J_smooth(k-1) + 0.3J(k). The relative improvement is calculated using the smoothed cost function to avoid premature termination caused by noise.

[0148] Further, the convergence information is recorded for diagnosis: the final cost function value J_final; the total number of iterations iter_total; the error components wTerr_T, wEerr_E, lambda1TV_M, lambda2TV_Theta. If err_T > 0.1 or err_E > 0.15, it is marked as weak convergence, and the data quality or weight adjustment needs to be checked.

[0149] The design of the convergence criterion balances accuracy and efficiency. Too strict convergence conditions will cause unnecessary calculations, and too loose will affect the quality of the results. Moving average eliminates small fluctuations in the cost function, making the convergence criterion more stable. Sub-error records help problem diagnosis and parameter tuning.

[0150] In some optional embodiments, an adaptive convergence threshold can be used. According to the noise level estimate: tol_adaptive = k_tol x sqrt(var_noise), where var_noise is the observation noise variance, and k_tol takes 3-5. Further, an early stopping mechanism can be introduced to terminate when the validation set error starts to rise, preventing overfitting.

[0151] Step S246, observability evaluation

[0152] In this embodiment, after solving the inverse problem, the sensitivity of the problem solution is evaluated to generate an observability index; when the observability index is lower than a preset threshold, an adaptive perturbation plan is generated, including: calculating the observability index J_obs: J_obs = gamma1 x ||dTsim_dM||_F + gamma2 x ||dTsim_dTheta||_F + gamma3 x cond(H) -1 ; where ||·||_F is the Frobenius norm; gamma1 = gamma2 = 0.4, gamma3 = 0.2; H is the approximate Hessian matrix, H_ij = sum_t[(dTsim_t / dxi) x (dTsim_t / dxj)], x is [M; deltaTheta]; cond(H) is the condition number, reflecting the degree of ill-conditioning of the problem.

[0153] Further, when J_obs < J_obs_min (J_obs_min takes 0.05), it is determined to be weak observability. At this time, a perturbation plan is generated: u_pert = beta0 + beta1 x v_max + beta2 x randn() x v_std; wherein, u_pert is the recommended perturbation signal; beta0 = 0.005 x P_nom, the basic perturbation amount; v_max is the perturbation direction that maximizes J_obs, which is obtained by solving max_v(||dTsim / dv||); beta1 = 0.01; v_std is a small random perturbation, beta2 = 0.002; randn() is a standard normal random number.

[0154] Further, the perturbation constraint ensures system safety: ||u_pert||_2 <= 0.01 x P_nom, the perturbation power does not exceed 1% of the rated value; the perturbation duration t_pert <= 200 ms; the perturbation frequency f_pert <= 1 times / hour.

[0155] Observability evaluation quantifies the reliability of inverse problem solution. Low observability means that the current data cannot effectively distinguish different parameter combinations, and active excitation is needed. The design of the perturbation plan takes into account the information gain and system safety, and improves the subsequent identification accuracy through short and small perturbation.

[0156] Alternatively, the Fisher information matrix can be used to evaluate the observability: FIM = E[grad x grad'], wherein E[] is the expectation operator. The minimum eigenvalue of FIM reflects the amount of information in the most difficult to identify direction. Further, an optimal experiment can be designed to determine the best excitation signal by solving max(det(FIM)) or max(lambda_min(FIM)).

[0157] The embodiment details how to perform path-dependent recursive update on the internal state variable Sedge representing the cumulative effect of the chip edge singular zone, including the following steps:

[0158] Step S331, thermal cycle segmentation and feature extraction

[0159] The modified temperature history is divided into a series of discrete thermal cycles; and for each thermal cycle, a recursive update is performed on an internal state variable representing the cumulative effect of the chip edge singular zone.

[0160] In this embodiment, the corrected temperature history Tj_corr obtained in Example Five is read and cycle identification is performed using rainflow counting method. Specifically, the temperature time series is treated as a stress-time curve, and cycles are extracted by the following rules: local extrema are detected to form a sequence of peaks and valleys; starting from the largest peak, paired valleys are searched on both sides to form primary cycles; the same rule is recursively applied to the remaining peak-valley pairs to extract secondary cycles; a minimum temperature difference threshold Delta_T_min = 5K is set to filter out pseudo cycles with too small amplitudes.

[0161] Further, for each identified thermal cycle k, statistical features are extracted: ; t_dwell_high(k) = dwell time in high temperature zone (> 0.8 x T_max), t_dwell_low(k) = dwell time in low temperature zone (< 0.2 x T_max + 0.8 x T_min); where t_rise and t_fall are the rise and fall times, respectively. These features comprehensively describe the intensity, rate and time characteristics of the thermal cycle.

[0162] The advantage of rainflow counting method is that it can correctly handle nested cycles and irregular load sequences, which are common in actual photovoltaic operation (such as intermittent cloud cover). The six extracted feature parameters cover the main factors affecting fatigue damage, providing complete input for subsequent state updating.

[0163] In some optional embodiments, wavelet decomposition can be used to preprocess the temperature signal to separate cycles of different time scales. For example, db4 wavelet is used for 5-level decomposition, and the 3rd-5th level reconstructed signals are retained to filter out high-frequency noise and ultra-low-frequency drift. Further, for scenarios with dense fast thermal cycles, equivalent damage method can be used to combine small cycles to reduce computational load.

[0164] Step S332, local stress calculation

[0165] The recursive update of the internal state variable is performed based on a path-dependent incremental function, and the inputs of the incremental function include: the rise and fall slopes and dwell times extracted from the thermal cycles; and also includes a local stress term that depends on the physical state parameters.

[0166] In this embodiment, M(t) and deltaTheta(t) output by Example Seven are read, and the local stress Olocal(k) of the kth cycle is calculated. The stress is composed of multiple components:

[0167] ​​​Thermal-mechanical stress component: sigma_thermal(k) = E_chip x alpha_CTE x Delta_T(k) x (1 + M_i(k); where E_chip is the chip elastic modulus, taken as 130 GPa; alpha_CTE is the thermal expansion coefficient mismatch, the difference between silicon and copper substrate is about 13 ppm / K; M_i(k) is the average power deviation during this cycle.

[0168] Interface stress concentration component: sigma_interface(k) = K_stress x sqrt(a_crack) x sigma_thermal(k) x (1 + 0.5 x deltaTheta_Rj(k) / Rj_nom); where K_stress is the stress intensity factor, a typical value is 1.2; a_crack is the equivalent crack length, related to Sedge; deltaTheta_Rj(k) / Rj_nom reflects the influence of interface degradation on stress concentration.

[0169] Creep stress component: sigma_creep(k) = sigma_0 x exp(-Q / (R_gas x T_mean(k))) x (t_dwell_high(k)) n_creep ; where sigma_0 is the reference stress, taken as 10 MPa; Q is the activation energy, taken as 40 kJ / mol for solder; R_gas is the gas constant 8.314 J / (mol.K); n_creep is the creep index, taken as 0.3.

[0170] Total local stress: Olocal(k) = w1 x sigma_thermal(k) + w2 x sigma_interface(k) + w3 x sigma_creep(k); weight coefficients w1 = 0.5, w2 = 0.3, w3 = 0.2, determined according to failure mode analysis.

[0171] The calculation of local stress converts power redistribution and thermal parameter changes into mechanical quantities, establishes the coupling relationship of electricity-heat-force, so that the update of Sedge not only depends on the temperature cycle itself, but also responds to the real degradation state inside the device, realizing physical closed loop.

[0172] Optionally, a plastic stress component can be introduced, and when sigma_thermal exceeds the yield strength, a Ramberg-Osgood model is used to calculate the plastic strain. Further, for the bonding wire, its bending stress can be calculated separately, considering the influence of wire diameter, arc height and span.

[0173] Step S333, directional gain and path dependence update

[0174] The path-dependent increment function differentiates the heating and cooling stages of the thermal cycle, specifically including: in the heating stage, the increment function adopts a first gain coefficient to weight the contribution of the heating and cooling slope to the damage increment; and in the cooling stage, a second gain coefficient different from the first gain coefficient is adopted.

[0175] Further, for the kth cycle, the internal state increment is calculated: Delta_Sedge(k) = f_base(k) + f_rate(k) + f_dwell(k) + f_stress(k) - f_relax(k); the base damage term: f_base(k) = p0 + p1 x T_mean(k) / T_ref + p2 x (Delta_T(k) / Delta_T_ref) m1 ; wherein, p0 = 1e-5, base damage rate; p1 = 2e-3, uniform temperature influence coefficient; p2 = 5e-3, temperature swing coefficient; T_ref = 100℃, Delta_T_ref = 50K, reference value; m1 = 1.5, nonlinear index. Rate-dependent term (reflecting directionality): f_rate(k) = p3 x [r_up x max(dT_dt_up(k), 0) m2 + r_down x |min(dT_dt_down(k), 0)| m2 ]; wherein, p3 = 1e-4; r_up = 1.3, heating gain; r_down = 0.7, cooling gain; m2 = 1.2. r_up > r_down reflects that the ratchet effect caused by compressive stress during heating is stronger than the tensile stress recovery during cooling. Dwell damage term: f_dwell(k) = p4 x [t_dwell_high(k) x exp(-t_dwell_high(k) / tau_h) + 0.1 x t_dwell_low(k)]; wherein, p4 = 5e-5; tau_h = 3600s, high-temperature dwell characteristic time. The exponential term indicates that the damage is fast at the beginning of dwell and tends to be saturated later. Stress coupling term: f_stress(k) = p5 x Olocal(k) / sigma_ref; wherein, p5 = 1e-3; sigma_ref = 100MPa, reference stress. Relaxation term (reflecting memory decay): f_relax(k) = beta_relax x Sedge(k) x exp(-T_mean(k) / T_relax); wherein, beta_relax = 0.01, relaxation coefficient; T_relax = 50℃. Relaxation is more obvious at low temperature and almost no relaxation at high temperature.

[0176] Furthermore, update the internal state: Sedge(k+1) = Sedge(k) + Delta_Sedge(k); Sedge(k+1) = min(max(Sedge(k+1), 0), 1), restricted to the interval [0,1];

[0177] The directional gain r_up ≠ r_down is key to simulating the thermal ratcheting effect. Experiments show that, under the same temperature swing, rapid heating followed by slow cooling causes 20-30% more damage than the reverse process. The relaxation term reflects the material's self-healing properties, especially the partial release of stress during long-term residence at low temperatures. This path dependence allows the model to distinguish the cumulative effects of different loading histories.

[0178] In some alternative implementations, a cyclic sequence effect can be introduced. The memory kernel function is defined as K(k,j) = exp(-(kj) / N_mem), where N_mem is the memory length (e.g., 10 cycles). State updates consider historical effects: Delta_Sedge(k) = sum_j[K(k,j)×f(j)]. Furthermore, the memory effect can be described using fractional derivatives, with the 0.7 derivative calculated using Caputo's definition.

[0179] Step S334, lift-off event weighted adjustment

[0180] From the transient electrical fingerprint, a list of lift-off events is identified and generated; and when performing recursive updates to the intrinsic state variables, if it is determined that the current thermal cycle is within the event neighborhood of the lift-off event list, the damage increment calculated by the path-dependent incremental function is adjusted by amplifying the weights.

[0181] In this embodiment, the lift-off event list generated in Embodiment 2 is read, which records the event time t_lo and the initial weight w_lo0. It is then determined whether the k-th cycle is within the event neighborhood: in_neighborhood(k) = (|t_cycle(k) - t_lo| < t_window); where t_cycle(k) is the center time of the k-th cycle; t_window is the neighborhood window, set to 3600s.

[0182] Further, the time-varying weights are calculated as follows: w_lo(k) = alpha_lo + K_lo × w_lo0 × exp(-|t_cycle(k) - t_lo| / tau_lo); where alpha_lo = 1.0, the base weight; K_lo = 0.5, the amplification factor; and tau_lo = 1800s, which affects the decay time.

[0183] Further, the state increment is adjusted by a weight: Delta_Sedge_adj(k) = Delta_Sedge(k) x w_lo(k) if the k-th cycle is within the event neighborhood; otherwise, Delta_Sedge_adj(k) = Delta_Sedge(k).

[0184] The lift-off event is detected based on a sudden change in the transient fingerprint, when delta_Vce_on > 10mV and lasts for more than 3 switching cycles. The impact of the event is propagated through an exponential decay function, reflecting the gradual stabilization of the current redistribution after the wire pull-off.

[0185] The identification and weighting of lift-off events is the key to connect transient diagnostics with long-term degradation, as the wire pull-off not only immediately changes the current distribution, but also accelerates the fatigue in the neighboring regions in subsequent cycles. The time-varying weight function quantifies the stress effect after trauma.

[0186] Optionally, the lift-off severity can be graded for different handling. Mild (delta_Vce_on < 15mV): K_lo = 0.3; Moderate (15-25mV): K_lo = 0.5; Severe (>25mV): K_lo = 0.8. Further, the superposition effect of multiple lift-off events can be considered, with a product form w_lo_total(k) = prod_i[w_lo_i(k)].

[0187] Step S335, damage mapping and thresholding

[0188] The state increment calculated by the increment function is mapped to the incremental damage for this thermal cycle, including: applying a pre-set damage threshold to the state increment, only when the state increment exceeds the threshold, it is confirmed to contribute to the damage; and the part exceeding the damage threshold is processed by a non-linear mapping function, including:

[0189] In this embodiment, the damage threshold is calculated: Delta_Sedge_eff(k) = max(Delta_Sedge_adj(k) - S_th, 0); where S_th is the fatigue threshold, taking 2e-5. Small cycles below this value do not produce permanent damage.

[0190] Further, the non-linear damage mapping is calculated: Delta_D(k) = q0 + q1 x Delta_Sedge_eff(k) + q2 x (Delta_Sedge_eff(k)) gamma; where q0=0 (no damage for no state change); q1=0.1, linear coefficient; q2=10, non-linear coefficient; gamma=2.0, acceleration exponent.

[0191] Further, a damage rate limit is introduced to prevent numerical problems: Delta_D(k) = min(Delta_D(k), D_max_per_cycle); D_max_per_cycle takes 0.01, i.e. no more than 1% damage per cycle.

[0192] Further, cumulative damage update: D_total(i) = D_total(i) + Delta_D(k); L_consumed(i) = D_total(i) / D_critical; where i is the chip index; D_critical=1.0 is the failure criterion; L_consumed is the consumed life proportion.

[0193] The physical meaning of threshold S_th corresponds to the fatigue limit of the material, i.e. there is a stress level below which the material can withstand infinite cycles. The non-linear mapping reflects the acceleration property of damage accumulation: as Sedge increases, the material enters the crack propagation region, and the damage rate rises sharply. Experimental data show that when Sedge > 0.6, the damage rate can reach more than 10 times that of the initial stage.

[0194] In some optional embodiments, a mapping in the form of S-N curve can be used: Delta_D(k) = (Delta_Sedge_eff(k) / S_f) -1 / b where S_f is the fatigue strength coefficient, and b is the fatigue strength exponent. Further, probabilistic damage can be introduced, and Weibull distribution can be used to describe the randomness of damage.

[0195] Step S337, parameter calibration process

[0196] This step details the calibration method of each parameter in the g1 and g2 functions. The calibration is carried out in two stages.

[0197] Exemplarily, stage one is threshold and directionality parameter calibration, including: building a thermal cycle test bench, applying different characteristic standard cycles to the same batch of devices: group 1: Delta_T=50K, dT / dt=±10K / s, symmetric cycle, used for calibrating p2 reference value; group 2: Delta_T=50K, 20K / s for heating, 5K / s for cooling, used for calibrating r_up; group 3: Delta_T=50K, 5K / s for heating, 20K / s for cooling, used for calibrating r_down; group 4: Delta_T=30K, high temperature residence for 10 min, used for calibrating p4 and tau_h; group 5: Delta_T=10K, continuous fast cycle, used for calibrating S_th.

[0198] In the embodiment, each group is cycled 1000-5000 times, and the first cracking time is detected by acoustic emission. Define the objective function: J_calib1=sum_g[(N_crack_pred(g)-N_crack_obs(g)) 2 / N_crack_obs(g) 2 ]; Wherein, g is the group index; N_crack_pred and N_crack_obs are the predicted and observed crack cycle numbers respectively. Grid search is used to optimize the parameters, search range: r_up∈[1.0,2.0], r_down∈[0.5,1.0], S_th∈[1e-5,1e-4].

[0199] In another example, stage two is damage evolution parameter calibration, including: continuing to cycle the sample with known initial cracks to failure, and monitoring the change of on-resistance. In the embodiment, Ron is measured every 500 cycles, and the Sedge-Ron relationship is established. Define the objective function:

[0200] J_calib2=sum_t[(Ron_pred(t)-Ron_meas(t)) 2 ]+lambda_smooth×TV(p_vector); Wherein, t is the measurement time; lambda_smooth=0.01, smoothing regular term; p_vector=[p0,p1,p3,p5,q1,q2,gamma]. Levenberg-Marquardt algorithm is used for optimization.

[0201] As a preferred embodiment, the verification of the calibration result uses leave-one-out cross-validation. The calibration samples are divided into N groups, and N-1 groups are used for calibration, and the remaining 1 group is used for verification, which is repeated N times. The average verification error is required to be <15%, and the standard deviation is <5%.

[0202] Exemplarily, typical calibration results are as follows: p0=1e-5±2e-6; p1=(2.0±0.3)e-3; p2=(5.0±0.5)e-3; p3=(1.0±0.2)e-4; p4=(5.0±1.0)e-5; p5=(1.0±0.2)e-3; r_up=1.3±0.1; r_down=0.7±0.1; S_th=(2.0±0.5)e-5; q1=0.10±0.02; q2=10±2; gamma=2.0±0.2; beta_relax=0.010±0.002.

[0203] Parameter calibration is the key to ensure the accuracy of model prediction. The two-stage method separates the triggering mechanism (threshold parameter) and evolution law (damage parameter), reducing the complexity of calibration. The calibration data should cover the typical working conditions of actual operation, including extreme cases. The uncertainty of parameters is quantified by confidence interval, which can be used for interval estimation of life prediction.

[0204] Optionally, a Bayesian calibration method can be used, which regards parameters as random variables and obtains posterior distribution through MCMC sampling. Further, a correlation model between parameters and device batch and package type can be established to realize fast migration calibration.

[0205] The present scheme solves the problem that the traditional method cannot utilize millisecond transient event information, leading to early degradation detection. Specifically, the method uses the fingerprint of transient electrical events as an effective input, and through inversion calculation, establishes a deterministic correlation between external electrical response and internal physical state (such as power redistribution matrix and thermal network parameter perturbation) that cannot be directly observed. It realizes online and quantitative perception of early degradation precursors such as wire lift-off, and captures the key aging accelerator-chip power imbalance caused by it, filling the observation blind area of the traditional method.

[0206] Further, the present scheme solves the problem that the traditional damage model lacks path memory and cannot simulate thermal ratchet effect through the path-dependent thermal ratchet incremental kernel of embodiment three. By introducing an internal state variable that can remember the history load, its update function not only considers the asymmetry of temperature rise and fall, but also couples the local stress corresponding to the real physical state inverted in the previous step, so that damage accumulation is no longer a simple superposition of isolated thermal cycles, but a more physically realistic, directional and history-dependent process, improving the accuracy of life prediction under complex and variable working conditions.

[0207] Further, to solve the problem of inaccurate system level instability early warning, the application solves it by using the electric-thermal coupling critical gain criterion. The criterion maps the microscopic cumulative damage to the change of macroscopic electric-thermal parameters, and calculates the closed-loop gain spectrum radius which describes and connects the positive feedback strength of the system, thereby quantifying the stability margin of the module tending to thermal runaway. A critical point prediction method with clear physical meaning is provided from device degradation to system bifurcation, which improves the reliability and advance of the early warning by promoting the judgment of life termination from the experience threshold to the level of system stability analysis.

[0208] The preferred embodiments of the application are described in detail above, but the application is not limited to the specific details of the above-described embodiments. Within the technical concept of the application, various equivalent transformations of the technical solutions of the application can be made, and these equivalent transformations all belong to the protection scope of the application.

Claims

1. A photovoltaic inverter lifetime prediction method, characterized by, The method comprises: acquiring operating waveform data of a photovoltaic inverter and device structure parameters, the operating waveform data containing transient electrical events; based on the transient electrical events and in combination with the device structure parameters, inversely determining physical state parameters representing internal characteristics of the power device, and generating a corrected temperature history; updating a damage state representing cumulative effects and having path dependence according to the corrected temperature history and the physical state parameters; wherein the updating of the damage state representing cumulative effects and having path dependence comprises: dividing the corrected temperature history into discrete thermal cycles; and for each thermal cycle, performing a recursive update on an internal state variable representing cumulative effects of the chip edge singular zone according to a path-dependent incremental function; the inputs of the incremental function include: the temperature rise and fall slope and the residence time extracted from the thermal cycle, and a local stress term dependent on the physical state parameters; based on the damage state having path dependence, evaluating the critical stability margin of the power module and determining the remaining life of the module; wherein the evaluation of the critical stability margin of the power module comprises: updating the time-varying equivalent electrical parameters and the time-varying equivalent thermal parameters of the power device according to the damage state having path dependence; wherein the time-varying equivalent electrical parameters at least include the on-resistance, and the time-varying equivalent thermal parameters at least include the thermal resistance; based on the updated time-varying equivalent electrical parameters and the time-varying equivalent thermal parameters, constructing a model describing the electrical-thermal coupling effect of the parallel chip system to determine the closed-loop gain matrix therefrom; calculating the spectral radius of the closed-loop gain matrix, which is used to determine the critical stability margin.

2. The method of claim 1, wherein, The step of inversely determining the physical state parameters representing the internal characteristics of the power device comprises: for any transient electrical event, extracting a set of transient electrical fingerprints from the corresponding operating waveform data for quantifying the electrical transient response during the event; constructing and solving an inverse problem with the transient electrical fingerprints as input, the solution of which is the physical state parameter and at least includes: a power redistribution matrix representing the imbalance of heating power among parallel chips, or a thermal network parameter perturbation representing the degradation of the device thermal path; wherein, before solving the inverse problem, the power redistribution matrix initial value and the thermal network parameter perturbation initial value are generated through a linear prior mapping and used as the starting point for iterative solution.

3. The method of claim 2, wherein, The solution of the inverse problem is achieved by minimizing a joint cost function, which at least includes: a temperature difference term representing the difference between the simulated temperature and the high-resolution temperature trajectory, and a total variation regularization term for constraining the physical reasonableness of the power redistribution matrix or the thermal network parameter perturbation; after solving the inverse problem, the sensitivity of the solution of the problem is evaluated to generate an observability index; when the observability index is lower than a preset threshold, an adaptive perturbation plan is generated for actively adjusting the operating conditions or data sampling strategy in the subsequent evaluation period.

4. The method of claim 1, wherein, The path-dependent incremental function differentiates the heating and cooling stages of the thermal cycle, including: in the heating stage, the incremental function adopts a first gain coefficient to weight the contribution of the temperature rise and fall slope to the damage increment; In the cooling phase, a second gain coefficient different from the first gain coefficient is used to weight the contribution of the temperature ramp-up and ramp-down slopes, forming an asymmetric damage accumulation.

5. The method of claim 1, wherein, The method further comprises: From the transient electrical fingerprint, a list of lift-off events is identified and generated; When performing the recursive update of the intrinsic state variable, if it is determined that the current thermal cycle is within the event neighborhood of the list of lift-off events, the damage increment calculated by the path-dependent increment function is adjusted by amplifying the weight to simulate the accelerated degradation effect caused by the event.

6. The method of claim 1, wherein, The method further comprises: calculating the damage increment determined by the path-dependent increment function, and calculating the sensitivity of at least one of the average junction temperature, the temperature swing, the temperature ramp-up and ramp-down slopes, or the dwell time of the thermal cycle; Based on the calculation result of the sensitivity, a singularity sensitivity report is generated; The damage increment determined by the path-dependent increment function is calculated, specifically, the state increment is calculated by the increment function, and is mapped to the incremental damage of the thermal cycle, wherein the mapping process includes: applying a preset damage threshold to the state increment, and only when the state increment exceeds the threshold, it is confirmed that it contributes to the damage; for the part exceeding the damage threshold, the non-linear mapping function is processed to determine the incremental damage.

7. The method of claim 1, wherein, The method further comprises: Monitoring the critical stability margin, when it narrows to a preset risk threshold, generating a sampling and window scheduling suggestion; Feedback sampling and window scheduling suggestions are used to adjust the data sampling density or window selection strategy for transient electrical events; After determining the critical stability margin, a hysteresis criterion with different entry threshold and exit threshold is used to process the critical stability margin to generate a critical imbalance point marker; wherein when the critical stability margin is lower than the entry threshold, the marker is generated; when the margin rises above the exit threshold after the marker has been generated, the marker is removed.

8. The method of claim 1, wherein, The step of determining the module remaining life comprises: Based on the accumulation of the incremental damage, the remaining life of each chip is determined; Combined with the weakest value in all chip remaining lives and the critical stability margin, the module remaining life is determined. The step of determining the module remaining life comprises: Based on the accumulation of the incremental damage, the remaining life of each chip is determined; Combined with the weakest value in all chip remaining lives and the critical stability margin, the module remaining life is determined.

Images