- #1
Mishal0488
- 16
- 0
Hi Guys
I am looking for some guidance with regards to a Lagrangian problem I am trying to solve.
Please refer to the attached documents.
Please neglect all the forcing functions for the time being. I am currently just trying to simulate the problem using initial conditions only
I have written Matlab code to simulate the problem, however my results are not as expected. the "theta2" variable always has a 0 solution, in some cases imaginary numbers form part of the solution. using initial conditions on the theta2 variable generally causes no solution as the ODE solvers on Matlab develop tolerance errors.
Note that the holonomic constraint whereby (b+r) is the subject of the formula is developed on Matlab using the symbolic package.
I believe the source of the issues is due to the holonomic constraint.
Can anyone provide me with some guidance on how to handle the system? and/or any mistakes which I have made. Could the system be represented in a more elegant manner?
Apologies for any funny terminology or poor representation of aspects, I am a structural engineer, Lagrangian mechanics isn't something that forms part of our undergraduate education hence this has been self taught.
I am looking for some guidance with regards to a Lagrangian problem I am trying to solve.
Please refer to the attached documents.
Please neglect all the forcing functions for the time being. I am currently just trying to simulate the problem using initial conditions only
I have written Matlab code to simulate the problem, however my results are not as expected. the "theta2" variable always has a 0 solution, in some cases imaginary numbers form part of the solution. using initial conditions on the theta2 variable generally causes no solution as the ODE solvers on Matlab develop tolerance errors.
Note that the holonomic constraint whereby (b+r) is the subject of the formula is developed on Matlab using the symbolic package.
I believe the source of the issues is due to the holonomic constraint.
Can anyone provide me with some guidance on how to handle the system? and/or any mistakes which I have made. Could the system be represented in a more elegant manner?
Apologies for any funny terminology or poor representation of aspects, I am a structural engineer, Lagrangian mechanics isn't something that forms part of our undergraduate education hence this has been self taught.
Matlab:
function xp = Statespace3(t,x,r,L,Cs,Cp,Fv1,Fv2,Fh,w,w2,w3,k,m1,m2)
A = (x(5)+r);
three = A*cos(x(1))*cos(x(3));
four = A*sin(x(1))*sin(x(3));
R1 = L^2*sin(x(3))^2;
R2 = -A^2*cos(x(1))^2*sin(x(3))^2;
R3 = -A^2*cos(x(3))^2*sin(x(1))^2;
R4 = -2*A*L*cos(x(3))^2*sin(x(1));
R5 = 2*A^2*cos(x(1))*cos(x(3))*sin(x(1))*sin(x(3));
R6 = 2*A*L*cos(x(1))*cos(x(3))*sin(x(3));
one = sqrt(R1+R2+R3+R4+R5+R6);
B = L*sin(x(3))-one+three+four;
three_dot = x(6)*cos(x(1))*cos(x(3))-A*x(2)*sin(x(1))*cos(x(3))-A*x(4)*cos(x(1))*sin(x(3));
four_dot = x(6)*sin(x(1))*sin(x(3))+A*x(2)*cos(x(1))*sin(x(3))+A*x(4)*sin(x(1))*cos(x(3));
R1_dot = 2*x(4)*sin(x(3))*cos(x(3))*L^2;
R2_dot = -2*x(6)*A*cos(x(1))^2*sin(x(3))^2+2*x(2)*A^2*cos(x(1))*sin(x(1))*sin(x(3))^2-2*x(4)*A^2*cos(x(1))^2*sin(x(3))*cos(x(3));
R3_dot = -2*x(6)*A*cos(x(3))^2*sin(x(1))^2+2*A^2*x(4)*cos(x(3))*sin(x(3))*sin(x(1))^2-2*A^2*cos(x(3))^2*x(2)*sin(x(1))*cos(x(1));
R4_dot = -2*x(6)*cos(x(3))^2*sin(x(1))^2+2*x(4)*A*cos(x(3))*sin(x(3))*sin(x(1))^2-2*x(2)*A*cos(x(3))^2*sin(x(1))*cos(x(1));
R5_dot = 4*x(6)*A*cos(x(1))*cos(x(3))*sin(x(1))*sin(x(3))-2*A^2*x(2)*sin(x(1))^2*cos(x(3))*sin(x(3))-2*A^2*cos(x(1))*x(4)*sin(x(3))^2+sin(x(1))+2*A^2*cos(x(1))*x(2)*cos(x(1))^2*sin(x(3))+2*A^2*cos(x(1))*sin(x(1))*x(4)*cos(x(3))^2;
R6_dot = 2*x(6)*L*cos(x(1))*cos(x(3))*sin(x(3))-2*A*L*x(2)*sin(x(1))*cos(x(3))*sin(x(3))-2*A*L*x(4)*cos(x(1))*sin(x(3))^2+2*A*L*x(4)*cos(x(1))*cos(x(3))^2;
X_dot = (R1_dot+R2_dot+R3_dot+R4_dot+R5_dot+R6_dot);
b_dot = L*x(4)*cos(x(3))+three_dot+four_dot-X_dot/(2*(one));
%theta1 coordinate
F_t1 = A*Fh*sin(w*t)*cos(x(1));
theta1_dotdot = (F_t1-2*m1*x(2)*x(6)*A-9.81*m1*sin(x(1))*A-Cp*x(2))/(m1*A^2);
%theta2 coordinate
F_t2 = B*Fh*sin(w*t)*cos(x(3));
theta2_dotdot = (F_t2-2*m2*x(4)*b_dot*B-9.81*m2*sin(x(3))*B-Cp*x(4))/(m2*B^2);
%a coordinate
Fa = Fh*sin(w*t)*sin(x(1))+Fv1*sin(w2*t);
a_dotdot = (Fa+m1*x(2)^2*A-k*x(5)+9.81*m1*cos(x(1))-Cs*x(6))/m1;
%b coordinate
b_dotdot = (m2*x(4)^2*B+k*(B-r)*9.81*m2*cos(x(3)))/m2;
%state vector
xp = [x(2);
theta1_dotdot;
x(4);
theta2_dotdot
x(5)
a_dotdot;];