Integrator State Port: Why is its use restricted?

Hello.

I'd like to know why the use of the Integrator Block's State Port is restricted (and I have some more questions, see below). Here's the long story:

What I want to do

I'm trying to write a Simulink block that computes the arithmetic mean value of a continuous input signal over some period of time. In general, beginning and end of a period are given by an external trigger signal, but in this minimal working example for simplification I use a Pulse Generator as a trigger signal (constant period every 1 second). Every rising edge marks the beginning of a new period and also the end of the previous one. Additionally, I use a Triggered Subsystem to sample the computed mean value and hold it until the next period is over.

To compute the mean value I integrate the given signal and divide it by the period length (which is equal to 1 here, so I left that out). The integrator is reset at the beginning of every period. The integrator state at the end of a period is the result I'm looking for.

First attempt

My first solution looked like this:

This doesn't work as naively expected because the integrator reset occurs before the Triggered Subsystem can copy the result. So the result is always zero:

I think I understand this behavior, therefore no questions so far ...

Second attempt

My second solution involved the integrator's state port. The Integrator Block documentation states:

"[...] If the block is reset in the current time step, the output of the state port is the value that would have appeared at the block's standard output if the block had not been reset. The state port's output appears earlier in the time step than the output of the Integrator block's output port. [...]"

This doesn't work either. Simulink throws an error:

"State ports can only be used to break algebraic loops or to "hand-off" states between systems. Use the output port rather than the state port of 'model/Integrator' as the source of the signal routed (either by direct or virtual connection) to 'model/Triggered Subsystem' "

The Integrator Block documentation also tells me that the use of the state port is restricted to only two special scenarios, stating:

"When updating a model, Simulink checks that the state port applies to one of these two scenarios. If not, an error message appears."

So my first set of questions is: Why is this an error? Is there a technical necessity? Is there a problem I don't see yet?

Third attempt

One of the mentioned two special scenarios where I'm allowed to use the state port is the case of a self-resetting integrator. So I tricked Simulink to believe that that's what I'm doing. I added two Gain blocks and a Sum block which should not change anything.

Now I don't get any error messages and the model works perfectly as expected (the mean value of a constant 1 signal over 1 second is equal to 1; the result is output for the first time when the first period ends, which is at t=1):

So my second set of questions is: Why does this work, although the former example doesn't? Does this solution have any hidden disadvantages, apart from its obvious ugliness? Is it more difficult for the variable-step solver, the zero-crossing detection or any other part of the simulation? Is there any reason why I shouldn't use this solution? Could I get wrong results in a "less minimal example" context?

And my third, only supplementary question: Is there an even better way to solve the problem? (I already tried manual (trapezoidal rule) integration using memory blocks but that's less accurate and also somewhat ugly)

Answer 1

[Paul]

Only answering the third, supplementary question ....

I think I was able to do what you want without using integrator reset at all. Just let the integrator run and use the triggered subsystem to sample, and store, the integrator output at the boundaries of your intervals and compute the average. I implemented your original block diagram, but without the reset on the integrator. The triggered subsystem is:

The delays and the digital clock all have sample time of -1. Both delays have an initial output of 0, though you might want to use a different value on the clock delay to avoid potential for divide by zero if the trigger hits at t=0. Set the initial output of the output port to 0.

Here's the result I got using a sine wave as the source input.

[creepydog]

Thank you for your answer to my supplementary question, Paul. I appreciate your efforts to help me! Your solution is much simpler than mine (which is good!). However, it uses an open integrator whose state/output will go towards +/– infinity (depending on the input signal). I am not a numerics expert, but generally I try to avoid calculating differences of infinitely growing signals. Of course, this might only be a problem with very long simulation times, so I think your solution is certainly suitable for many cases.

In the meantime I've come up with a solution that uses two integrators in parallel and resets and reads them alternatingly. This unfortunately involves some overhead, so it is not ideal, but at least it avoids the risk mentioned above.

[Paul]

If you really want to reset the integrator at the start of every window, you could modify the above solution as follows, assuming your trigger signal is really an impulse. You may be able to adapt this depending on the actual form of your trigger signal. Here is the top level model:

The signal to be averaged is the sine wave input to the integrator. That business in the upper left generates an example trigger signal that defines the edges of the windows. Note that the integrator resets on the falling edge of the triger and the subsystem triggers on the rising edge. This is key as it allows the subsystem to use the integrator output before it gets reset. The subsystem is basically the same as before, except now there is no need to sample the integrator and take a difference, because the integrator output is always zero at the leading edge. Here is the subsystem.

And here is the result:

Top plot shows the window edges, middle is the average over the interval between edges, and the bottom is the integrator, showing the reset at the start of every window.

[creepydog]

Understood.

Although my minimal example uses a Pulse Generator, my real model uses the output of a Hit Crossing block exactly as in your example.

I had tried a similar approach before by using a Memory block between the Integrator output and the Triggered Subsystem (and resetting the integrator on the rising edge). Both approaches do not integrate over the whole period. Mine doesn't cover the last time step of the period, yours starts late and misses the first time step. Strictly speaking, the result is not correct. Quantitatively that might be negligible if (!) Simulink places some time steps very closely before and after the trigger events, which it kindly seems to do in my case. Example extract of "tout" (trigger at t=1):

                      0.85
                      0.93
         0.999999999999993
                         1
          1.00000000000001
                      1.05
                      1.13

I know you were aware of this restriction because you wrote "assuming your trigger signal is really an impulse"!

But is that behavior guaranteed by Simulink whenever I use a Hit Crossing block output to generate my trigger signal? I don't think so. I just tried this

and got these time steps (ode23tb, default settings):

                      0.85
                      0.93
                         1
                      1.05
                      1.13

and, of course, grossly wrong integration results (2nd plot should stay at exactly 1):

I'm aware, though, that feeding pulses into a Hit Crossing is suboptimal ...

[Paul]

Putting a memory block or delay between the integrator and Triggered Subsystem is not a good appraoch IMO.

Based on the the doc page, the output of the hit crossing block is well defined based on the input and the block parameters. Having said that, I'm not sure what's going on in your example. AFAICT, the only way to get your output from the hit crossing block is

• If the input signal is exactly the value of the offset value after the hit crossing is detected in the specified direction, the block continues to output a value of 1.

I used the pulse generator, bias, and hit crossing and got a nearly identical result as using the sine wave to generate the trigger for a sine wave input to the integrator (using the default solver settings, maybe that has something do with it?). In the result below, y is the output of the Triggered Subsystem. out1 is using the pulse generator and bias input to the hit crossing and out2 is using the sine wave input to the hit crossing.

>> [out1.y.Time(7:13) out1.y.Data(7:13) out2.y.Time(7:13) out2.y.Data(7:13) out1.y.Data(7:13)-out2.y.Data(7:13)]

ans =

   0.650000000000000                   0   0.650000000000000   0.000000000000000  -0.000000000000000
   0.850000000000000                   0   0.850000000000000   0.000000000000000  -0.000000000000000
   1.000000000000000   0.459697693754258   1.000000000000000   0.000000000000000   0.459697693754258
   1.000000000000014   0.459697693754258   1.000000000000014   0.459697693754263  -0.000000000000005
   1.000000000000028   0.459697693754258   1.000000000000028   0.459697693754263  -0.000000000000005
   1.050000000000000   0.459697693754258   1.050000000000000   0.459697693754263  -0.000000000000005
   1.250000000000000   0.459697693754258   1.250000000000000   0.459697693754263  -0.000000000000005

As you can see, the pulse generator hit triggered exactly at 1.000 but the sine wave hit triggered 0.00**14 seconds later. I can only assume that this happens because sin(2*pi*t) hasn't crossed zero at t = 1 due to rounding, or other vagaries in how Simulink hit crossing detection works (maybe in these cases it takes advantage of knowing exactly what the input is to the hit crossing block?) In any case, I'm not terribly surprised by this small and subtle difference. And the output y will be very slightly different because of that small difference in when the hit crossing is detected and the integrator resets.

I'm not sure exactly what you mean when you say that this approach won't integrate over the whole period. Are you concerned about missing time between 1.000000000000000 and 1.000000000000014? If so, I doubt that you can get much better with the constraint that the integrator be reset.

[creepydog]

Did you simulate the Sine and Pulse cases simultaneously in one model?

• That would easily explain why you get nearly identical results. Due to the rounding etc. the Sine case seems to force a time step at 1.000000000000014, which of course applied to the whole model, so the output of the pulse-driven Hit Crossing quickly steps back to 0. But without the Sine case in the same model, it would not have created that time step on its own.

Perhaps there's a little misunderstanding: In my example, the output of the Hit Crossing block is equal to 1 only for one time step (at t=1, t=2, ...) and returns to 0 "immediately", that is: at the next time step (t=1.05, t=2.05, ...), as expected. So the "input exactly equal to offset" case never occurs here (*).

I'm not "concerned about missing time between 1.000000000000000 and 1.000000000000014", of course. However, I'm concerned about missing one time step if that time step is much larger. In my example, the variable step solver chooses to take a large time step immediately after the zero crossing (from t=1 to t=1.05), which is 5% of my averaging period. I cannot neglect that, of course.

I don't know why the solver decides to take "1e-14" steps around the zero crossing in one case (Sine) and not in the other (Pulse). Apparently, there's no guarantee for these negligibly small steps.

But, I'm afraid, our discussion has drifted quite far away from my original questions. ;-)

(*) The "Hit Crossing offset" parameter is set to 0. Pulse Generator parameters: amplitude=1, period=1, pulse width=5, phase delay=0. Bias parameter=-0.1. The input of the Hit Crossing can only be -0.1 or +0.9, but never equal to the offset parameter. I just reimplemented the whole thing in a new model from scratch, on another PC, using R2020a instead of R2018b, using ode45 instead of ode23tb, default parameters. Same result. Changing the simulation time can change the result slightly, of course, as the auto max step size is derived as simtime/50, which is 4/50=0.08 in my case. With pulse width=5[%], the next time step after t=1 occurs at t=1.05 (because the Pulse Generator wants to have a time step there). When I change the pulse width to values larger than 8[%] (which is irrelevant for the rising edge triggered Hit Crossing), then the max step size becomes effective and I get the next time step at t=1.08. Okay, in a larger real world model there might be a lot of reasons for the solver to take smaller steps, but not in the 1e-14 region.

[Paul]

"Did you simulate the Sine and Pulse cases simultaneously in one model?"

As it turns out I did. I thought I had it set up so that the sine trigger signal would have no impact when running with the pulse+bias, but apparently not. I ran with just the pulse+bias input into the hit crossing and got the same result that you did. I played around with some other sources driving the Hit Crossing, and it seems to behave like I want (spike output) if the inputs are continuous at the hit crossing. At this point, after rereading the doc page and trying all of this, I'm not sure what to expect from that block, gurantee or otherwise. Sorry can't be of more help. This is an interesting problem and what you are trying to do shouldn't be this difficult.

[Paul]

Here's a modification that might have more appeal, though I'm sure it has some limitations somewhere. The idea is as follows. Use the hit crossing to do two things:

Trigger the average system as has already been done.

Trigger a separate subsystem that computes a future time (from the time of the hit crossing detection) that you can use to reset the integrator.

This way the trigger only depends on the solver catching the hit crossing and you can control how much of the interval that you're going to miss at the start of the window.

Here's the model:

The Average block triggers on the hit crossing as does the Reset Time block, which is:

All the reset time block does is compute a small offset from the time of the trigger that is the time to reset the integrator. The inital output is set to zero.

When sim clock >= reset time, the integrator resets (rising). As you can see, I made the time offset 1e-8. I'd not be surprised if problems arise if you make that offset too small, but I didn't experiment. Anyway, here's the output:

So far the only odd thing I've noticed is that the solver takes a lot more steps to integrate through some reset times than others. For example, it seems to get through t = 3, pretty quickly:

>> out.tout(45:55)

ans =

     2.482494474052150e+00
     2.682494474052150e+00
     2.882494474052150e+00
     3.000000000000000e+00
     3.000000009999972e+00
     3.000000010000000e+00
     3.000000010000028e+00
     3.050000000000000e+00
     3.099999989999971e+00
     3.199999969999914e+00
     3.399999929999800e+00

But a lot more steps to get through t = 4:

 >> out.tout(55:85)

ans =

     3.399999929999800e+00
     3.599999929999800e+00
     3.799999929999800e+00
     3.999999929999801e+00
     4.000000000000000e+00
     4.000000009999943e+00
     4.000000010000000e+00
     4.000000010000057e+00
     4.000000360001054e+00
     4.000000710002051e+00
     4.000001410004045e+00
     4.000002810008033e+00
     4.000005610016010e+00
     4.000011210031963e+00
     4.000022410063869e+00
     4.000044810127682e+00
     4.000089610255308e+00
     4.000179210510558e+00
     4.000358411021060e+00
     4.000716812042064e+00
     4.001433614084070e+00
     4.002867218168084e+00
     4.005734426336111e+00
     4.011468842672165e+00
     4.022937675344273e+00
     4.045875340688489e+00
     4.050000000000000e+00
     4.070623296557556e+00
     4.141246583115055e+00
     4.282493156230053e+00
     4.482493156230053e+00    

I don't know what's going on here, but I did notice that when I made the pulse generator period 1.1 sec the solver took about the same, small number of steps to integrate through each reset time.

[creepydog]

I still think the approach using the integrator state port (see me original question) should be better (cleaner, more robust ...) than all our workaround solutions that rely on sufficiently small time steps.

Apart from that, I've posted a new question about zero-crossings, partially based on our discussion. See here: https://de.mathworks.com/matlabcentral/answers/553960-number-of-necessary-time-steps-to-handle-a-zero-crossing

Answer 2

[Jonas]

"This doesn't work as naively expected because the integrator reset occurs before the Triggered Subsystem can copy the result. So the result is always zero:"

Maybe you can add a Unit Delay in front of the Integrator reset, such that the triggered subsystems is triggered one sample earlier than the Integrator is reset.

[creepydog]

Thank you, Jonas, for your proposal regarding my supplementary question. I'm afraid, adding a Unit Delay block doesn't work in a Continuous Sample Time environment because it only supports discrete sample time. So it would delay the signal for a whole discrete step, not only for a single time step of the variable step solver.

In case you instead meant adding a Memory block: This only works under some assumptions. See Paul's similar proposal and my reply.

[Jonas]

Ah indeed, my mind was fixed on discrete time simulation. I have too little experience with continuous time simulation to answer your questions.

Good luck!