Spectral Steady-State Analysis of Inverters With Temperature-Dependent Losses Using Harmonic Balance

Accurate calculation of semiconductor losses and temperature is the foundation of any design methodology for a power electronic converter. Computation accuracy and speed play a vital role if a large set of parameters needs to be considered. Averaged loss models often neglect the temperature dependence of transistors, leading to fast, but inaccurate results. In contrast, iterative methods and simulation tools, which can include temperature dependence, take significantly longer to compute, but yield more precise results. This paper presents a best of both worlds approach, by using the harmonic balance method to obtain the steady-state solution for any inverter topology including temperature dependent conduction and switching losses. The proposed method solves for the discrete Fourier series of the device temperature, by expressing the temperature dependence and operating parameters in the frequency domain. The set of equations for each coefficient is solved by a single matrix inversion, resulting in very fast computation for steady-state temperature cycles. The steady-state operation of over one thousand possible inverter designs is calculated within less than one minute, matching iterative simulation in device temperature, conduction and switching losses, at a fraction of the computation time. In addition, the method shows good agreement with temperature measurements of a three-phase silicon-carbide inverter.


I. INTRODUCTION
Optimization of power electronic converters often aims at multiple goals, such as efficiency, power density, lifetime of components or cost [1], [2]. This procedure often starts with a large number of potential operating parameters and components, in order to meet the desired specifications. For each set of parameters, the losses and performance of the converter components are either calculated based on averaged states or iterated in time for more precise results [3] [4]. Usually, averaged models are much faster than iterative ones, but iterative solutions can include non-linear effects often neglected in averaged models [5]. Depending on the severity of the nonlinearity, iterative simulation might be the only commonly available option towards correct results.
In case of power electronic devices, such as silicon or silicon-carbide (SiC) MOSFETs, accurate calculation of the total device losses is severely impeded by the temperature dependence of both conduction and switching losses [6], [7]. The temperature dependent losses can either be obtained through device characterization [8], or analytic modelling [9]. Given the half-bridge inverter in Fig. 1, each transistor experiences a temperature swing because of varying load conditions during each ac cycle. Because the device losses depend on temperature, an algebraic loop is formed: As the device heats up, the losses increase, which in turn further increases the temperature [4]. The convergence of this loop depends on the output frequency of the inverter as well as the thermal interface of the device. Given proper thermal design, the device will reach a periodic steady-state within the allowed operating parameters, after a certain number of cycles. Once in steady-state, the device temperature alternates between a constant minimum and maximum. In this steady-state, the average device losses determine the operating efficiency of the inverter. Furthermore, the temperature swing experienced by the transistor is relevant for lifetime estimation, as thermal expansion of the device, and subsequent bond-wire lift-off or solder-cracks are dominant failure mechanisms in power converters [10]. With non-linearity affecting efficiency and lifetime, precise models and tools are crucial for correct optimization.
One state-of-the-art tool for simulating non-linear thermal behaviour is PLECS [11] [12]. Starting from a given initial condition, it solves the differential equations for all systemstates step by step. To do this, PLECS describes the circuit as a set of state-space matrices denoting the change in electric quantities of all passive components, for all possible switching states of semiconductors. To reduce the matrix size and accelerate computation, PLECS can decouple parts of circuits where possible. More detailed information is given in [13]. While PLECS is faster than traditional time-based simulation, a search through a vast solution space still takes a significant amount of time, as several thermal cycles need to be simulated per parameter set until a steady state can be determined. This paper presents a closed-form solution for the periodic junction-temperature and losses of devices in a power electronic converter including temperature dependence. The method is based on harmonic balance, which thus far, has been applied in power electronics for calculation of distortion caused by non-linear loads [14]. However, the methodology can be adapted to include the non-linearity of losses in calculation of the device temperature. The solution is obtained through a single matrix inversion, making the proposed method orders of magnitude faster than any time-based iterative calculation. To include second-order temperature dependencies, the method can be applied twice to further increase the accuracy. The derivation of the proposed method is covered step-by-step in Section II. Afterwards, a potential use case is demonstrated in Section III. Here, the operating point of a three-phase grid-inverter is calculated for over 1000 possible parameter combinations. The results of the proposed method are compared to simulation done in PLECS, both in accuracy as well as computation speed. In Section IV, the proposed method is compared to measurements on a three-phase SiC inverter. Finally, Section V concludes the paper.

II. TEMPERATURE CALCULATION USING HARMONIC BALANCE
Harmonic balance has been used in circuit simulation to deal with non-linear components in otherwise linear networks [15]. The underlying assumption for harmonic balance is that a given input frequency only produces responses that can be expressed as integer multiples of itself. This allows to express all waveforms in the circuit as discrete Fourier series, and solve the response of each coefficient to the given input frequency. Because harmonic balance works in the frequency domain, differential equations from the time domain transform into algebraic operations. The resulting set of equations can be solved using standard solving methods.

A. POWER TRANSISTOR IN THE FREQUENCY DOMAIN
The basis for the presented method is the equation for the junction temperature of a device T j (t ), given varying conduction losses P c (t ) and switching energies E on (t ), E off (t ) over time. The device is cooled by a thermal network with a given thermal impulse response z th (t ), connected to the ambient temperature T 0 . The device's junction temperature is given by (1) where ' * ' represents a convolution between the thermal impulse response z th (t ) and the losses in the transistor. Transformed to the frequency domain, the convolution becomes a multiplication of the thermal impedance Z th ( jω) and the loss spectrum, leading to the temperature spectrum The goal is to obtain the discrete Fourier series of the device temperature where each component corresponds to the temperature spectrum at discrete frequencies ω 0 k. The underlying assumption is that any sinusoidal current will only cause harmonics at integer multiples of its own frequency. This assumption allows to express each Fourier coefficient of the temperature as a linear combination of the dissipated power and thermal impedance. Only a finite set of N coefficients needs to be considered, as the devices thermal network acts as a low pass and higher frequency coefficients will inevitably tend towards zero. The procedure will produce a system of linear equations for this finite set of coefficients, involving the thermal network, conduction losses and switching losses of the device. By modelling each of these in the frequency domain, including their temperature dependence, a solution for the device temperature can be obtained.

B. VECTOR NOTATION OF DISCRETE FOURIER SERIES
To solve the set of linear equations, the discrete Fourier coefficients of each series are written as a vector such that the center element corresponds to the average component at zero frequency. Components to the right correspond to positive k and components to the left to negative k. The same representation can be used for other quantities of the inverter, such as the thermal impedance Z th or the device losses P c and E sw . Using this vector notation, and by writing the thermal impedance as a diagonal matrix of the discrete coefficients Z th = diag( Z th ), the discrete form of (2) can be written as

C. THERMAL NETWORK
Each transistor has a thermal path from junction to ambient, which determines the thermal response to the power dissipated in the junction. A Cauer network can be used to model this path, as depicted in Fig. 2. Each of the consecutive RC elements can model individual layers within the device or external layers, such as the electric insulation from the device to the heatsink or the heatsink's resistance to ambient temperature. As described in [16], the transfer function of this network can be expressed as

D. CONDUCTION LOSSES
The basis calculation of the conduction losses is a twodimensional lookup table V ds (i ds , T j ), which returns the device forward voltage for a given temperature and load current. The table is used to obtain the losses for ambient temperature through direct evaluation, and to model a temperature dependent component through fitting of an arbitrary curve. For simplicity of the loss calculation, the device temperature is assumed to be constant within one switching cycle, which is reasonable given that typical switching frequencies are significantly above the thermal networks cut-off frequency. To begin, the device current is expressed as the sum of the modulated output current I L and the superimposed current ripple, as depicted in Fig. 3. The load current can be of arbitrary shape, controlled by the modulating the duty ratio d (t ) of the half-bridge. This duty ratio is expressed as where T on (t ) is the timespan during which the device is conducting current and T sw is the switching period. The amplitude of the current ripple can now be written as Within one switching cycle, as depicted in Fig. 4, the device current at turn-on and turn-off can be expressed as I on = I L − I and I off = I L + I . At first, the conduction losses at ambient temperature are obtained from the device's look-up table by integrating the voltage current product. The current i ds is known from the operating parameters of the inverter, and the voltage is obtained by evaluating the look-up table V ds (i ds , T 0 ) at ambient temperature T 0 , as shown in Fig. 3(a). The integral is found using Simpson's rule, which gives an exact solution for polynomials provided single values at the beginning, center and end of the conduction interval. The calculation is done in two individual sections, to account for a potential zero crossing of the current, as indicated in Fig. 4, marked in blue and yellow. Either one becomes zero if both turn-on and turn-off current are above or below zero. The total losses at ambient are given in (9) at the bottom of this page. To now include the temperature dependence, the lookup table is fitted to a second-order temperature-dependent polynomial, with arbitrary functions f n (i ds ) for the known device current, such that (9), the linear and quadratic dependencies can be formulated as P c,1 and P c,2 , proportional to the linear or square of the temperature. The conduction losses can then be expressed as From this expression in the time domain, the spectrum of the conduction losses can be obtained through a Fourier transform, yielding Using the vector notation described in subsection II-B, the coefficients of the conduction loss spectrum can be formulated as where P c,T is a Toeplitz matrix 1 built from Multiplication of a vector with a Toeplitz matrix carries out the convolution of two discrete Fourier series. P c,T still contains the unknown junction temperature resulting from the second-order terms. This issue is addressed by assuming T j = T 0 for the initial calculation of P c,T and T j , and a subsequent second iteration with the resulting T j . Further iteration is possible to improve the results, but little benefit is found for more than two iterations.

E. SWITCHING LOSSES
The switching losses are the sum of the switching energy dissipated in the device each switching cycle at turn-on and turn-off. Again, due to the low pass behaviour of thermal network, the nearly impulse-like switching energy can be averaged over one cycle, yielding For brevity, this section only covers the derivation of the turnon losses E on (t ). The turn-off losses E on (t ) are obtained using the same methodology and the resulting model is included. The turn-on losses depend on the device current at the beginning of the conduction interval, previously established as I on . Depending on the direction of this current, the device will either experience a hard-switched commutation, dissipating a significant amount of energy, or a soft commutation with little to no energy loss.
PLECS deals with this behaviour through a look-up table E on (I on , T j ) and interpolation between the discrete values. The PLECS loss model for a MOSFET is depicted in Fig. 5(b). For the calculation with harmonic balance, similar as with the conduction losses, the turn-on losses are formulated as an ambient component E on,0 = E on (I on , T 0 ), obtained directly through the loss table, and a part with linear temperature dependence E on,T = f on (I on ) T j . The arbitrary function f on (I on ) is chosen, such that E on (I on , T j ) = E on (I on , T 0 ) + B on (I on ) f on (I on ) T j . (16) As most loss models have zero losses for soft-commutation, the fit f on (I on ) only needs to be correct for positive turn-on currents. To eliminate the wrong values when f on is evaluated at negative turn-on currents, the resulting energy is multiplied with a boolean function blanking any values when the inverter is operating in softcommutation.
The turn-on losses over time can then be expressed as which can be transformed into the frequency domain as Same as with the conduction losses, the vector form of the turn-on loss spectrum is obtained by construction of a Toeplitz matrix E on,T from the temperature dependent spectrum E on,T ( jω), resulting in The same steps can be followed to obtain the expression for the turn-off losses, using the turn-off current I off , its corresponding look-up table and boolean function B off . With all loss components modelled in vector format, the initial equation for the junction temperature is solved.

F. SOLVING FOR THE JUNCTION TEMPERATURE
Using the vector notation, (2) can be rewritten using the conduction losses (13), turn-on (20) and turn-off losses (21), resulting in The terms proportional to T j can be collected and T 0 subtracted to obtain Subsequent subtraction and inversion of these terms yield the solution for the increase in junction temperature. Fig. 6 shows the obtained temperature and losses for the high-side switch of the inverter specified in Table 1    calculation are uploaded in conjunction with this paper [18]. The PLECS simulation requires an active PLECS standalone license, the harmonic balance method is implemented in python3, which is open source. The semiconductor loss models are available at the manufacturer's website.
Already a very low number of harmonics N = 16 yields good accuracy in power losses and temperature. For a higher number of harmonics the calculated temperature approaches the shape of the simulated reference, but at a slight offset. The difference in average temperature can be explained by the slightly lower peak conduction losses. Still, the rest of the curve matches well with simulation results. The switching losses also show very good agreement, with a slight improvement at larger N for the discontinuous sections, where the device enters soft-commutation. To show how the proposed method can be used for fast evaluation and design of inverters, the next section will explore a multitude of operating points for a three-phase inverter, again using PLECS as reference, for both accuracy as well as computation time.

III. DESIGN SPACE EXPLORATION OF A THREE-PHASE SPWM GRID INVERTER
To demonstrate the speed and accuracy of the proposed method, a three-phase grid-inverter, as depicted in Fig. 7, is analysed by search through a vast design space. This space, given in Table 3, spans different load currents, switching frequencies, as well as different filter sizes, yielding a total of 1014 possible combinations. As an application example, the load current is stepped through the profile of the European Efficiency for solar inverters, given by η euro = 0.03η 5% + 0.06η 10% + 0.13η 15%

TABLE 3. Design Space and Runtime of Simulation
where η xx is the inverter efficiency at the indicated percentage of maximum load current. For simplicity, the inductive component is treated as lossless, but its impact on the semiconductor losses is considered. The results obtained through harmonic balance are compared to simulations done in with the PLECS, serving as reference for state-of-the-art steady-state-solvers. Both the harmonic balance method and PLECS use the same look-up table for the device losses, available from the device manufacturer. Thus any discrepancy can be entirely addressed to the proposed method.
For now, both simulation and calculation do not include dead-time, which can become a noticeable contributor to losses at higher switching frequencies [19]. Next to additional losses, dead-time leads to distortion in the output voltage of the inverter [20], which would result in different current ripple and modulation of the device, hindering a direct comparison between calculation and simulation.

A. SIMULATION AND CALCULATION SETUP
Within PLECS, the steady-state analysis tool is used to obtain the device temperature and losses of a single grid cycle. For a fair comparison, both PLECS and the harmonic balance method are run on the same machine (Intel Core i5-7500), as single threaded applications. In a practical scenario, paralleled computing can be used to significantly reduce the absolute computation time of both methods. To evaluate each method's computation time, and not the machine's disk speed, both methods are run without saving of the time domain waveforms, as this significantly affects the time spent by each method. Moreover, the time PLECS takes is extremely dependent on tolerances and other settings. Thus, PLECS is kept at the default parameters preset after installation. The analysis settings and scripts are included within the attached model files.

B. RESULTS
The runtime of each program is listed in Table 3, with PLECS taking 18 min 23 s to analyze all 1014 possible circuit variations. In contrast, harmonic balance run at N = 32 harmonics finishes after 11 s. For slightly higher accuracy, the number of harmonics can be increased to N = 128, which leads to a computation time of 1 min 10 s, still considerably faster than PLECS.
The calculation results are shown in Fig. 8. The upper set of plots show the average junction temperature for load currents from 0.8 A to 8 A over one fundamental cycle obtained with harmonic balance. To quantify how close the calculation is to simulation, the plots below show the rms error in temperature, calculated by (26) For all operating points, the error stays well below 1 • C. Fig. 9 shows the temperature, conduction and switching losses at a load current of 16 A and their respective rms error to PLECS. Again, the error stays below 1 • C and 1 W respectively. At low inductance and low switching frequency, the high current ripple leads to increased conduction losses. Increasing the switching frequency reduces the conduction losses at first, at a cost of increased switching losses. By further increasing the switching frequency, the feedback effect mentioned in section I sets in: The conduction losses begin to rise, as the increased switching losses heat the device, resulting in higher on-state resistance, and thus larger conduction losses.
For real-world verification, the following section will compare the proposed harmonic balance method to a three-phase inverter, where the junction-temperature of each device can obtained through a high-resolution measurement of the drainsource voltage during conduction.

IV. EXPERIMENTAL RESULTS
The harmonic-balance method is further validated by experiment with a 400 V three-phase SiC inverter, depicted in Fig. 10. The setup parameters are listed in Table 4. The temperature of the transistors is obtained by measurement of the drain-source voltage during conduction. The detailed setup is described in [21]. The measurement is limited by a  clamp circuit to obtain sufficient resolution. In combination with a measurement of the load current, the device temperature is obtained through a look-up table. To cycle the device temperature, the inverter supplies a 50 Hz sinusoidal load current with varying amplitude to a 10 resistive load. The amplitude is stepped from 5 A to 10 A in 5 s intervals. The resulting temperature cycle is depicted in Fig. 11, alongside the calculation using harmonic balance. Due to the vast difference in load cycle and output frequency, the matrices become immense and calculation for one cycle takes FIGURE 11. Measured temperature cycle compared to calculation using harmonic balance. The calculation result has an offset and amplitude error of about 0.5 • C. Due to imprecise timing of the load pulse generated by the microcontroller, the pulse duration slightly deviates from the programmed 5 s interval. As a result, the early pulses slightly overshoot the calculated response, as the device has less time to cool down. The delayed pulses, on the other hand, undershoot the calculated response, as the device is given slightly more time to cool down. several seconds. The matrices now contain frequency components up to several 100 Hz at a step size of 100 mHz, which significantly slows down the matrix inversion. Still, the calculation shows good agreement with the measured device temperature. The mismatch between measurement and calculation is within of 2 • C. Next to possible discrepancies between device and the manufacturer's loss model, the neglect of dead-time losses does affect the calculation result.

V. CONCLUSION
The calculation of device temperature using harmonic balance offers a fast, precise alternative to iterative steady-state solvers, which includes temperature dependence and saturation effects. The results match simulations obtained with commercial state-of-the-art steady-state solvers. Furthermore, good agreement is shown between calculation and measurement results of a 400 V three-phase SiC inverter. The steady state is calculated given a temperature dependent look-up table for the device losses available from the manufacturer, as well as the device current and duty ratio. For well known topologies, these operating parameters are easily written analytically. Otherwise, harmonic balance could serve as an addition to steady-state solvers, producing a solution after the device operating parameters are obtained through a first simulation cycle.
Several aspects of the presented method can be improved for more accurate results and universal application. As previously addressed, dead-time was omitted to simplify the analytic description of the inverter's output voltage and to limit potential error sources. For future work, dead-time losses need to be added as a term proportional to the switching frequency, including load current and temperature dependence. Furthermore, the conduction losses are calculated assuming a triangular output current, as is the case for PWM inverters or buck/boost converters. To extend application of the harmonicbalance calculation method to other topologies, different load current shapes have to be modelled, such as sinusoidal waveforms for resonant converters, or trapezoidal waveforms for active-bridge converters.
In combination with analytic loss models such as [9], it could enable fast, fully analytic evaluation of a power stage solely from the transistor datasheet. Next to thermal design of the inverter, the obtained temperature cycles can serve as input for life-time estimations [22]. State-of-the-art calculation methods for inductor and capacitor losses could be added to optimize the complete system [23], [24]. Beyond design of inverters, a possible application of harmonic balance could be over-temperature protection. Given that harmonic balance only requires standard matrix functions, it could be implemented on capable DSPs used to control the inverter. As the calculation produces the inverter's steady-state response, the controller could predict over-temperature of the transistors, before their thermal network reaches the undesired operating point and limit the load current to prevent this.
Harmonic balance performs best at small timescales, such as 50 Hz to 60 Hz grid cycles, up to several seconds for step responses of the thermal network, due to the required matrix inversion. For longer periods, such as drive cycles for electric vehicle inverters, the system matrix becomes extremely large, as it would contain frequency components in the kHz range for the load current, down to components in the mHz range for the driving profile. As the computation time for matrix inversion grows exponentially with size, iterative methods would become the faster option. Further research is necessary to see if simplified, reduced timescale models would provide sufficient accuracy.