Modeling piecewise linear diodes (discontinuous behaviors) with Dyad?

Hello,

In the light of Dyad announcement, and perhaps as a continuation of the Ideal Diode and Ideal Switch instability problem or Diode model not working as expected with ModelingToolkit discussions, I’d like to know how Dyad (or I guess just MTK down the road) deals with the following case which isn’t well managed by the present OpenModelica Compiler (OMC): when a particular parameter value makes the problem singular because it changes the order of the ODE.

Here is below the Modelica code (taken from this exercise), which I’m assuming can be translated quite directly to Dyad (if the if/then/else piecewise function definition syntax is supported, which I didn’t check). The diode is modeled almost like the Modelica Standart Library diode, a resistance which value depends in the sign of the current.

From the way OMC works (first model compilation for “any” parameter value, then simulation as separate stage for a given parameter value), this models compiles fine and run fine as long as Ron>0. However, the simulations makes a runtime “division by zero” error for Ron=0.

In particular, I’ve seen that Dyad has the structural parameter and final parameter syntax, while Modelica doesn’t have structural. Would structural be the way to model the case Ron=0?

To make things clearer, perhaps I should split my questions into two:

  1. forgetting about the piecewise linear aspect (no more if/else definition, instead replacing the diode by a constant resistor: vd=Ron*i), how does the Dyad toolchain handles the case Ron=0?
  2. about the piecewise linear aspect, is the vd = if i>0 then Ron*i else Roff*i syntax supported?
    • and if yes, is there a zero-crossing event generation (it seems in was not the case in Feb 2023)?

model DiodeRectifier
  import Modelica.SIunits.{Voltage, Current, Resistance, Conductance, Capacitance, Frequency};
  import Modelica.Constants.pi;
  
  Voltage u, vd, vr, vc;
  Current i, ir, ic;
  parameter Resistance R=1;
  parameter Resistance Ron=1e-2 "R closed diode";
  parameter Resistance Roff=1e3 "R open diode";
  parameter Capacitance C=1;
  parameter Frequency f=3;
  
equation
  /*KCL, KVL*/
  vr = vc "R//C";
  u = vd+vc "mesh";
  i = ir+ic "node";
  
  /*Components*/
  u = sin(2*pi*f*time);
  vr = R*ir;
  C*der(vc) = ic;
  /*Diode: long syntax version*/
  /*if i>0 then
    vd = Ron*i;
  else
    vd = Roff*i;
  end if;*/
  /*Diode: compact syntax alternative*/
  vd = if i>0 then Ron*i else Roff*i;

end DiodeRectifier;
2 Likes

Yeah I have a solution coming here, it’s actually not symbolic but instead to build solvers capable of handling the switching behavior. So it’s an MIT project for DifferentialEquations.jl, but something that will be targeted through Dyad. Give it about a year though since it’s a bit more researchy, but we have a concrete plan on the exact algorithm.

Great news! Let’s see in a year then.
Only, I keep in mind that back in 2001, Michael Tiller wrote in his first Modelica book that modeling discrete events such as mechanical backlash was complicated but was about to get solved in the near future :wink:

3 Likes

The modelica community kept trying to do symbolic things while keeping DASSL. The fundamental issue is DASSL.

2 Likes

OK, so there is ongoing work on a switched system solver then. Good to know.

Some months ago, a nearby researcher @benoit_caillaud pointed me to the work of one of his colleague: the nonsmooth ODE solver Siconos. Very nice contact mechanics simulation videos. I’ve never tried it though. From what I got, the downside is that it is a first order ODE solver.

We have something similar to that going, but that’s a different project. Yeah that drops to first order and makes some other concessions. It’s a bit heavy handed to go all the way to non-smooth, but that is a fun idea. The BVP/NonlinearSolve work is building towards something of that sort though.