Can we use Mathematica for electronics design?

Question

A spoiler alert: I intend to answer my own question. I am posting this because I would like to share my excitement regarding the control theory tools in Mathematica, and also to promote further discussion.

I invite any and all to criticize what I post, correct any errors, and add to the discussion!

Answer 1

An example of electronics design is the design of amplifiers using operation amplifiers (op amps). One common configuration is the inverting amplifier depicted below in a schematic entered in LTSpice (a free spice tool from Linear Technology).

The amplifier makes use of an LT1365 op amp and is designed for a gain of -2. In this implementation, a 10pF capacitor Cp represents parasitic capacitance as may be found with some implementations. The op amp has a single pole in its transfer function. The effect of the added capacitance is to produce a second pole in the feedback path. The two poles therefore appear in the loop gain and threaten to add 180 degrees to the phase delay turning what was negative feedback into positive feedback -- producing instability and even threatening oscillation. What can be done? One method is to introduce a zero in the feedback path to improve the phase margin. We do that with Cf in parallel with Rf. How to choose its value?

First, we use the control theory tools to develop a transfer function model for the closed-loop system:

(* convenient shortcuts *)

(* circuit impedances in the s-domain *)
xl[l_] := s l; xc[c_] := 1/(s c); par[z1_, z2_] := (z1 z2)/(z1 + z2);

(* prefixes *)
k = 1000.; M = 1.*^6; u = 1.*^-6; p = 1.*^-12;

(* nominal values *)
nominals = {cp -> 10 p, ri -> 100 k, cf -> 0 p, rf -> 200 k, 
   av -> -4500., pole -> 2 Pi 30000.};

(* the op amp has a one-pole roll off at 30kHz *)
lt1365 = av pole/(s + pole);

(* derive a closed-loop transfer function *)

(* node currents sum to zero into the inverting input node with \
voltage v *)
eq1 = (vin - v)/ri + (0 - v)/xc[cp] + (vout - v)/par[rf, xc[cf]] == 0;

(* equation for op amp gain *)
eq2 = vout == lt1365 v;

(* determine transfer function as vout/vin *)

temp = Eliminate[{eq1, eq2}, v];

tf = (vout /. Solve[temp, vout][[1]])/vin // Simplify;

tfm = TransferFunctionModel[tf, s];

(* nominal response of transfer function *)

nominalTFM = tfm /. nominals // Simplify;

The nominal transfer function has a value of zero for Cf -- that is no compensation. We see the effects of instability in both the frequency and transient response. The circuit is close to oscillation.

SetOptions[BodePlot, FeedbackType -> None,
  ScalingFunctions -> {{"Log10", "dB"}, {"Log10", "Degree"}},
  GridLines -> Automatic, ImageSize -> 300,
  FrameLabel -> {{"Frequency", "dB"}, {"Frequency", "Degrees"}},
  PhaseRange -> {-\[Pi], \[Pi]}, PlotRange -> Automatic];

plot[1] = 
 BodePlot[nominalTFM, {2 Pi 10^4, 2 Pi 1*^8}, 
  PlotLabel -> "Nominal Response"]

(* the response to a negative pulse shows severe ringing *)
(* brought about by the parasitic capacitance *)

out = OutputResponse[
   nominalTFM, -UnitStep[t - 5 u] + UnitStep[t - 10 u], {t, 0, 15 u}];

plot[2] = 
 Plot[out /. t -> tt u, {tt, 0, 15}, PlotRange -> Automatic, 
  Frame -> True, FrameLabel -> {"\[Mu]s", None}, 
  PlotLabel -> "Pulse Response", PlotLegends -> {"Nominal"}]

We want to determine a value of Cf that will quench the ringing but still provide a fast transient response. The root-locus plot is the perfect way to do that.

(* The ringing can be reduced by placing a zero in the feedback using \
cf *)
(* check the root locus plot *)
(* without compensation the systen is almost oscillating *)
(* a value of cf = 0.478p brings the poles to the real axis *)
plot[3] = RootLocusPlot[tf /. cf -> cf1 /. nominals, {cf1, 0, 1 p},
  FeedbackType -> None,
  PlotRange -> {{-.3*^8, .1*^8}, All},
  AspectRatio -> .7, PlotLabel -> "Closed-Loop Poles",
  PoleZeroMarkers -> {"", Automatic, "", 
    "ParameterValues" -> {0, .25 p, .45 p, .478 p, .5 p, .8 p}}]

We can see on the root-locus plot that with no compensation the two poles are nearly on the imaginary axis -- close to oscillation. Using a set of values for Cf we find the value that brings the poles down to the real axis.

We now have a smooth frequency response:

improvedTFM = tfm /. cf -> .478 p /. nominals;

plot[4] = 
 BodePlot[improvedTFM, {2 Pi 10^4, 2 Pi 1*^8}, 
  PlotLabel -> "Improved Response"]

The transient response is just what we want:

out2 = OutputResponse[
   improvedTFM,
   -UnitStep[t - 5 u] + UnitStep[t - 10 u],
   {t, 0, 15 u}];

plot[5] = 
  Plot[out2 /. t -> tt u, {tt, 0, 15}, PlotRange -> Automatic, 
   Frame -> True, FrameLabel -> {"\[Mu]s", None}, PlotStyle -> {Red}, 
   PlotLabel -> "Pulse Response", PlotLegends -> {"Improved"}];

plot[6] = Show[plot[2], plot[5], PlotLabel -> "Pulse Response"]

It is interesting to compare this with an analysis in LTSpice. It will not be exact, Spice will use the actual op amp model, not just a 1 pole roll-off. But it is quite close. (And note that in Mathematica we get the use of root-locus as a tool.)

Uncompensated:

Compensated:

So I hope my fellow Mathematica fanatics find this interesting. For me -- who started out doing this stuff with slide rule and graph paper -- it is fantastic!

EDIT:

In response to the comment by andre314 that the label on the frequency axis on the bode plot should be radian frequency -- I agree. But I would really prefer frequency. I have tried to implement that using ScalingFunctions like this:

fscale = {Log10[#/2/Pi] &, 2 Pi 10^# &};
SetOptions[BodePlot, FeedbackType -> None, 
  ScalingFunctions -> {{fscale, "dB"}, {fscale, "Degree"}}, 
  GridLines -> Automatic, ImageSize -> 300, 
  FrameLabel -> {{"Frequency", "dB"}, {"Frequency", "Degrees"}},
  PhaseRange -> {-\[Pi], \[Pi]}, PlotRange -> Automatic];

However, when I evaluate BodePlot it returns an error that the value of option ScalingFunction is invalid. I suspect that BodePlot does not accept the user-defined value for ScalingFunction as described in the documentation.

Does anyone know how to get BodePlot to plot in frequency rather than radians/sec?

**** This is answered in the comments! BodePlot[tfm[2 Pi s]] plots in Hz ****

EDIT 2:

In continuing to look at this model I have come on an issue that concerns me greatly. To plot the response of the nominal system to a pulse, I converted the transfer function (s-domain) to a TransferFunctionModel and substituted into it the nominal values. I then used OutputResponse to plot the response to a pulse composed of UnitStep functions.

I have now done this another way. I start with the same transfer function and substitute into it the nominal values. I then multiply it by a pulse in the s-domain. That is the response in the s-domain. I then do an inverse Laplace transform to the time domain to obtain the response.

You see the results below. The output does not agree with the output obtained in the first method. But it does agree with the simulation in LTSpice. (The frequency of ringing is slightly different, but so are the op amp models.)

It might be that I have made some error. Or it might that the system is close enough to unstable that a slight difference in the numerical methods show a burst into oscillation (like after the falling edge in plot[2] ) in one method and not in the other.

Here is the work:

(* a negative pulse in the s domain *)
pulseS[t0_, t1_] := (-Exp[-s t0] UnitStep[t0] + 
  Exp[-s t1] UnitStep[t1])/s

(* the same 10us pulse used for OutputResponse *)
pulse = pulseS[5 u, 10 u];

(* the nominal transfer function in the s domain *)
nominalTF = tf /. nominals // Simplify;

(* the output in the s domain *)
outSD = pulse nominalTF // Simplify;

(* the output as a function in the time domain *)
outTD[t_] = InverseLaplaceTransform[outSD, s, t];

plot[7] = 
 Plot[Re[outTD[t u]], {t, 0, 15}, PlotRange -> Automatic, 
  Frame -> True, FrameLabel -> {"\[Mu]s", None}, 
  PlotLabel -> "Pulse Response Using Laplace Transform"]

EDIT 3:

And here is a fix for this problem thanks to the comment from Andre 314:

(* And a fix thanks to Andre314 *)

out3 = OutputResponse[
   nominalTFM, -UnitStep[t - 5 u] + UnitStep[t - 10 u], {t, 0, 15 u}, 
   Method -> {"NDSolve", MaxStepSize -> 10^-8}];

plot[8] = Plot[out3 /. t -> tt u, {tt, 0, 15}, PlotRange -> Automatic, 
  Frame -> True, FrameLabel -> {"\[Mu]s", None}, 
  PlotLabel -> "Pulse Response", PlotLegends -> {"Nominal"}]

The damping appears slightly stronger, but the overall characteristics are quite similar to the Laplace transform method.

I reduced the step size until the output no longer changed.

Here is the result at MaxStepSize->10^-10

EDIT 4: I have posted a summary of the convergence issue to Wolfram Community here.

Answer 2

I've been exporting schematic designs from geda-gaf to Mathematica for many years using its "mathematica" exporter together with code that has a new home at https://github.com/noqsi/gEDAmath. This also works with Lepton-eda.

The notebook there shows how to create linear transfer function models for the kind of analysis you are doing. While this development predates Mathematica's system modeling machinery, it is easy to apply the new functionality here.