ode45 with vector of inputs using symbolic toolbox
Show older comments
My system of equations are setup like:
So in my system I have a single input (delta_r/Rudder_Angle) but it is a vector of discrete constant values (U_new) with a time stamp (t_new) for each input value that is there is a new input value at each time stamp.
Also the "Actuator_node.fins" which is the data file with inputs and time is attached
As you can see in the image below:
I am using Symbolic toolbox for my code !!
clc
clear all
close all
syms v r psi delta_r
PSV1 = [v; r; psi]; % Yaw Axis State Variables vector
% A and B matrix
PA_A1 = [-0.736243229602963,-0.593872968776296,0;0.105400050677545,-0.887577743421107,0;0,1,0];
PA_B1 = [0.168078575433693;-0.565331828240409;0];
%% Input Data set and Time Span
Input_Data_2019_06_13 = readmatrix('actuators_node.fins.dat'); % Actuator fins data file contains the elevator angle and rudder angle input values in form of Starboards side and Bottom fins respectively
Input_Data_2019_06_13(:,1) = Input_Data_2019_06_13(:,1) - Input_Data_2019_06_13(1,1) ; % Trying to generalise the time stamp
% Convert the Inputs from degree to radians since the whole script model is
% in radians
Rudder_Angle = (-Input_Data_2019_06_13(:,4) + Input_Data_2019_06_13(:,5))/2; % Rudder angle is the subtraction of 2 fins divided by 2
Rudder_Angle = Rudder_Angle(768:2616)*(pi/180); % Selecting Input data when the fins start changing i.e. after the thrust is being provided
U = Rudder_Angle; % Our Simulation input for Yaw Axis Model is the Rudder Angle
U(isnan(U)) = 0; % To convert all the NaN elements in U matrix to 0
% Note: Time span is in accordance to the input data set, i.e. at each time
% stamp we have a input data
t_old = transpose(linspace(0,100,length(Input_Data_2019_06_13(768:2616,:)))); % t = linspace(0, 10, 100); % Want the time stamp of Input to be from 0 to 100, with no of input values be equal to the leght of input data we have selected
t_new = transpose(linspace(0,100,length(Input_Data_2019_06_13(768:2616,:)))) ;
U_new = interp1(t_old,U,t_new);
U_new(isnan(U_new)) = 0;
% Setting up first order differential equations v' and r' and psi'
dv = diff(v);
dr = diff(r);
dpsi = diff(psi);
eqn1 = dv == PA_A1(1,1)*v + PA_A1(1,2)*r + PA_A1(1,3)*psi + PA_B1(1,1) * U_new;
eqn2 = dr == PA_A1(2,1)*v + PA_A1(2,2)*r + PA_A1(2,3)*psi + PA_B1(2,1) * U_new;
eqn3 = dpsi == PA_A1(3,1)*v + PA_A1(3,2)*r + PA_A1(3,3)*psi + PA_B1(3,1) * U_new;
%Initial Conditions for v, r and psi
PSV1_0(1) = 0; %v
PSV1_0(2) = 0; %r
PSV1_0(3) = 0; %psi
% Reducing a higher order differential equation to an equivalent first
% order differential equations
[eqs, vars] = reduceDifferentialOrder([eqn1, eqn2, eqn3], [v, r, psi]);
% Converting differential equation to mass matrix form
% Way of rewriting the equations in form of variables and the constants
% being inserted, Note 'f' gives the final v' and r' adn psi' eq after computation
[M,F] = massMatrixForm(eqs,vars);
f = M\F;
% odeFunction is used to convert symbolic expressions in our case 'f' to
% function handle in our case 'odefun' to be used for ODE solvers
odefun = odeFunction(f,vars);
disp(odefun)
% Using ode45 in time period of 0 to 140 second and with initial conditions
% ic
[t,y] = ode45(odefun, t_new, PSV1_0);
% Plotting
figure
subplot(2,2,1)
plot(t_new,y(:,1),'-r');
grid on;
ylabel('v(t) in m/s');
xlabel('t in sec');
subplot(2,2,2)
plot(t_new,y(:,2),'-r');
grid on;
ylabel('r(t) in m/s');
xlabel('t in sec');
subplot(2,2,3)
plot(t_new,y(:,3),'-r');
grid on;
ylabel('psi(t) in m/s');
xlabel('t in sec');
I am trying to find how to put those input values in my code since my input are a bunch of discrete constant values at each time ?
I want to eventually plot v(t) vs t, r(t) vs t and psi(t) vs t graphs !!!
Accepted Answer
Star Strider
on 5 Nov 2019
The Symbolic Math Toolbox is not the best approach to this straightforward linear control problem. If you do not want to use the Control System Toolbox, just use the expm function and a loop:
A = [-0.736243229602963,-0.593872968776296,0;0.105400050677545,-0.887577743421107,0;0,1,0]; % Define State Space Matrices
B = [0.168078575433693;-0.565331828240409;0]; % Define State Space Matrices
C = eye(3); % Define State Space Matrices
D = zeros(3,1); % Define State Space Matrices
t = linspace(0, 5, 50); % Use Your Own Time Vector
uv = exp(-linspace(0, 5, 50)); % Use Your Own Input Vector
for k = 1:numel(t)
y(:,k) = (C*expm(A*t(k))*B + D)*uv(k);
end
figure
plot(t, y)
grid
legend('\nu', 'r', '\psi', 'Location','SE')
Experiment to get different results.
The assignment for ‘y’ is derived from:
Laplace transform:
inverting:
combining:
13 Comments
Japnit Sethi
on 5 Nov 2019
Star Strider
on 5 Nov 2019
The Symbolic Toolbox will make these easier. However I do not know what nonlinear equations you want to use with it, so I cannot suggest a specific spproach. I most often use odeToVectorField and then matlabFunction to create anonymoius functions to use with the ODE solvers, particularly ode45. Different numeric ODE solvers are used for different problems.
If you have a particular input vector (function of time) that you want to use with the numeric ODE solvers, see the documentation section on: ODE with Time-Dependent Terms to understand how to use it with your ODE anonymous function or function file. For an illustration of how to do this with an anonymous function ode function, see: ODE Solver with array value parameter. It is necessary to do this because the MATLAB ode solvers are adaptive, not fixed-step, so your ODE function will need to know the value of your input function at whatever time ode45 integrates it, not only the values in the vector.
Japnit Sethi
on 5 Nov 2019
Star Strider
on 5 Nov 2019
The ode45 function only wants ‘t’ and ‘v’.
Try this:
[T, Y] = ode45(@(t,v) Eqns(t,v,r,psi,t_new,U_new), t_new, PSV1_0);
I cannot run your code. Note that ‘r’ is not defined, and neither are most of the other constants.
Japnit Sethi
on 5 Nov 2019
Star Strider
on 5 Nov 2019
I downloaded the file and opened it as in your code.
When I ran it (with the file), it threw:
Unrecognized function or variable 'r'.
witrh the ode45 call.
Japnit Sethi
on 5 Nov 2019
Japnit Sethi
on 5 Nov 2019
Star Strider
on 5 Nov 2019
I do not understand what you are doing. I have no idea what system you want to integrate, except the one you originally described.
I already provided a workable method using expm to integrate them, as well as odeToVectorField and matlabFunction to use with the numerical ODE solvers, including a method to use your input data file. You can certainly create a ODE function that the ODE solvers can then use to integrate your system of ordinary differential equations.
It should be fairly straightforward to define a system of ordinary differential equations, then use the appropriate Symbolic Math Toolbox functions to create anonymous functions or function files that the ODE solvers can use to integrate them.
Japnit Sethi
on 5 Nov 2019
Edited: Japnit Sethi
on 5 Nov 2019
Star Strider
on 5 Nov 2019
Edited: Star Strider
on 5 Nov 2019
The problem appears to be in the continuation ellipses in ‘Eqns’. When I eliminated them, it worked.
Try this:
Eqns = @(t,PSV1,t_new,U_new) [PA_A1(1,1)*PSV1(1) + PA_A1(1,2)*PSV1(2) + PA_A1(1,3)*PSV1(3) + PA_B1(1,1) * interp1(t_new(:), U_new(:), t, 'linear', 'extrap');
PA_A1(2,1)*PSV1(1) + PA_A1(2,2)*PSV1(2) + PA_A1(2,3)*PSV1(3) + PA_B1(2,1) * interp1(t_new(:), U_new(:), t, 'linear', 'extrap');
PA_A1(3,1)*PSV1(1) + PA_A1(3,2)*PSV1(2) + PA_A1(3,3)*PSV1(3) + PA_B1(3,1) * interp1(t_new(:), U_new(:), t, 'linear', 'extrap')];
I also included extra arguments to interp1 in case your code attempted to ‘colour outside the lines’, and request values that are not in your data. They may not be necessary, however interp1 will not use them if it does not need to.
EDIT —
The plot —
Japnit Sethi
on 6 Nov 2019
Star Strider
on 6 Nov 2019
As always, my pleasure!
More Answers (0)
Categories
Find more on Spline Postprocessing in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
