CODEWITHSUNDEEP
modelicadymolaopenmodelica
EngineerInProgress
EngineerInProgress

Reputation: 109

Dymola/Modelica - How can I calculate the mean of a signal only in the steady state?

I’m struggling with calculating the mean of a continuous signal in my Modelica model. The issue is that I’m only interested in calculating the mean of the steady state. For instance, in the model (Minimal, Reproducible Example) attached I’d like to calculate the mean between 40 and 60 seconds. Therefore, firstly I’ve tried to use the Modelica.Blocks.Math.ContinuousMeanbut the t_0 is protected, so I wrote a similar code in my personal library (it’s also attached). However, there is an error and I can’t figure out what is happening. Could someone shed some light on this simple issue? Thanks in advance!!!

Model:

model MeanExampleStackOverFlow

  parameter Modelica.SIunits.Time period=3;
  final parameter Real angularFrecuency=(2*Modelica.Constants.pi)/period;
  parameter Real amplitude=2;
  Real signal;
  Real mean;

  PersonalLibrary.Utilities.ContinuosMeanStartTime meanBlock(t_0=40);

equation 

  signal = amplitude*sin(angularFrecuency*time) + log(time*100 + 1);
  meanBlock.u = signal;
  meanBlock.y = mean;

end MeanExampleStackOverFlow;

Block saved in PersonalLibrary:

block ContinuosMeanStartTime "Similar block as Continuos Mean in MSL but with unprotected Start time (t_0)"
  extends Modelica.Blocks.Icons.Block;
  parameter Modelica.SIunits.Time t_eps(min=0.0) = 1e-7 "Mean value calculation starts at startTime + t_eps";
  parameter Modelica.SIunits.Time t_0(min=0.0) = 0 "Start time";

  Modelica.Blocks.Interfaces.RealInput u "Noisy input signal";
  Modelica.Blocks.Interfaces.RealOutput y "Expectation (mean) value of the input signal";
protected 
  Real mu "Internal integrator variable";
initial equation 
  mu = u;
equation 
  der(mu) = noEvent(if time >= t_0 + t_eps then (u - mu)/(time - t_0) else 0);
  y = noEvent(if time >= t_0 + t_eps then mu else u);

end ContinuosMeanStartTime;

Finally, here there are some results that show what I'd like to achieve (it has been done editing the results with other software). Best regards.

Note: By the way, is there any other way to access to the protected parameter in MSL instead of writing a new block?

enter image description here

Upvotes: 1

Views: 381

Answers (2)

MMeissner
MMeissner

Reputation: 61

When I run your code, I get an error that a singularity is present at t=40.

I think you basically have a numerical problem at t= t_0 + t_eps, because you then suddenly have a differential equation on mu with eigenvalue of time-t_0, which right after the if-condition becomes true is t_eps and hence very small, which is a problem for the solver.

Here are some ways to solve that:

  • use another solver (for me DASSL did not run through, but LSODAR worked), some might work better for your specific numerical challenge
  • increase t_eps (for me the model runs for t_eps > 1e-6 or so with DASSL), so the eigenvalue at time= t_0 + t_eps becomes bigger
  • remove the noEvent (I am not quite sure why this solves the issue, my guess is that, then the event triggers a reset of DASSL, and reduces its step-size, so the fast eigenvalue can be handled. This also would explain why the original model Modelica.Blocks.Math.ContinuousMean does work at time = t_0 + t_eps, because here if-statement becomes true right after the start of the simulation, where DASSL is in a similar state as after an event and uses very small step-sizes. But these are only educated guesses.)

The other thing with your code is, that if you want to create the green line, you have to change the output of your ContinuosMeanStartTime block to

y =  mu;

otherwise you will get the signal as your mean before time= 40 + t_eps.

Upvotes: 0

Hans Olsson
Hans Olsson

Reputation: 12517

You cannot use t_eps to achieve that at the moment, it's just to guard against division by zero.

If t_0 were exposed it would be possible change that.

Upvotes: 1

Related Questions