Verifying Simulink generated C code

I generated embedded C code for the plant in the following MATLAB Simulink model:

To verify that the generated code behaves the same as the Simulink model, I could save the controller outputs to a file, create a C project with a main() function to read those outputs, and feed them to the plant C code at each time step—essentially converting a closed-loop simulation into an open-loop one. This approach would work if I use a single-step numerical integrator/solver like ode1 (Euler method). However, I usually prefer ode4 (Runge-Kutta 4) due to its better accuracy with less computation. Despite RK4 having 4 intermediate computations, it has a global error of O(h^4) and Euler O(h). If RK4 uses a step size of h, Euler would need h^4 to achieve similar accuracy. So, in total Euler would require more function evaluations than RK4—specifically, (1/h^3)/4 times more. For example, if hRK4 = 0.1, Euler would need 250x more computations.

Unfortunately, RK4 is a multi-step method with minor step sizes. Since the Simulink model saves controller outputs only at major steps, reading them from a file to drive the plant might result in divergent behavior—especially if the system has high-frequency dynamics—because controller outputs could differ at minor steps, while reading from the file only accounts for major step outputs.

You can only fully verify the system by integrating the C code into your embedded setup, which generates controller outputs based on real time plant-sensor data. Since testing on an embedded system takes time, to catch issues such as uninitialized variables in the plant code, you can create a C project, build and run your plant model with constant controller inputs. If it compiles and runs without NaN values, you can proceed to the embedded testing phase.

If you want to be more rigorous, you could save the plant outputs together with the controller outputs. In your C project, you read the controller output, feed it to the plant, run it for only one time step, read the Simulink plant output for that time instance from a file, and compare it with the C code plant output. Even this comparison is not perfect and only works for the first time step because Simulink integrator block outputs depend on the system state in the previous time step. To check the whole time span, you also have to save the system state and initialize the system with the correct state at each time step. But this in effect would be equal to running the system for one time step, because you initialize the system, apply a constant control input and check the output.

Denormalized floating-point numbers

Recently, when debugging with Visual Studio 2019, I noticed that a variable had the unusual value of 6.324...e-322#DEN. This is "almost" equal to zero. In Visual Studio's debugger, the #DEN label refers to a denormalized (or subnormal) floating-point number. This occurs when the value of a floating-point variable is so close to zero that it cannot be represented in the normalized format of the IEEE 754 standard.

A floating-point number is typically represented as:

This is called a normalized number because the leading digit of the fraction (mantissa) is assumed to be 1.

When the exponent is at its smallest possible value (the minimum exponent), and the number is still too small to be represented in the normalized form, it is stored as a denormalized number. For denormalized numbers, the leading digit is 0 instead of 1, so the representation becomes:

Denormalized numbers fill the gap between the smallest normalized number and zero. 

By the way, it turned out that the strange 6.3...e-322#DEN value I was observing during debugging was due to an unitialized variable.

You can use the following C++ code to investigate further:

Avoiding catastrophic cancellation

Catastrophic cancellation occurs when subtracting two nearly equal numbers, leading to a significant loss of precision due to limited precision of floating-point numbers in a computer. Below are some examples and their better versions (reference):

log(a)−log(b) = log(a/b)

Computing cos(x)−1 for small values of x --> cos(x) ≈ 1−(x^2)/2 

Applying derivative and integral on noisy data

Let's say you have a motion with sinusoidal acceleration, i.e. xddot(t) = sin(t). The analytic solution for speed and position will be xdot(t) = -cos(t)+1 (assuming zero initial speed) and x(t) = -sin(t) + t.  Let's say you don't know the analytical expressions of position, speed and acceleration and only have their values at a fixed sample rate. If you want to increase the samples, you can use Matlab's spline function.

If you have position as input and calculate speed and acceleration, if there is any noise in the signal, it will be amplified. If you have acceleration as input and calculate speed and position via integration, noise will be attenuated. Another advantage of using integration instead of derivation is that we can use linear interpolation since when integrating, we don't need smooth second derivatives. This both decreases computational load and lag due to two points in linear interpolation compared to four points in spline. See below Matlab code for an example: