Coriolis Matrix of a 14 DOF bipedal robot won't compile its matlab function

Hi, I have this 14 DOF bipedal robot, programmed on matlab. The trajectory works (which is for now to elevate the legs and put them back), the Inverse Kinematics works (I have already tested it).

Now I'm trying to solve the dynamic equations to obtain all the torques I need.

The dynamic script provides me the Mass matrix, the coriolis matrix and the gravitational force matrix. I should put them all togheter in this way M(q)q¨+C(q,q˙)q˙+G(q)=τ(+J⊤λ) into a matlab function on simulink and get the torques.

Unfortunately my script makes all the matlab functions with the exeption of the Coriolis matrix (CC).

I'm really trying everything to smooth the code, but now I'm starting to think that there is some error that I cannot recognize, that is not only slowing the CC_MORM down, but maybe is blocking it. I can't tell.

% Second derivative of the generalized coordinates
syms ddN ddE ddD dde0 dde1 dde2 dde3 real 
ddeta = [ddN; ddE; ddD; dde0; dde1; dde2; dde3; ddqtot];

Unrecognized function or variable 'ddqtot'.

syms NN(t) EE(t) DD(t) ee0(t) ee1(t) ee2(t) ee3(t) q0(t) q1(t) q2(t) q3(t) q0s(t) q1s(t) q2s(t) q3s(t) t
eeta = [NN(t) EE(t) DD(t) ee0(t) ee1(t) ee2(t) ee3(t) q0(t) q1(t) q2(t) q3(t) q0s(t) q1s(t) q2s(t) q3s(t)].';
deeta = diff(eeta,t);
ddeeta = diff(deeta,t);

% Substitute symbolic variables with their time-dependent counterparts.
% Perform the substitutions in this order to obtain the correct result
LL = subs(L, eta, eeta);
LL = subs(LL, deta, deeta);
LL = subs(LL, ddeta, ddeeta);

% Lagrangian derivative w.r.t. eta: dL/deta
dLL_deeta = jacobian(LL,eeta).';
% Lagrangian derivative w.r.t. deta: dL/ddeta
dLL_ddeeta = jacobian(LL,deeta).';
% Euler–Lagrange dynamic equations set
DynEq_eeta = diff(dLL_ddeeta) - dLL_deeta - tau;
phi_q = qT.' * qT - 1;          % constraint (no time derivatives)
EulerLagrangeEq = [DynEq_eeta; phi_q];

% Back-substitute to the original symbolic variables
EulerLagrangeEq = subs(EulerLagrangeEq, ddeeta, ddeta);
EulerLagrangeEq = subs(EulerLagrangeEq, deeta, deta);
Euler_Lagrange_Eq = subs(EulerLagrangeEq, eeta, eta);

etastar = [eta; lambda];
detastar = [deta; lambda];
ddetastar = [ddeta; lambda];
M = jacobian(Euler_Lagrange_Eq, ddetastar);

%% CALCULATIONS OF DERIVATIVES FOR CORIOLIS MATRIX ELEMENTS
% Compute the Coriolis matrix C
% Initialization
n = length(detastar);
dM = sym(zeros(n,n,n));
CC = sym(zeros(n,n));
% --- PROGRESS BAR ---
waitbarHandle = waitbar(0, 'Symbolic derivative of matrix M...');
total = n^3;
counter = 1;
% Symbolic derivative calculation
for i = 1:n
    for j = 1:n 
        for k = 1:n
            dM(i,j,k) = diff(M(i,j), etastar(k));
            waitbar(counter / total, waitbarHandle, ...
                sprintf('Derivative progress: %3.0f%%', 100 * counter / total));
            counter = counter + 1;
        end
    end
end
close(waitbarHandle);

%% COMPUTE CORIOLIS COMPONENTS 
waitbarHandle = waitbar(0, 'Computing symbolic Coriolis components...');
total = n^2;
counter = 1;
for i = 1:n
    for j = 1:n
        CC_ij = sym(0);  % initialize the j-th element of the Coriolis matrix
        for k = 1:n
            C_ijk = 0.5 * (dM(i,j,k) + dM(i,k,j) - dM(j,k,i)) * detastar(k);
            CC_ij = CC_ij + C_ijk;               
        end
        CC(i,j) = CC_ij;
        waitbar(counter / total, waitbarHandle, ...
            sprintf('Coriolis components: %3.0f%% completed (%d/%d)', 100 * counter / total, counter, total));
        counter = counter + 1;
    end
end
close(waitbarHandle);

%% CREATE GRAVITATIONAL FORCE VECTOR 
G = sym(zeros(16, 1));
G = jacobian(-U, etastar);

%% CREATE JACOBIAN FOR QUATERNIONS
phi = qT.' * qT - 1;
J_phi_q = jacobian(phi, qT);   % 1x4

%% CREATE MATLAB FUNCTIONS (fast export)
vars_pack = {eta, deta, Masse, Ib, I0, I1, I2, I3, L0, L1, L2, L3, L_offs, CG, g};

if exist('M_MORM.m','file')==2
    fprintf('M_MORM already exists. Skip.\n');
else
    matlabFunction(M,  'Vars', vars_pack, 'File', 'M_MORM',  'Optimize', true);
    fprintf('M_MORM exported.\n');
end

if exist('CC_MORM.m','file')==2
    fprintf('CC_MORM already exists. Skip.\n');
else
    matlabFunction(CC, 'Vars', vars_pack, 'File', 'CC_MORM', 'Optimize', false);
    fprintf('CC_MORM exported.\n');
end

if exist('G_MORM.m','file')==2
    fprintf('G_MORM already exists. Skip.\n');
else
    matlabFunction(G,  'Vars', vars_pack, 'File', 'G_MORM',  'Optimize', true);
    fprintf('G_MORM exported.\n');
end

% J_phi 1x16 with respect to etastar
phi   = qT.' * qT - 1;
J_phi = jacobian(phi, etastar);

if exist('J_phi_func.m','file')==2
    fprintf('J_phi_func already exists. Skip.\n');
else
    matlabFunction(J_phi, 'Vars', {etastar}, 'File', 'J_phi_func', 'Optimize', true);
    fprintf('J_phi_func exported.\n');
end   

[Torsten]

ddqtot, L, ... are undefined.

Answer 1

[Umar]

Hi @Giovanni Maria,

After spending several hours researching MATLAB documentation and robotics literature online, I believe I've identified the root cause of your CC_MORM export failure and can offer some practical solutions that don't require additional toolboxes.

Your issue isn't a coding error - it's computational complexity. For your 14-DOF system with quaternion constraints, the symbolic Coriolis computation involves 4,352 operations on complex trigonometric expressions, causing memory bottlenecks that prevent matlabFunction from completing the export. So, my recommendations would be:

Option 1 - Enable Optimization: Change your Coriolis export line to:

matlabFunction(CC, 'Vars', vars_pack, 'File', 'CC_MORM', 'Optimize', true);

You currently have 'Optimize', false which makes the expressions much larger.

Option 2 - Simplify Before Export: Add simplification before creating the function:

CC_simplified = simplify(CC, 'Steps', 100);
matlabFunction(CC_simplified, 'Vars', vars_pack, 'File', 'CC_MORM', 
'Optimize', true);

Option 3 - Memory Management: Add between your derivative calculations:

clear dM;
pack;

Option 4 - Numerical Coriolis at Runtime: Keep your symbolic M and G (which already work), but compute Coriolis numerically in your Simulink function block for the full dynamic equation M(q)q_ddot + C(q,q_dot)q_dot + G(q) = tau.

Now, as a quick test, temporarily reduce your system to 8-10 DOF to verify the approach works, then implement one of the above optimizations for the full system.

Based on my research, this complexity issue is well-documented for multi-DOF systems. The symbolic expressions become exponentially complex, which is exactly what you're experiencing.

References from my research:

• MATLAB Symbolic Math Documentation: https://www.mathworks.com/help/symbolic/

• Robotics dynamics papers consistently recommend numerical methods for systems >10 DOF

[giovanni maria]

thank you for your time and effort. I've tried simplifying and clearing the dM but it's still to slow and it saturates the ram. I will proceed in two ways.

First I'm going to compute CC numerically as you suggested.

Then if it this won't work I'll try to use a computer with a 32 or 64GB ram if I find one xD, since my pc have only 16.

[Umar]

Hi @Giovanni Maria,

For a 14-DOF robot, calculating the Coriolis matrix symbolically is too heavy for your PC and will always be slow or fail.

The best approach is to:

1. Compute the Coriolis matrix numerically during runtime instead of generating a symbolic function. 2. Keep the mass and gravity matrices symbolic—they work fine. 3. Optionally, test on a smaller system first to make sure your method works.

More RAM might help, but numerical computation is the practical solution.