In the previous Interpolation example, we saw external C functions could be used to manage and interpolate data. In this section, we will explore how to integrate C code for an embedded controller into a Modelica model.
When building mathematical models of physical systems in Modelica, it is sometimes useful to integrate (external) control strategies into these models. In many cases, these strategies exist as C code generated for use with an embedded controller. This example revisits the topic of Hysteresis, but with an interesting twist. Instead of implementing the hysteresis behavior in Modelica using discrete states, we will implement it in an external C function. Although this example is extremely simple, it demonstrates all the essential steps necessary to integrate external control strategies.
Let’s start by looking at our “physical model”. In this case, it is essentially just a reimplementation of the model we created previously during our discussion on Hysteresis. Our revised implementation looks like this:
model HysteresisEmbeddedControl "A control strategy that uses embedded C code"
type HeatCapacitance=Real(unit="J/K");
type Temperature=Real(unit="K");
type Heat=Real(unit="W");
type Mass=Real(unit="kg");
type HeatTransferCoefficient=Real(unit="W/K");
parameter HeatCapacitance C=1.0;
parameter HeatTransferCoefficient h=2.0;
parameter Heat Qcapacity=25.0;
parameter Temperature Tamb=285;
parameter Temperature Tbar=295;
Temperature T;
Heat Q;
initial equation
T = Tbar+5;
equation
when sample(0, 0.01) then
Q = computeHeat(T, Tbar, Qcapacity);
end when;
C*der(T) = Q-h*(T-Tamb);
end HysteresisEmbeddedControl;
Let’s take a closer look at the equation
section:
equation
when sample(0, 0.01) then
Q = computeHeat(T, Tbar, Qcapacity);
end when;
C*der(T) = Q-h*(T-Tamb);
end HysteresisEmbeddedControl;
The function computeHeat
is called every 10 milliseconds to
determine the amount of heat to be used. As we will see in a moment,
the controller implements a “bang-bang” control strategy. That means
it flips between zero heat generation and full heat generation. As we
saw in our previous discussion on Hysteresis, this kind of
approach can lead to “chattering”. For that reason, we put the
calculation of Q
inside a when
statement that is executed
every 10 milliseconds. This 10 millisecond interval is essentially
implementing the behavior of what is normally called a “scheduler”
which decides when different control strategies are executed.
The Modelica function computeHeat
used to determine how much heat
should be applied to the system at any given time is implemented as:
impure function computeHeat "Modelica wrapper for an embedded C controller"
input Real T;
input Real Tbar;
input Real Q;
output Real heat;
external "C" annotation ...
end computeHeat;
Note the presence of the external
keyword. This time, however, we
don’t see the name of the external C function like we did in the
previous examples. This means that the external C function has
exactly the same name and arguments as its Modelica counterpart.
Looking at the source code for the C function, we see that is the
case:
#ifndef _COMPUTE_HEAT_C_
#define _COMPUTE_HEAT_C_
#define UNINITIALIZED -1
#define ON 1
#define OFF 0
double
computeHeat(double T, double Tbar, double Q) {
static int state = UNINITIALIZED;
if (state==UNINITIALIZED) {
if (T>Tbar) state = OFF;
else state = ON;
}
if (state==OFF && T<Tbar-2) state = ON;
if (state==ON && T>Tbar+2) state = OFF;
if (state==ON) return Q;
else return 0;
}
#endif
In other words, we can save ourselves the trouble of specifying how the input and output arguments of our Modelica function map to those of the underlying C function by defining them in such a way that no mapping is truly necessary.
Note the presence of the static
variable state
in the C
implementation of computeHeat
. The use of the static
keyword
here indicates that the value of the variable state
is preserved
from one invocation of computeHeat
to another. This kind of
variable is quite common in embedded control strategies because they
need to preserve information from one invocation of the scheduler to
the next (e.g., to implement hysteresis control, as we see here).
The presence of this static
variable is a significant problem
because it means that the function computeHeat
can return
different results for the same input arguments. Mathematically
speaking, this is not a true mathematical function since a
mathematical function can only depend on its input arguments. In
computer science, we say such a function is “impure”. This
means that each invocation of the function changes some internal
memory or variable which affects that value returned by the function.
Given that such impurity is implemented in embedded control strategies by design, we need to be careful when using them in a mathematically oriented environment like Modelica. This is because the Modelica compiler assumes, by default, that all functions are pure and side effect free and the presence of impurity or side effects can result in very inefficient simulations, at best, or completely erroneous results, at worst.
These problems occur because the underlying solvers must compute many “candidate” solutions before they compute the “real” solution. If generating candidate solutions requires the solver to invoke functions with side effects, the solver will be unable to anticipate the effects triggered by the changes to variables it is not aware of.
It is for precisely this reason that the impure
qualifier is
applied to the definition of computeHeat
:
impure function computeHeat "Modelica wrapper for an embedded C controller"
input Real T;
This informs the Modelica compiler that this function has side effects or returns a result that depends on something other than its inputs and that it should not be invoked when generating candidate solutions. At first, this seems like it would completely prohibit calling the function, but that isn’t the case. Recall our integration of the control strategy:
equation
when sample(0, 0.01) then
Q = computeHeat(T, Tbar, Qcapacity);
end when;
C*der(T) = Q-h*(T-Tamb);
end HysteresisEmbeddedControl;
In particular, note that computeHeat
is invoked only within the
when
statement and not as part of a “continuous” equation. As a
result, we can be certain that computeHeat
will only be invoked in
response to an event but not when evaluating candidate solutions for the
continuous variables.
In the C function computeHeat
, we see that these two statements
implement a +/- 2 degree deadband around the setpoint:
if (state==OFF && T<Tbar-2) state = ON;
if (state==ON && T>Tbar+2) state = OFF;
It is this functionality that prevents chattering and which can be clearly observed in the simulated results for our example: