# Solving Linear Systems for Multibody Systems

10 views (last 30 days)
Tiago Carvalho on 21 Jul 2022
Commented: Tiago Carvalho on 22 Jul 2022
Good afternoon,
I'm currently coding a Multibody Foward Dynamcis Simulator in MATLAB, for my master thesis, and I am having issues with solving linear systems due to low RCOND values, ill-condition of the matrices.
This problem is expected in this type of problems since each column entry in a line corresponds to the elements of a joint equation for two bodies resulting in a highly sparse matrix and consequently ill-conditioned matrices.
What my algorithm does is it calculates the initial accelerations of bodies in order for them to be integrated through an ode solver for position and velocity. My problems arise in solving the linear system for the initial accelerations that will be integrated. The system I have to solve is the following:
Where: M is a mass matrix, Phiq is a Jacobian, alpha/omega/mu are scalars, Phi, gamma and upsilon are vectors.
For my algorithm I tried using pinv(A)*B or lsqr(A,B), but I am not able to get the correct results (only get constant values from initial time until the last integration). On the other hand I get good results by using mldivide or \, but i get the following message:
Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 1.394430e-19.
I wanted to know if there is anyway to avoid this issue, since I know it can cause significant noise in my results, I usually use pinv for this but it seems to not be working this time.
Thank you for your time and attention.
Tiago
(Attached files are the left and right hand side of the linear problem)
Torsten on 21 Jul 2022
Edited: Torsten on 21 Jul 2022
MATLAB's ODE solvers are designed to solve systems of the form
M(t,y)*y' = f(t,y)
Why do you invert your matrix in the ODE function routine and supply M^(-1)*f and don't let the ODE solver do the job by defining two functions in which you separately define the mass matrix M(t,y) and the right-hand side vector f(t,y) ?
Tiago Carvalho on 21 Jul 2022
Edited: Tiago Carvalho on 21 Jul 2022
Hello Torsten,
I don't think I can do that in this case in particular, since my algorithm follows the flowchart in the Picture. The process of this function it self is iterative inside the ode, in order to ensure that the constraints of my specific problem are solved, I don't see how I could integrate that outside the ode to allow him to do the inversion.
(EDIT) What the ode integrates is the product of the equation that I previously shared and the velocity that is an input of the odefun.

Bjorn Gustavsson on 22 Jul 2022
As best I can interpret your flow-chart it seems that it ought to be "reasonably straightforward" to follow Torsten's advice by converting the equation for into one for both and . If I get it right it should be something like:
function Mout = modified_massmatrix(M)
Mout = [eye(size(M,1)),zeros(size(M));
zeros(size(M)),M];
end
or if the mass-matrix is a function of t, q and :
function Mout = modified_massmatrix(t,qqdot,M)
Mnow = M(t,qqdot);
Mout = [eye(size(Mnow,1)),zeros(size(Mnow));
zeros(size(Mnow)),Mnow];
end
Then you should be good to go with the ODE-integrating functions.
Tiago Carvalho on 22 Jul 2022
Yes my issue is very similar with what you have described! What I am trying to solve is the equations of motion for constrained multibody systems in relation to the bodies accelerations. Where i use the Lagragian method:
However, these equations of motion are only dependent on the Lagrangian Mutipliers and in the Acceleration vector, and do not explicitly use the position equations, written using kinematic constraints. This leads to an accumulation of a quadratic error in the position solution (the integration of velocity qd).
The only way I can avoid this is if I implement a "stabilization constraints algorithm" (like a feedback loop in control theory), in this case I use the Augmented Lagrange Formula, a penalty method, to stabilize the constraints in each iteration ( the while loop).
Nevertheless, I will try to implement th suggestions that you gave me to see if I can solve my issues. I have already follow your advise in simplifying the problem and I have reduce it to a basic one (a body which actuated only by gravity) and in this case I don't have any issues, so I think maybe the problem comes when I introduce other variables to the problem (like springs and dampers, which are part of the g vector).
Anyway, thank you for your help. I understand that probably I haven't been 100% explicit in my questions, but I think you can understand that my problem is a bit complex and is difficult to me to explain these in simpler terms when I have been working on this code for months already. For this I want to apologize.
Thank you again for your help, it is much appreciated.
Tiago Carvalho
Tiago Carvalho on 22 Jul 2022
Hello @Sam Chak, no I use the standard Newton-Euler's equation and Lagrange Method, not Kane's.
My main source of information was the P.E Nikravesh book titles Computer Aided Analysis of Mechanical Systems, which suggested the use of the Newton-Euler's, to be 100% honest I am not familiriazed with Kane's method but I will research it to see if it can simplify my issues!