Experimenting with Diodes and Non-Linear Curve Fitting

How much does a diode's characteristics actually diviate from diode to diode? I measured a small number of diodes, and then recovered the SPICE DC parameters. Due to the limitations of my equipment, I only did this in the DC forward biased region. I'll also describe how to accurately recover the SPICE parameters using numerical methods.

Forward bias SPICE Diode Model at DC

SPICE models a diode using the Shockley diode equation with a series resistor. That is,

$$I = I_s \left(e^{\frac{V_D - I R_s}{n k_b T}}-1\right)$$ Where:

• $I~[A]$ is the current through the diode.
• $V_D~[V]$ is the voltage across the entire diode (i.e. voltage across diode junction and internal series resistance)
• $I_s~[A]$ is the reverse biase saturation current/scale current.
• $R_s~[\Omega]$ is the series resistance.
• $n$ is the emission coefficient.
• $k_b~[\frac{eV}{K}]$ is Boltzmann's constant.
• $T~[K]$ is the absolute temperature of the junction.

Unfortunately, many of these parameters are not given in diode datasheets, let alone any usable tolerances. However, they do usually have a V-I diagram (plot of $V_D$ vs. $I$). Now inverting this equation is a bit tricky... as far as I'm aware, there's not a closed form solution that can be used for curve fitting. However, we can use Newton's method with least squares to get an accurate curve fit.

Non-linear curve fitting (aside)

Re-arranging the model equation in the form of a root equation (note: it's not necessary to do this level of re-arrangement; simply moving I over to the right would have worked, too),

$$0 = I R_s + \ln\left(\frac{I}{I_s} + 1\right) n k_b T - V$$ The first step is to take multiple measurements of $V_D$ vs. $I$, either with a real diode by picking data from a datasheet. Then, write a set of $m$ equations (where $m$ is the number of datapoints): $$\begin{matrix} 0 = I_1 R_s + \ln\left(\frac{I_1}{I_s} + 1\right) n k_b T - V_1 \\ 0 = I_2 R_s + \ln\left(\frac{I_2}{I_s} + 1\right) n k_b T - V_2 \\ \vdots \\ 0 = I_m R_s + \ln\left(\frac{I_m}{I_s} + 1\right) n k_b T - V_m \end{matrix}$$ Newton's method in 1D (1 independent variable) is pretty simple: $$x_{i+1} = x_i + \frac{y_{i}}{y'_{i}}$$ All we need to do is generalize this to multiple variables. This is done by writing a Jacobian matrix: $$\mathbf{J} = \begin{bmatrix} \frac{\partial y_1}{\partial x_1} & \frac{\partial y_1}{\partial x_2} & \ldots & \frac{\partial y_1}{\partial x_j}\\ \frac{\partial y_2}{\partial x_1} & \frac{\partial y_2}{\partial x_2} & \ldots & \frac{\partial y_2}{\partial x_j}\\ \vdots & \vdots & \ddots & \vdots\\ \frac{\partial y_m}{\partial x_1} & \frac{\partial y_m}{\partial x_2} & \ldots & \frac{\partial y_m}{\partial x_j}\\ \end{bmatrix}$$ For $j$ is the number of independent variables. In our case, the independent variables we have are $n$, $R_s$, and $I_s$, so $i = 3$. Now we can't "divide" matrices, but we can multiple by the inverse (note, this derivation is a slight lie because there isn't an inverse for non-square matrices), so Newton's method becomes: $$\vec{x}_{i+1} = \vec{x}_i + \mathbf{J}_i^{-1} \vec{y}_i$$ Rather than compute the inverse of $\mathbf{J}_i^{-1}$, we can re-arrange and solve a system of equations. $$\mathbf{J}_i \vec{x}_{i+1} = \mathbf{J}_i \vec{x}_i + \vec{y}_i$$ Everything on the right is known, so simply compute this and solve the system. The careful reader might want to know what do if $\mathbf{J}_i$ is not square. We can use least squares. This turns the system into a system $\mathbf{A} \vec{x} = \vec{b}$, where now we do have a square matrix. We solve the normal equation: $$\mathbf{J}_i^T \mathbf{J}_i \vec{x}_{i+1} = \mathbf{J}_i^T \mathbf{J}_i \vec{x}_i + \mathbf{J}_i^T \vec{y}_i$$ However, this is poorly conditioned and if you use a computer to try and solve this, it's likely your answer will be off due to rounding. We can improve our answer using QR factorization (note: here, I'm using the "square" version of QR factorization, where $\mathbf{R}$ gets zero-padded). $$\mathbf{Q_i R_i} = \mathbf{J}_i$$ $$\mathbf{R}_i^T \mathbf{Q}_i^T \mathbf{Q}_i \mathbf{R}_i \vec{x}_{i+1} = \mathbf{R}_i^T \mathbf{Q}_i^T \mathbf{Q}_i \mathbf{R}_i \vec{x}_i + \mathbf{R}_i^T \mathbf{Q}_i^T \vec{y}_i$$ Due to some great properties of $\mathbf{Q}_i$ (it's orthonormal), $$\mathbf{R}_i^T \mathbf{R}_i \vec{x}_{i+1} = \mathbf{R}_i^T \mathbf{R}_i \vec{x}_i + \mathbf{R}_i^T \mathbf{Q}_i^T \vec{y}_i$$ However, everything has the same factor of $\mathbf{R}_i^T$ at the front, so we can remove this: $$\mathbf{R}_i \vec{x}_{i+1} = \mathbf{R}_i \vec{x}_i + \mathbf{Q}_i^T \vec{y}_i$$ Which along with being numerically much better (it's condition number is $O(\kappa)$ instead of $O(\kappa^2)$ for the naive implementation), is also easier to solve because $\mathbf{R}_i$ is upper triangular (granted, the QR factorization is not necessarily easy, so it's not actually much easier to solve if you have to write the QR factorization algorithm). Incidentally, the number of required data points is $m \ge j$.

Implementation

The Jacobian matrix for our diode model is: $$\mathbf{J}_i = \begin{bmatrix} \ln\left(\frac{I_1}{I_{s,i}} + 1\right) k_b T & I_1 & \frac{-I_1 n_i k_b T}{I_1 I_{s,i} + I_{s,i}^2}\\ \ln\left(\frac{I_2}{I_{s,i}} + 1\right) k_b T & I_2 & \frac{-I_2 n_i k_b T}{I_2 I_{s,i} + I_{s,i}^2}\\ \vdots & \vdots & \vdots\\ \ln\left(\frac{I_m}{I_{s,i}} + 1\right) k_b T & I_m & \frac{-I_m n_i k_b T}{I_m I_{s,i} + I_{s,i}^2} \end{bmatrix}$$ Even though I've explained the mathematical theory for how to perform a non-linear curve fit, the code I used actually uses the SciPy leastsq function. This does effectively what I described, though it may do it in a slightly different fashion. The only inputs are a function for evaluating $\vec{y}_i$ and $\mathbf{J}_i$ (technically you don't need to provide a function for evaluating $\mathbf{J}_i$, but it will converge quicker if you do).

# Yes, I know. I'm lazy.
from numpy import *
from scipy.optimize import leastsq

def diode_res(args, T, V, I):
    '''
    Shockley diode equation with series resistor residual

    @param n emission coefficient [unitless]
    @param Rs series resistance [ohms]
    @param Is saturation current [A]
    @param T reference temperature [K]
    @param V voltage across diode [V]
    @param I current through diode [A]
    
    @return I*Rs+log(I/Is+1)*n*k*T-V
    '''
    n, Rs, Is = args
    kb = 8.617332478e-5
    return I*Rs+log(I/Is+1)*n*kb*T-V

def jac_diode(args, T, V, I):
    '''
    Jacobian of the diode equation.
    Independent vars: n, Rs, Is
    
    @param n emission coefficient [unitless]
    @param Rs series resistance [ohms]
    @param Is saturation current [A]
    @param T reference temperature [K]
    @param V voltage across diode [V]
    @param I current through diode [A]

    @return Jacobian matrix transposed
    '''
    n,Rs,Is = args
    
    kb = 8.617332478e-5
    res = zeros([3,len(V)])
    res[0,:] = log(I/Is+1)*kb*T
    res[1,:] = I
    res[2,:] = -I/(I*Is+Is**2)*n*kb*T
    return res

def diode_coeffs(V, I, T):
    '''
    Find diode coefficients given V vs. I at a specified T.
    @param V
    @param I
    @param T in K
    @return n, Rs, Is
    '''
    return leastsq(diode_res, array([1,0,1e-14]),args=(T,V,I), Dfun=jac_diode, col_deriv=True,xtol=1e-15)[0]

Experimental Setup

The setup couldn't be any easier:

Figure 1. Test ciruit

My power supply (Korad KA3005D-3S) is capable of measuring $V_s$ to $\pm 0.01V$, and my multimeter (Extech MN36) is capable of measuring current to $\pm 0.01 mA$ or $\pm 0.1 mA$, depending on the range. You may immediately recognize that I am not using a second meter to simultaneously measure $V_D$, instead relying on the power supply meter! This means I will have to take the burden voltage of the internal shunt resistor into account. The reason I don't use a second multimeter to measure $V_D$ is because I don't have one. It's not an ideal setup, but it does work. The power supply was set to the following voltages for every diode: 0.6V, 0.7V, 0.8V, 1V, 1.5V, 2V, 2.5V, and 3V.

The diodes tested were the 1N4148 made by Fairchild Semiconductor.

Data

The data was measured at $T = 19^\circ C \pm 1^\circ C$. Table 1 shows the measured currents at the various voltages for each diode tested. The measured lead resistance (using a 4-wire measurement with the power supply panel meter and multimeter) was $17.319 \Omega$.

Table 1. Measured currents [mA]

	0.6V	0.7V	0.8V	0.9V	1V	1.5V	2V	2.5V	3V
Diode 1	0.52	2.31	5.58	9.82	14.46	39.8	66.5	93.8	121
Diode 2	0.56	2.43	5.75	9.98	14.58	39.7	66	93	119.8
Diode 3	0.54	2.4	5.7	9.91	14.54	39.8	66.1	93.3	120.3
Diode 4	0.49	2.27	5.54	9.75	14.38	39.98	66.1	93.3	120.3
Diode 5	0.54	2.38	5.69	9.95	14.6	40	66.7	94	121.1
Diode 6	0.57	2.46	5.8	10.06	14.7	40	66.5	39.7	120.8
Diode 7	0.52	2.32	5.61	9.85	14.5	39.9	66.4	93.6	120.9

Analysis

Figure 2 shows the reconstructed I-V curves, demonstrating that the technique appears to work.

Figure 2. Reconstructed I-V curves at $25^\circ C$.

Figure 3 shows the relative deviation of the physically tested diodes vs. the manufacturer's SPICE model. This indicates there is likely some sort of systematic error in the test setup. Possible reasons for this error include assuming the given constants are constant vs. temperature, the accuracy issues associated with using a low-resolution panel meter for measuring voltage, and doing a round-about reconstruction for $V_D$ using a constant measured shunt resistance. Also, because these are relative differences small errors become amplified for small input currents. Ignoring the relative issues at really low currents, the deviations really are quite small, being withing 10% down to 0.1 mA.

Figure 3. Relative deviation from manufacturer's SPICE model.

Figure 4 comparing just the physical diodes vs. the mean measured parameters. This removes most of the systematic contributions to deviation, and the batch I tested was within 10% deviation all the way down to $0.6~\mu A$, due to one outlier. Without this outlier the error region is much lower at 2 nA.

Figure 4. Relative deviation from mean measured parameters.

Table 2 shows statistics on the measured model parameters. This shows there are significant relative deviations in $R_s$ and $I_s$. Deviations are due to physical placement of the leads (may not have been the same everytime, leading to different lead resistances), as well as the issues due to the current shunt resistance. Self heating could also play a factor as the measured currents went up to 120 mA, leading to ~111 mW dissipated in the diode.

Table 2. Statistics on measured diode parameters.

	$n$	$R_s$	$I_s$
mean	1.7875	$0.636~\Omega$	1.081 nA
std. deviation	0.022	$0.0913~\Omega$	0.1924 nA
relative error ($\frac{std. dev.}{mean}$)	1.23%	14.4%	17.8%

Conclusion

This test does provide a convenient rough bound on what should be relatively easy to achieve in practice. The parameters are certainly within an order of magnitude of each other. Unfortunately, it's hard to determine if most of the deviation errors in the diode model parameters are due to experimental setup or intrinsic to the diodes. Additionally, 7 test samples is hardly what could be called a statistically significant quantity. The test setup was quite limited because I lacked a second DMM to directly measure the diode voltage, instead relying on a fixed shunt resistance that was measured before-hand. Comparisons to the datasheet and manufacturer SPICE model had some sort of systematic error. This could have been associated with assuming the model parameters are constant vs. temperature, where in reality they could drift. There was significant self heating and getting an accurate junction temperature was not possible.