How do I solve forward dynamics 2nd order ODE with time-dependent parameters using ODE45?
7 views (last 30 days)
Show older comments
I'm trying to solve a forward dynamics problem using ODE45, but the problem is tau is a vector in part (c) and so ODE45 is showing error when I'm using tau vector,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916234/image.jpeg)
I'm not able to solve (c) part (denoted by red dot) of my assignment. here is my code,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916239/image.jpeg)
This is the function I'm using, here tau = Torque, I've already calculated tau in first practical and I've to use it in this practical but I don't know how to do it,
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916244/image.jpeg)
5 Comments
Sam Chak
on 6 Mar 2022
Edited: Sam Chak
on 6 Mar 2022
Hi @Akash Vyas
The symbol tau (τ) is the control torque.
You've got the error because you put tau into an algebraic loop error that ODE45 cannot effectively solve:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916674/image.png)
Either it is a feedback control law (a formula) to compute the desired torque value that depends on θ and
, or a fixed torque value (or a series of torque values at certain time intervals where you need to inject them in to the system) that you obtained from Practical 1.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916614/image.png)
If it is a feedback control law, then tau looks like this:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916619/image.png)
The closed-loop system then becomes
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916624/image.png)
which can be simplified to
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916629/image.png)
So, the desired response of the 2nd-order system is determined by the selection of
the fundamental proportional gain,
and the fundamental derivative gain,
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916634/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916639/image.png)
Since the operation time only lasts for 3 seconds, I assume that the robot arm shall achieve the steady-state (
) at 1 second from the initial time. If that is acceptable to you, then you can try for
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916679/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916644/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916649/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916654/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916659/image.png)
Don't forget to multiply the fundamental control gains with the inertia
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/916669/image.png)
kp = 36;
kd = 12;
tau = - Izz*kp*theta - Izz*kd*dtheta + m*g*Lc*cos(theta);
Accepted Answer
Sam Chak
on 7 Mar 2022
Edited: Sam Chak
on 7 Mar 2022
With the new info on the desired joint angle trajectory profile
, the following simulations are based on 3 approaches on the control torque, τ.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917269/image.png)
Approach #1: Based on the computed torque from your Practical #1 when the desired
,
,
are known.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917269/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917279/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917284/image.png)
Approach #2: Based on the computed torque from Approach #1, but
is replaced by
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917289/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917294/image.png)
Approach #3: Based on the computed torque from Approach #2 and PD feedback control law.
function dydt = odefcn(t, y)
dydt = zeros(2,1);
% desired angular position profile
yd = ((120*pi/180)/3)*(t - (3/(2*pi))*sin((2*pi/3)*t));
% desired angular velocity profile (time derivative of yd)
dyd = (2/9)*pi*(1 - cos((2*pi/3)*t));
% desired angular acceleration profile (time derivative of dyd)
ddyd = (4/27)*(pi^2)*sin((2*pi/3)*t);
% parameters
m = 1; % mass of robot arm
g = 9.81; % gravity
L = 1; % length of robot arm
l = 0.5; % length of link
Iz = (1/3)*m*L^2; % moment of inertia
omega = 6; % natural frequency
zeta = 1; % damping ratio
kp = omega^2; % fundamental proportional control gain
kd = 2*zeta*omega; % fundamental derivative control gain
% Control Torque (3 approaches)
tau1 = Iz*ddyd + m*g*l*cos(yd);
tau2 = Iz*ddyd + m*g*l*cos(y(1));
tau3 = Iz*ddyd + m*g*l*cos(y(1)) - Iz*kp*(y(1) - yd) - Iz*kd*(y(2) - dyd);
dydt(1) = y(2);
dydt(2) = (tau1 - m*g*l*cos(y(1)))/Iz;
end
Results
tspan = linspace(0, 3, 3001)';
y0 = [120*pi/180; 0]; % initial condition (based on Problem 3.a)
[t, y] = ode45(@(t,y) odefcn(t, y), tspan, y0);
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917299/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917304/image.png)
Approach #1 does not work because the initial joint angle is
, but
is computed from
when
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917309/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917314/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917269/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917324/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917329/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917334/image.png)
Despite both trajectories have the same pattern, Approach #2 is unsuccessful because there is an offset of
.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917339/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917344/image.png)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917349/image.png)
By introducing the PD control to correct the offset error, Approach #3 can drive the joint angle θ to follow the desired profile
after 1 second.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/917269/image.png)
More Answers (1)
Bjorn Gustavsson
on 7 Mar 2022
Following from the comments above it seems as if you have a smoothly varying time-dependent torque that you know at a set of points in time. For that you could use simple interpolation of the torque at any given point in time. To do that you need to modify your ODE-function to allow for both a torque_of_t and a t_for_torque variables, instead of a single torque-input. Perhaps something like this:
function dy = ODE_fcn(..., tau_of_t,t_for_tau)
tau_now = interp1(t_for_tau,tau_of_t, t,'pchip');
% Then plug in tau_now where you use tau.
% The rest should be fine
end
Or you could wrap the interpolation into a function-handle:
tau_fcn = @(t) interp1(t_for_tau,tau_of_t,t,'pchip');
% then send this into your ODE
% But then you need to change the use of tau to a call to the function:
% Inside ODE:
tau_now = tau(t);
HTH
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!