I am using ode15s to resolve a DAE system. Since there are 300 variables x1(t)~x300(t) in my case, so I feel it is a little bit inconvenient to define them by writing
syms x1(t) … x300(t)
one by one. Are there any efficient codes to define them?
qd_t( i ) = str2sym( char( ( 'qd_t' )' , num2str( i )' , ('(t)')')' ) ;
end
which can define a column vector that can be applied in the next lines in my file. So are there any more compact or concise methods to define a column vector?
Thanks so much for your nice answer and suggestion!
Frankly, I am trying to solve a very typical problem, i.e., the dynamics of a constrained mechanical system. I can format the model as
where q is a the genialized coordinate vector, andare its first and second time derivatives, respectively, and is the Lagrange multiplier vector.
From many publications, I know it is an index-1 DAE system and many researchers use ode15s to numerically resolve this equation. However, I don’t know how to format it into a suitable form that can be used by ode15s, such as the mass matrix, the f vector on the right side of the equation, and the state variables. Frankly, I have the numerical values of q ,,, and at the time t=0. Could you pls give me any hints to rewrite this equation into a form suitable for ode15s?
The fact that there are specialized codes for Mechanical Systems indicates that rewriting the equations and solving them with a standard solver like "ode15s" won't be that easy.
As said: I strongly recommend using a code that is especially designed for the solution of mechanical systems. I only know the FORTRAN code by the well-known experts Hairer/Wanner I linked to, but maybe there are other free codes available.
In the odeset of my DAE problem using ode15s, I supplied the Jacobian as a function handle in odeset. However, the statistics show that the number of times that the partial derivative matrix is not zero. What is the reason for it? I think with the supplied Jacobian matrix, this number should be zero.
PHI is the 18 * 1 algebraic constraint equations, and par_mtx_2_rank_vtr is a small function I wrote by myself to compute the derivative of a symbolic matrix with respect a rank vector.
Thus, PHI is the function of the symbolic vector q_d,
and PHI_q_dot is the functio of both the symbolic vector q_d and q_v.
Hi Torsten, please see the attached file. Pls firstly run I_symbolic_FD.m to get the symbolic results and then I_numeric_FD.m to get the numerical results.
cnt: 1
23 successful steps
1 failed attempts
170 function evaluations
2 partial derivatives
8 LU decompositions
32 solutions of linear systems
Elapsed time is 2.309397 seconds.
cnt: 2
20 successful steps
0 failed attempts
96 function evaluations
1 partial derivatives
6 LU decompositions
26 solutions of linear systems
Elapsed time is 0.144984 seconds.
cnt: 3
19 successful steps
0 failed attempts
95 function evaluations
1 partial derivatives
6 LU decompositions
25 solutions of linear systems
Elapsed time is 0.115081 seconds.
cnt: 4
18 successful steps
0 failed attempts
94 function evaluations
1 partial derivatives
6 LU decompositions
24 solutions of linear systems
Elapsed time is 0.118874 seconds.
cnt: 5
19 successful steps
1 failed attempts
96 function evaluations
1 partial derivatives
7 LU decompositions
27 solutions of linear systems
Elapsed time is 0.139268 seconds.
cnt: 6
22 successful steps
1 failed attempts
99 function evaluations
1 partial derivatives
7 LU decompositions
30 solutions of linear systems
Elapsed time is 0.103593 seconds.
cnt: 7
20 successful steps
1 failed attempts
100 function evaluations
1 partial derivatives
7 LU decompositions
28 solutions of linear systems
Elapsed time is 0.094056 seconds.
cnt: 8
20 successful steps
0 failed attempts
98 function evaluations
1 partial derivatives
6 LU decompositions
26 solutions of linear systems
Elapsed time is 0.083510 seconds.
cnt: 9
19 successful steps
0 failed attempts
97 function evaluations
1 partial derivatives
6 LU decompositions
25 solutions of linear systems
Elapsed time is 0.083450 seconds.
cnt: 10
23 successful steps
0 failed attempts
101 function evaluations
1 partial derivatives
6 LU decompositions
29 solutions of linear systems
Elapsed time is 0.085572 seconds.
cnt: 11
46 successful steps
12 failed attempts
719 function evaluations
9 partial derivatives
21 LU decompositions
111 solutions of linear systems
Elapsed time is 0.559771 seconds.
cnt: 12
44 successful steps
19 failed attempts
1467 function evaluations
20 partial derivatives
29 LU decompositions
122 solutions of linear systems
Elapsed time is 0.890452 seconds.
cnt: 13
47 successful steps
13 failed attempts
1060 function evaluations
14 partial derivatives
23 LU decompositions
118 solutions of linear systems
Elapsed time is 0.562474 seconds.
cnt: 14
26 successful steps
9 failed attempts
814 function evaluations
11 partial derivatives
16 LU decompositions
72 solutions of linear systems
Elapsed time is 0.359789 seconds.
cnt: 15
44 successful steps
17 failed attempts
1266 function evaluations
17 partial derivatives
26 LU decompositions
122 solutions of linear systems
Elapsed time is 0.438848 seconds.
cnt: 16
30 successful steps
3 failed attempts
405 function evaluations
5 partial derivatives
11 LU decompositions
66 solutions of linear systems
Elapsed time is 0.116305 seconds.
cnt: 17
23 successful steps
1 failed attempts
176 function evaluations
2 partial derivatives
8 LU decompositions
37 solutions of linear systems
Elapsed time is 0.049336 seconds.
cnt: 18
21 successful steps
0 failed attempts
103 function evaluations
1 partial derivatives
6 LU decompositions
31 solutions of linear systems
Elapsed time is 0.030169 seconds.
Elapsed time is 0.030273 seconds.
function [mass_matrix_handle,F_f_handle,Jacobian_handle] = symbolic_FD
yeah, allthough we specify the Jacobian matrix handle, the partial derivatives are still needed in each count.
No. You can supply them, but if not, ode15s will compute them by finite difference approximations. So you can delete the computation of the Jacobian matrix handle.
Why do you think the partial derivatives are needed ?
How did you arrive at the results you compare with (q_dv_lambda_val) ?
The reason to specify the Jacobian matrix is, our DAEs are very stiff: some components of y are pretty tiny while some others are extremely huge. From the help center of ode15s, it is critically important to provide the Jacobian.
And do you have an explanation for the differences between the q_dv_lambda_val and the q_dv_lambda_FD results ? Maybe that you keep torque_LA constant in each subinterval of length t_h ? Did you try whether the integration still works for even smaller values of RelTol and AbsTol ?
Out of interest, I solved the pendulum problem in cartesian coordinates directly according to the equations you listed above:
Maybe it's worth trying to solve your equations this way because you don't need to put lambda as a solution variable and you don't need to work with a mass matrix.
% Solves the pendulum equation in cartesian coordinates
Instead of solving the DAE, have you considered alternatives, such as Kane's Method, to instead formulate the problem as an ODE that can be integrated with an ode solver?
one only needs an ODE integrator to solve for position (q) and velocity (q_dot).
This is accomplished by solving the above linear system of equations for q_doubledot (and lambda, but this is not important in the sequel) using backslash (\) and then setting up the ode system as
dq/dt = q_dot
dq_dot/dt = q_doubledot
I doubt there is any method that does not need to integrate for q and q_dot in some way.
I agree that the approach you've described might be viable for the specific problem under discussion.
Assuming I understand the setup of the problem under discussion .... two things to keep in mind.
1. The size of the matrix to be inverted increases as the number of constraints in the system increases. The dimension of q in the problem is not clear to me, though the original question suggests it might be on the order of hundreds (6 coordinates for each of 50 elements in the system?). So we might be starting with M quite large to begin with, and then the size of the matrix to be inverted grows with number of constraints (i.e., the dimension of lamda). With Kane's Method the size of the matrix to be inverted decreases as the number of constraints increases, which results in computational savings (which might or might not be significant).
2. The suggested approach will inherently be influenced by computational innaccuracy, which can result in the integrated solution drifting and not satisfying the constraint equation Phi*qdot = 0 (or its integral if the constraints apply at the configuration level). I know that you're aware of this issue as shown in your "constraint error" plot, but I thought it worth mentioning again. With Kane's method, the constraints are inherently incoporated into the EOM for the minimal* number of coordinates and therefore "constraint drift" is not a concern.
*I believe that there are formulations of Kane's Method for a non-minimal number of coordinates, but I'm not familiar with that approach nor its motivations.
I cannot recommend a certain method because I have no experience with the constrained mechanical system equations. That's why I gave a link to a FORTRAN package that solves this special type of DAE system. Since the authors of the package are experts in the solution of differential-algebraic problems, I guess they took care of the drift problem you described above.
But as @Tony Cheng already prepared everything in MATLAB, I feel it's worth testing the method I used in the solution of the pendulum equations for his problem because I guess that only little changes to his code will be necessary. Further, comparison between the results from his formulation and the new one could be of interest.
Looking again at @Tony Cheng 's code, I'm surprised about his right-hand side of the DAE-system. In particular, I cannot make sense of the line
PHI_q * q_v_dot = - ( PHI_q_dot + 2 * ALPHA * PHI_q ) * q_v - BETA ^ 2 * PHI , because I use the Baumgarte stablization method (well-known in constrained mechanical systems) to try the best to reduce the position, velocity, and acceleration drift errors.
I gather that the DAE is shown in equation (31), though the upper-right term on the LHS should be transpose(phi_q).
From what I can tell after a quick skim, the authors intentionally formulate the problem to yield a DAE to avoid "sophisticated velocities transformation between the joint-space and task-space." But that was a choice, not a requirement of the method.
Though the authors are solving a DAE, they explicitly say that they use Matlab ode45. Perhaps they are using the method suggesed by @Torsten.
The coefficient matrix on the LHS in this comment is singular. Is that expected?
you choose q_dv_lambda_val as initial condition. Shouldn't this be q_dv_lambda_FD because you want to continue an integration with solution variable q_dv_lambda_FD ? Otherwise, the integrations in the loop are all independent from another since q_dv_lambda_val is prescribed externally - thus there is no temporal continuation of a process.
At the instant cnt*delta_h, the actuating torques fed into the generalized force vector are updated.
That's right, but you must continue the computation with the solution so far obtained with ode15s, not with a solution obtained by an external software (or whereever the q_dv_lambda_val array stems from).
Yes, u are rite. Now I see the meaning of long time simulations. if q_dv_lambda_val is used as initial condition for each loop, then the integrations in one loop are all independent from another, having nothing to do with the long time simulations.
However, if the q_dv_lambda_FD obtained from the end of the last loop is used as the initial condition of the next loop, the violation in terms of displacement, velocity, and acceleration is very huge in the entire time span.
However, if the q_dv_lambda_FD obtained from the end of the last loop is used as the initial condition of the next loop, the violation in terms of displacement, velocity, and acceleration is very huge in the entire time span.
I noticed this, too. So either the models used to computeq_dv_lambda_FD and q_dv_lambda_val are not comparable or you or the one responsible for the q_dv_lambda_val results made a mistake in the implementation.
That could be. From the articles I find that the index-1 DAEs have no acceleration violations. Violations only exist in coordinates and velocity levels. However, from my codes, the acceleration violation is the largest among the three.
Probably I will turn to Simscape multibody for some simulations @Paul.
These days I am reading relevant papers and find that in the form
the largest contraint violations occur in displacement and velocity levels, while the violations in acceleration are the smallest. I agree with what you suspect.
Displacement and velocity are the solution variables, acceleration is a deduced quantity that is computed from these solution variables. Thus violations in acceleration only follow from violations in displacement and velocity (and force terms) - they should not occur per se.
This can easily be seen if you look at the underlying system:
Acceleration and Lagrage parameters can directly be computed from displacement and velocity (and force terms) by solving a linear system of equations. Thus they are only "wrong" if q, qdot, F, C and/or G are "wrong".
Using the above linear system to compute acceleration and Lagrange parameters from displacement, velocity and forces, I suggest computing acceleration and Lagrange parameters with your code using the displacements and velocities saved in "q_dv_lambda_val". They should be similar to the accelerations and Lagrange parameters of the other code if you expect similar results for the two simulations.