Why does the eig function return different values in a Simulink 'MATLAB Function' block versus the MatLab command line?

Hello,

I have been having problem with a function I've written which relies on the zeros of a system. The problem has been narrowed down to the eig function that is being called. For the given system matrices, the computation in MatLab returns 4 finite values, corresponding to the zeros of the system. However, when packaged up into a Simulink 'MATLAB Function' block, the same code returns only 3 finite values, and they are quite different. Why do they return different answers? Are different algorithms being used for the computation? Even so, I would expect the same answer, up to numerical precision. I will paste the specific matrices I was using when I encountered this issue. It seemed to work in some other cases.

 Atemp =     [ 0    1.0000         0         0         0         0         0         0;
   -0.5291   -1.3801    0.3567   -0.3759    0.1828   -0.0428    0.0677   -0.3570;
         0         0         0         0         0    1.0000         0         0;
         0         0         0         0         0         0    1.0000         0;
         0         0         0         0         0         0         0    1.0000;
         0   -1.0000  -18.2230   -0.2500    0.2500  -20.4795         0         0;
         0   -1.0000    0.2500  -15.8050    0.2500         0  -17.4062         0;
         0    0.2918   -0.1782    0.1393   -9.6273         0         0  -11.5506];
Btemp =         [0         0         0         0;
    0.5291         0         0         0;
         0         0         0         0;
         0         0         0         0;
         0         0         0         0;
         0   18.4730         0         0;
         0         0   15.5585         0;
         0         0         0    9.5866];
Ctemp =      [1     0     0     0     0     0     0     0;
     0     0     1     0     0     0     0     0;
     0     0     0     1     0     0     0     0;
     0     0     0     0     1     0     0     0];
% coder.extrinsic('eig')
zcl = zeros(size(Atemp,1)+1,1)*1i;
zcl(1:size(Atemp,1)+1,1) = eig([Atemp Btemp(:,3); Ctemp(2,:) 0],diag([ones(1,size(Atemp,1)) zeros(1,1)]));
out = nnz(abs(zcl)<1e6);

Placing the above code into a MATLAB Function block in Simulink and commenting/uncommenting the coder.extrinsic command will cause the output to switch between 4 and 3. To be clear, none of the finite computed values exceeds 1e6, so it is not being lost there. The computed values are:

zcl =
    Inf + 0.0000i
   -Inf + 0.0000i
 0.0423 + 0.0000i
 0.7100 + 0.0000i
-0.8209 + 0.0000i
   -Inf + 0.0000i
   -Inf + 0.0000i
    Inf + 0.0000i
    Inf + 0.0000i

versus

      Inf + 0.0000i
      Inf + 0.0000i
 -10.6663 + 0.0000i
   0.0510 + 0.6434i
   0.0510 - 0.6434i
  -0.8891 + 0.0000i
      Inf + 0.0000i
      Inf + 0.0000i
      Inf + 0.0000i

Any help in understanding this would be appreciated. Thank you

Answer 1

[Andrew Schenk]

Since the MATLAB Function Block uses code generation, the same limitations and notes apply. For code generation, these notes found in the following documentation can provide some additional inside as to the limitations and differences: http://www.mathworks.com/help/simulink/ug/functions-supported-for-code-generation--alphabetical-list.html

When you use coder.extrinsic('eig'), code is not generated for the eig function, but instead Simulink directly calls into the MATLAB engine. As opposed to code generated from the eig function, MATLAB uses a number of different algorithms as well as numerical libraries to solve the Eigenvalue problem.

The problem you are trying to solve has a singular matrix (same as this example). Because of this singularity, the algorithm and numerical library differences between generated code and the MATLAB engine cause the difference in calculated Eigenvalues.

Using coder.extrinsic('eig') will ensure you always have the same answer in Simulink as well as in MATLAB. However, perhaps a better solution is to reformulate the problem you are trying to solve so you do not need to directly solve for Eigenvalues of a singular matrix.