Speed Measurement¶

Baseline System¶

There are many applications where we need to model the interaction between continuous behavior and discrete behavior. For this section, we’ll look at techniques used to measure the speed of a rotating shaft. For our discussion here, we will reuse the mechanical example we discussed previously in our discussion of Basic Equations:

model SecondOrderSystem "A second order rotational system"
  type Angle=Real(unit="rad");
  type AngularVelocity=Real(unit="rad/s");
  type Inertia=Real(unit="kg.m2");
  type Stiffness=Real(unit="N.m/rad");
  type Damping=Real(unit="N.m.s/rad");
  parameter Inertia J1=0.4 "Moment of inertia for inertia 1";
  parameter Inertia J2=1.0 "Moment of inertia for inertia 2";
  parameter Stiffness c1=11 "Spring constant for spring 1";
  parameter Stiffness c2=5 "Spring constant for spring 2";
  parameter Damping d1=0.2 "Damping for damper 1";
  parameter Damping d2=1.0 "Damping for damper 2";
  Angle phi1 "Angle for inertia 1";
  Angle phi2 "Angle for inertia 2";
  AngularVelocity omega1 "Velocity of inertia 1";
  AngularVelocity omega2 "Velocity of inertia 2";
initial equation
  phi1 = 0;
  phi2 = 1;
  omega1 = 0;
  omega2 = 0;
equation
  // Equations for inertia 1
  omega1 = der(phi1);
  J1*der(omega1) = c1*(phi2-phi1)+d1*der(phi2-phi1);
  // Equations for inertia 2
  omega2 = der(phi2);
  J2*der(omega2) = c1*(phi1-phi2)+d1*der(phi1-phi2)-c2*phi2-d2*der(phi2);
end SecondOrderSystem;

We will reuse this model by adding an extends clause to our model. This essentially imports everything from the model we are extending from. We’ll talk more about the extends clause later when we discuss Packages. For now, just think of it as copying the contents of another model into the current model.

Recall the solution for the SecondOrderSystem model looks like this:

In this case, we are simply plotting the solution that we computed. But in a real system, we can’t directly know the rotational velocity of a shaft. Instead, we have to measure it. But measurement introduces error and each measurement technique introduces different kinds of errors. In this section, we’ll look at how we can model different kinds of measurement techniques.

Sample and Hold¶

The first type of measurement we will examine is a sample and hold approach to measurement. Some speed sensors have circuits for measuring the rotational speed of the system. But instead of providing a continuous value for the speed, they sample it at a given point in time and then store it somewhere. This is called “sample and hold”. The following model demonstrates how to implement a sample and hold approach to measuring the angular velocity omega1:

model SampleAndHold "Measure speed and hold"
  extends BasicEquations.RotationalSMD.SecondOrderSystem;
  parameter Real sample_time(unit="s")=0.125;
  discrete Real omega1_measured;
equation
  when sample(0,sample_time) then
    omega1_measured = omega1;
  end when;
end SampleAndHold;

Note the presence of the discrete qualifier in the declaration of omega1_measured. This special qualifier indicates that the specified variable does not have a continuous solution. Instead, the value of this variable will make (only) discrete jumps during the simulation. It is not required to include the discrete keyword but it is useful because it provides additional information about the intent of the model that the compiler can check (e.g., making sure we never request the derivative of that variable).

Let us now examine the solution generated by this model:

The important thing to note in the solution is how the measured value is piecewise-constant. This is because the value of omega1_measured is set only when the when statement becomes active. The sample function is a special built-in function that first becomes true at the time indicated by the first argument (0 in this case) and then at regular intervals after that. The duration of these regular intervals is indicated by the second argument (sample_time in this case).

Interval Measurement¶

In the previous example, we weren’t actually making any estimates for the speed, we were simply reporting the value of the variable omega1 only at specific times. In other words, at the moment that we sampled omega1 our sample was completely accurate. But by “holding” our measured value (instead of continuing to track omega1), we introduced some artifact in the measurement.

In the remaining examples, we will focus on techniques used to estimate the speed of a rotating shaft. In these cases, we will never make direct use of the actual speed in our measurement. Instead, we will respond to events generated by the physical system and attempt to use these events to reconstruct an estimate of the rotational speed.

The events that we will be responding to are generated by the discrete elements attached to the spinning shaft. For example, a typical way to produce these events is to use a “tooth wheel encoder”. A tooth wheel encoder includes a gear on the rotating shaft. On either side of the gear, we place a light source and a light sensor. As the gear teeth pass in front of the light source, they block the light. The result is that the signal from the light sensor will include an approximate square wave. The leading edge of these square waves are the events we will be responding to.

The first estimation method we will examine computes the speed of the shaft by measuring the time interval that passes between events. Knowing that these events occur whenever the shaft has rotated by an angle of , we can estimate the speed as:

This technique for speed estimation can be represented in Modelica as:

model IntervalMeasure "Measure interval between teeth"
  extends BasicEquations.RotationalSMD.SecondOrderSystem;
  parameter Integer teeth=200;
  parameter Real tooth_angle(unit="rad")=2*3.14159/teeth;
  Real next_phi, prev_phi;
  Real last_time;
  Real omega1_measured;
initial equation
  next_phi = phi1+tooth_angle;
  prev_phi = phi1-tooth_angle;
  last_time = time;
equation
  when {phi1>=pre(next_phi),phi1<=pre(prev_phi)} then
    omega1_measured = tooth_angle/(time-pre(last_time));
    next_phi = phi1+tooth_angle;
    prev_phi = phi1-tooth_angle;
    last_time = time;
  end when;
end IntervalMeasure;

where tooth_angle represents . Note how tooth_angle is not something the user needs to specify. Instead, the user specifies the number of teeth using the teeth parameter. The tooth_angle parameter is then computed using the value of teeth (note that while we have hand coded the value of here, we’ll learn how to avoid this later in the book when we talk about Constants).

Let’s take a close look at the when statement in this model:

equation
  when {phi1>=pre(next_phi),phi1<=pre(prev_phi)} then
    omega1_measured = tooth_angle/(time-pre(last_time));
    next_phi = phi1+tooth_angle;
    prev_phi = phi1-tooth_angle;
    last_time = time;
  end when;

Here, we use the vector expression syntax used previously in the Bouncing Ball example. Recall that the when statement becomes active if any of the conditions become true. In this case, the when statement becomes active if the angle, phi1, becomes greater than next_phi or less than prev_phi.

Another thing to note is the use of the pre operator throughout the when statement. When an event occurs in a model, there is a chance that the value of some variables may change discontinuously. During an event, while we are trying to resolve what values all the variables should have as a result of the event, the pre operator allows us to reference the value of a variable prior to the event. The pre operator is used in this model to refer to the previous (pre-event) values of next_phi, prev_phi and last_time. The pre operator is necessary because all of these variables are affected by the statements inside the when statement. So, for example, last_time (without the pre operator) refers to the value of last_time at the conclusion of the event while pre(last_time) refers to the value of last_time prior to any event occurring.

Let’s take a look at the speed estimates provided by this approach:

There are two important properties of this estimation algorithm that we can immediately see in these results. The first is that the estimate is unsigned. In other words, we cannot tell from a device like a tooth wheel encoder which direction the shaft is rotating. Also, low speeds and changes in rotational direction can degrade the accuracy of the estimate significantly. The results are also very sensitive the number of teeth involved. If we were to reduce the number of teeth used in our encoder by setting teeth to 20, we’d get very different results.

To understand exactly why the measured signal is so inaccurate, it helps to consider the following plot which shows how the angle, phi1 compares to the angles associated with the adjacent teeth, next_phi and prev_phi.

In this plot, we can clearly see how relatively low speeds and speed reversals create irregular events that introduce significant estimation error.

Pulse Counting¶

The interval measuring technique mentioned above requires hardware that can perform speed calculations on hardware interrupts. Another approach to estimating speed is the count how many events occur within a given (fixed) time interval and use that as an estimate of speed. Using this method, only the summation of events occurs when the events occur and the calculations are deferred to a regularly scheduled update.

Building on the previous examples in this section, the following seems like a natural way to create a model that implements this estimation technique:

model Counter "Count teeth in a given interval"
  extends BasicEquations.RotationalSMD.SecondOrderSystem;
  parameter Real sample_time(unit="s")=0.125;
  parameter Integer teeth=200;
  parameter Real tooth_angle(unit="rad")=2*3.14159/teeth;
  Real next_phi, prev_phi;
  Integer count;
  Real omega1_measured;
initial equation
  next_phi = phi1+tooth_angle;
  prev_phi = phi1-tooth_angle;
  count = 0;
equation
  when {phi1>=pre(next_phi),phi1<=pre(prev_phi)} then
    next_phi = phi1+tooth_angle;
    prev_phi = phi1-tooth_angle;
    count = pre(count) + 1;
  end when;
  when sample(0,sample_time) then
    omega1_measured = pre(count)*tooth_angle/sample_time;
    count = 0 "Another equation for count?";
  end when;
end Counter;

However, there is a problem in this model. Note that there are actually two equations for count. Trying to compile such a model will lead to a situation where there are more equations than variables (i.e., the problem is singular).

So what can we do about this? We need two different equations because the updates to count occur in response to different events. We could try to formulate everything under a single when statement, like this:

when {phi1>=pre(next_phi),phi1<=pre(prev_phi),sample(0,sample_time)} then
  if sample(0,sample_time) then
    omega1_measured := pre(count)*tooth_angle/sample_time;
    count :=0;
  else
    next_phi := phi1 + tooth_angle;
    prev_phi := phi1 - tooth_angle;
    count :=pre(count) + 1;
  end if;
end when;

But this kind of code quickly becomes hard to read. Fortunately, we can address this situation by placing all the when statements in an algorithm section.

The nature of an algorithm section is that it is treated as one single equation for any variables that are assigned within it. This allows multiple assignments to count, for example. When using an algorithm section, it is very important to understand that the order of assignment becomes important. If a conflict should arise (e.g., a variable is assigned two values within the same algorithm section), the last one is the one that will be used. Another thing to note about algorithm sections is that you cannot write general equations. Instead, you must write assignment statements.

In this way, an algorithm section is very much like the way most programming languages work. The statements in the algorithm section are executed in order and each statement isn’t interpreted as an equation, but rather as an assignment of an expression to a variable. The familiarity of assignment statements may make using algorithm sections attractive to people with a programming background who find the otherwise equation oriented aspects of Modelica disorienting and unfamiliar. But be aware that one big reason to avoid algorithm sections is because they interfere with the symbolic manipulation performed by the Modelica compiler. This can result in both poor simulation performance and a loss of flexibility in how you compose your models. So it is best to use an equation section whenever possible.

In our case, there are no significant consequences to using the algorithm section. Here is an example of how the previous estimation algorithm could be refactored using an algorithm section:

model CounterWithAlgorithm "Count teeth in a given interval using an algorithm"
  extends BasicEquations.RotationalSMD.SecondOrderSystem;
  parameter Real sample_time(unit="s")=0.125;
  parameter Integer teeth=200;
  parameter Real tooth_angle(unit="rad")=2*3.14159/teeth;
  Real next_phi, prev_phi;
  Integer count;
  Real omega1_measured;
initial equation
  next_phi = phi1+tooth_angle;
  prev_phi = phi1-tooth_angle;
  count = 0;
algorithm
  when {phi1>=pre(next_phi),phi1<=pre(prev_phi)} then
    next_phi := phi1 + tooth_angle;
    prev_phi := phi1 - tooth_angle;
    count := pre(count) + 1;
  end when;
  when sample(0,sample_time) then
    omega1_measured := pre(count)*tooth_angle/sample_time;
    count :=0;
  end when;
end CounterWithAlgorithm;

The simulated results of this estimation technique can be seen in the following plot:

Again, we see that this approach cannot determine the direction of rotation. With the following plot, we can get a sense of how many events occur within each sample interval:

In general, the higher the count gets in an interval, the more accurate the estimate.

Conclusion¶

This section demonstrates how we can use the when construct to respond to physical events that occur in our system. These kinds of events and the impact they have on our system are just as important as the continuous dynamics we’ve covered previously. The ability to capture and respond to these physical events is an important part of why Modelica is so well suited to model complete systems since those system frequently include both continuous and discrete behavior.