Saturation with hysteresis

Hysteresis are necessary when some numerical noise appears in simulation and perturbs the saturated signal to a value close to one of the (xmin, xmax) limits. In this case, the solver will lose a lot of time switching from the “saturated” case to the “free” case, and the benefit of discontinuity handling on simulation CPU time will be lost.

Figure 16: saturation with hysteresis

The example of saturation can be reused to demonstrate the use of a hysteresis.

The main modification in the code concerns the settings of the ic region. The part “flag==0” of the code contains 2 different cases.

• the first case corresponds to the first time the submodel is called. In such a case, the ic element is set as if there were no hysteresis. This case is detected by the test "if (ic[0] == -2)". To make it work, the ic[0] element must obviously be initialized to -2 (or any other value not used for discontinuity region) in the initialization function of the submodel, which is not described in the following piece of code.

• the second case covers the rest of the simulation, and sets a new value to ic[0] only when the x value is out of a hysteresis region, i.e. out of the intervals [xmin, xmin+hyst] and [xmax-hyst, xmax].

Copy

void saturation_(int *n, double *y, double *x, double rp[3],
                 int ic[1], int *flag)
{
   int logi;
   double xmax, xmin, hyst;

   xmax = rp[0];
   xmin = rp[1];
   hyst = rp[2];  
   logi = 0;

   /* Discontinuity restart: check current state. */

   if(*flag == 0)
   {
      if (ic[0] == -2) /* Initialization of ic : no hysteresis */
      {
         if (*x > xmax) /* Limited to higher value. */
         {
            ic[0] = 1;
         }
         else if (*x < xmin) /* Limited to lower value. */
         {
            ic[0] = -1;
         }
         else /* Unlimited. */
         {
            ic[0] = 0;
         }
      }
      else /* Change ic only if not in hysteresis region. */
      {
         if (*x > xmax) /* Limited to higher value. */
         {
            ic[0] = 1;
         }
         else if (*x < xmin) /* Limited to lower value. */
         {
            ic[0] = -1;
         }
         else if ((*x <= xmax – hyst) || (*x >= xmin + hyst))
	      {
            ic[0] = 0;
         }
      }
   }

   /* Calculation part. */

   if (ic[0] == 1) /* Limited to higher value. */
   {
      *y = xmax;

      /* Check if becoming free of upper limit. */
      logi = (*x <= xmax - hyst);
   }
   else if (ic[0] == -1) /* Limited to lower value. */
   {
      *y = xmin;

      /* Check if becoming free of lower limit. */
      logi = (*x >= xmin + hyst);
   }
   else /* Unlimited. */
   {
      *y = *x;

      /* Check if becoming limited at upper value ... */
      logi = (*x > xmax) || (*x < xmin);
   }

   if (logi) /* logi is true (1) when using wrong equation. */
   {
      /* Send a warning to the solver. */
      disloc_(&logi);
   }
/* <<<<<<<<<<<<End of Calculation Executable Statements. */
}

Note that the hysteresis is not compulsorily stuck on one side of a limit (e.g. xmin) and can also be centered on the limit, as on the following figure. Note that the hysteresis will have an impact on the submodel results. On the saturation example, the result of the saturation can then be greater than xmax and equal to xmax + hyst / 2.

Figure 17: saturation with hysteresis centered on xmaxsaturation with hysteresis centered on xmax