Function 'odeToVectorField' introduces undesired 'free' variable

2 次查看（过去 30 天）
N/A
N/A 2022-5-21
评论： N/A 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: 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!

请先登录，再回答此问题。

采纳的回答

Paul
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Electrical Block Libraries 的更多信息

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by