Function 'odeToVectorField' introduces undesired 'free' variable
5 views (last 30 days)
Show older comments
Good day everyone,
I am having a problem with Matlab built in function odeToVectorField and i cannot find the solution on the available documentation.
I am tryng to derive using the Lagrangian approach the Equations of Motion (EoMs) of a multibody (rigid) system characterized by several Degrees of Freedom (DoFs). To do so I am using Matlab symbolic tool as reported, for instance in the example at: https://it.mathworks.com/help/symbolic/odetovectorfield.html.
My problem is that the command odeToVectorField introduces a free variable 'x3' which I did not define nor expect. Morover, I cannot understand the phisical meaning of it.
In order to make my problem clear I add below the code of a simplified problem, consisting in the derivation of the EoMs of a 3 DoFs system (a body in the 2D cartesian space) with only 1 DoF controlled by a torque.
% Code:
clear all; close all; clc
%% State Variables
% Position [m] and Orientation [rad] of the base
syms x(t) y(t) phi(t) dx(t) dy(t) dphi(t) ddx(t) ddy(t) ddphi(t)
% Base Torque [Nm]
syms TAU(t)
% Mass [kg]
syms m0
% Geometry of base
syms a b
% Inertia Moments in Principal Axes Frame {base-frame}
syms I0
% Lagrangian variables
syms L K V
%% Position and Velocity of the Masses in inertial frame {N}
% Position of Mass 0 (base)
x0 = x; y0 = y;
xx0 = [x0, y0]';
% Velocities [m/s]
vv0 = diff(xx0,t);
%% Moments of Inertia in {base-frame} [kg*m^2]
% Total Inertia
Itot = I0;
%% Lagrangian
% Kinetic Energy
K = 0.5*(m0*(vv0'*vv0) + Itot*dphi^2);
% Substitution of variables
K = subs(K, diff(x, t), dx);
K = subs(K, diff(y, t), dy);
K = subs(K, diff(phi, t), dphi);
% Potential Energy
V = 0;
% Lagrangian
L = K - V;
%% Partial Derivatives of the Lagrangian
% Partial wrt state
pL_x = simplify(diff(L, x));
pL_y = simplify(diff(L, y));
pL_phi = simplify(diff(L, phi));
% Partial wrt state (time) derivatives
pL_dx = simplify(diff(L, dx));
pL_dy = simplify(diff(L, dy));
pL_dphi = simplify(diff(L, dphi));
%% Generalized Equations of Motion:
EoM1 = diff(pL_dx, t) - pL_x == 0; % null ext. force
EoM2 = diff(pL_dy, t) - pL_y == 0; % null ext. force
EoM3 = diff(pL_dphi, t) - pL_phi == TAU; % actuated base
%% Substitution of Variables in EoMs:
% EoM1
EoM1 = subs(EoM1, dx, diff(x, t));
EoM1 = subs(EoM1, ddx, diff(x, t, t));
EoM1 = subs(EoM1, dy, diff(y, t));
EoM1 = subs(EoM1, ddy, diff(y, t, t));
EoM1 = subs(EoM1, dphi, diff(phi, t));
EoM1 = subs(EoM1, ddphi, diff(phi, t, t));
% EoM2
EoM2 = subs(EoM2, dx, diff(x, t));
EoM2 = subs(EoM2, ddx, diff(x, t, t));
EoM2 = subs(EoM2, dy, diff(y, t));
EoM2 = subs(EoM2, ddy, diff(y, t, t));
EoM2 = subs(EoM2, dphi, diff(phi, t));
EoM2 = subs(EoM2, ddphi, diff(phi, t, t));
% EoM3
EoM3 = subs(EoM3, dx, diff(x, t));
EoM3 = subs(EoM3, ddx, diff(x, t, t));
EoM3 = subs(EoM3, dy, diff(y, t));
EoM3 = subs(EoM3, ddy, diff(y, t, t));
EoM3 = subs(EoM3, dphi, diff(phi, t));
EoM3 = subs(EoM3, ddphi, diff(phi, t, t));
%% Substitution of geometric and mass values
% Lentghs [m]
a = 2;
b = 1;
% Masses[Kg]:
m0 = 100;
m1 = 2;
m2 = 1;
% Inertia of cuboid [Wikipedia: https://en.wikipedia.org/wiki/List_of_moments_of_inertia#Moments_of_inertia]
I0 = (1/12)*m0*4*(a^2*b^2); % 4* due to a,b being HALF of width and depth
% Substitution in EoMs
EoM1 = subs(EoM1);
EoM2 = subs(EoM2);
EoM3 = subs(EoM3);
%% Convertion to Vector-Field
[V,S] = odeToVectorField(EoM1, EoM2, EoM3);
%% Convertion to Matlab Function (GIVES PROBLEM)
% ode_sys = matlabFunction(V,'vars',{'t', 'Y'});
Thank you for time and help!
0 Comments
Accepted Answer
Paul
on 21 May 2022
By default, the Symbolic Toolbos assumes all variables are complex. So using the ' operator, which is the complex transpose operator, introduces some issues. Change these two lines to just use the tranpose operation
%xx0 = [x0, y0]';
xx0 = [x0, y0].';
% Kinetic Energy
%K = 0.5*(m0*(vv0'*vv0) + Itot*dphi^2);
K = 0.5*(m0*(vv0.'*vv0) + Itot*dphi^2);
Also, it's good to get in the habit of using assumptions as appropriate, e.g.,
syms m0 positive % assumes mass is positive
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!