Function 'odeToVectorField' introduces undesired 'free' variable

조회 수: 2(최근 30일)
Alessandro Scalvini
Alessandro Scalvini 2022년 5월 21일
댓글: Alessandro Scalvini 2022년 5월 21일
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:
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:]
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!

채택된 답변

Paul 2022년 5월 21일
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
  댓글 수: 1
Alessandro Scalvini
Alessandro Scalvini 2022년 5월 21일
Thank you very much for the answer and the useful tip!

댓글을 달려면 로그인하십시오.

추가 답변(0개)




Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by