I'm implementing a power electronics circuit simulator. Currently I have all linear components including capacitors/inductors/transformers, switches, and as the only nonlinear component Shockley diode. I'm planning to add BJT and MOSFET models. BJT is simple: two Shockley diodes and two current sources controlled by the current through the Shockley diodes.
However, MOSFET is more complex. The simplest way I could model it is for linear region:
$$ I_D = K \left( (V_{GS} - V_{th})V_{DS} - \frac{V_{DS}^2}{2} \right) $$
and for saturation region:
$$ I_D = \frac{K}{2} (V_{GS} - V_{th})^2 $$
depending on the operating region. At the boundary \$V_{DS} = V_{GS} - V_{th}\$ the two equations agree.
However, although this model is continuous, it is not complete. It ignores subthreshold leakage. According to Wikipedia, it's modeled by
$$ I_D = I_{D0} e^{\frac{V_{GS} - V_{th}}{n V_T}} $$
but this model is both absurd (it claims there is current even if \$V_{DS} = 0)\$ and noncontinuous (at \$V_{GS} = V_{th}\$ it doesn't agree with the other equations although it should for continuity) so I'm planning not to include it, especially since for a power electronics simulator we are mainly interested in large currents, not small currents.
However, one thing I'm troubled by is that in the saturation mode, the current is not dependent on the drain voltage, which would create a constant current source with no parallel resistance, not optimal for simulation (for example two MOSFETs in series with slightly different drain current would be two current sources fighting with each other, resulting in a paradox). So essentially the resistance of the MOSFET is infinite. I think the simulation would be far more stable if I included Early effect in it, to have a finite output resistance. This Early effect according to Wikipedia is modeled by in the saturation region:
$$ I_D = \frac{K}{2} (V_{GS} - V_{th})^2 (1 + \lambda V_{DS} )$$
which would create yet another discontinuity. At \$V_{DS} = V_{GS} - V_{th}\$ (the boundary between linear and saturation region), we have on one hand:
$$ I_D = \frac{K}{2} (V_{GS} - V_{th})^2 (1 + \lambda (V_{GS} - V_{th}) )$$
And on other hand:
$$ I_D = K \left( (V_{GS} - V_{th})^2 - \frac{(V_{GS} - V_{th})^2}{2} \right) $$ $$ I_D = \frac{K}{2} (V_{GS} - V_{th})^2 $$
...which are not equal. To make them equal, I could either modify the Early effect equation for the saturation region:
$$ I_D = \frac{K}{2} (V_{GS} - V_{th})^2 (1 + \lambda (V_{DS} - (V_{GS} - V_{th})) )$$
Or alternatively include Early effect in the linear region too:
$$ I_D = K \left( (V_{GS} - V_{th})V_{DS} - \frac{V_{DS}^2}{2} \right) (1 + \lambda V_{DS}) $$
Both approaches would make the equations continuous. However, I'm uncertain which approach would model a real MOSFET in a better manner.
I can't include noncontinuous equations in the circuit simulation, because the discontinuities could create problems where the simulation is jumping between two operating regions, not finding a good solution, because the two regions are fighting with each other.
So, how to model Early effect in a MOSFET in a continuous manner?
