ode45 with vector of inputs using symbolic toolbox

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 !!!

Answer 1

[Star Strider]

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:

[Japnit Sethi]

I have already used the lsim function and wanted to try out the ode45 solver because I want to eventually go and check the non linear equations I used to create this linear system of equations in future !!

Is there another way I can use the ode45 solver if I do not use the symbolic toolbox ?

Thanks

[Star Strider]

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]

So I went ahead and followed your approach to solve my system of equations with my input vector (function of time) using ode45.

My code is below:

% 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;


Eqns = @(t,v,r,psi,t_new,U_new) [PA_A1(1,1)*v + PA_A1(1,2)*r + PA_A1(1,3)*psi + PA_B1(1,1) * interp1(t_new(:), U_new(:), t);...
                                  PA_A1(2,1)*v + PA_A1(2,2)*r + PA_A1(2,3)*psi + PA_B1(2,1) * interp1(t_new(:), U_new(:), t)...
                                  PA_A1(3,1)*v + PA_A1(3,2)*r + PA_A1(3,3)*psi + PA_B1(3,1) * interp1(t_new(:), U_new(:), t)];

%Initial Conditions for v, r and psi
PSV1_0(1) = 0;  %v
PSV1_0(2) = 0;  %r
PSV1_0(3) = 0;  %psi

[T, Y] = ode45(@(t,v,r,psi) Eqns(t,v,r,psi,t_new,U_new), 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');

But I am getting an error:

Not enough input arguments.

Error in testing_10_21_2019>@(t,v,r,psi)Eqns(t,v,r,psi,t_new,U_new) (line 480)
[T, Y] = ode45(@(t,v,r,psi) Eqns(t,v,r,psi,t_new,U_new), t_new, PSV1_0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in testing_10_21_2019 (line 480)
[T, Y] = ode45(@(t,v,r,psi) Eqns(t,v,r,psi,t_new,U_new), t_new, PSV1_0);

[Star Strider]

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]

I am attaching the data file "actuator_node.fin" if that is the reason you are not able to run my code !!

What do you mean by r is not defined ? Also if i just use @(t,v) as you said above !! it throws me an error

Unrecognized function or variable 'r'.

Error in testing_10_21_2019>@(t,v)Eqns(t,v,r,psi,t_new,U_new) (line 480)
[T, Y] = ode45(@(t,v) Eqns(t,v,r,psi,t_new,U_new), t_new, PSV1_0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in testing_10_21_2019 (line 480)
[T, Y] = ode45(@(t,v) Eqns(t,v,r,psi,t_new,U_new), t_new, PSV1_0);

Just for clarification my state variables are v,r and psi and the input is U_new

[Star Strider]

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]

That is the same error I am getting which is because We did not define it as @(v,r,psi) in the function but still use it in the function like :

[T, Y] = ode45(@(t,v) Eqns(t,v,r,psi,t_new,U_new), t_new, PSV1_0);

Like over here in my Eqns equation I am still using r and psi but we do not define it for the ode45

[Japnit Sethi]

Now assuming I replace the 3 state variable v, r and psi by a single variable PSV1 like:

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);...
                                  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)...
                                  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)];

%Initial Conditions for v, r and psi
PSV1_0(1,1) = 0;  %v
PSV1_0(2,1) = 0;  %r
PSV1_0(3,1) = 0;  %psi

[T, Y] = ode45(@(t,PSV1) Eqns(t,PSV1,t_new,U_new), t_new, PSV1_0);  

I am getting a different error like :

Error using vertcat
Dimensions of arrays being concatenated are not consistent.

Error in
testing_10_21_2019>@(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);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),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)]
(line 471)
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);...

Error in testing_10_21_2019>@(t,PSV1)Eqns(t,PSV1,t_new,U_new) (line 480)
[T, Y] = ode45(@(t,PSV1) Eqns(t,PSV1,t_new,U_new), t_new, PSV1_0);

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in testing_10_21_2019 (line 480)
[T, Y] = ode45(@(t,PSV1) Eqns(t,PSV1,t_new,U_new), t_new, PSV1_0);
 

[Star Strider]

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]

I already said that there is no use of the expm function in my scenario as I am trying to solve the problem using an ode solver.

I have followed the workable method that you shared but it is causing some errors and I was simply asking help to solve the errors. I dont think I can explain the problem in any better way than I have up where I have posted the question !!

You and me were getting the same errors but I was simply asking of how to proceed in order to clear out the error which not sure if you know how to !!

https://www.mathworks.com/matlabcentral/answers/489355-ode45-with-vector-of-inputs-using-symbolic-toolbox#comment_763981

I appreciate the help though !! thanks for the time

[Star Strider]

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]

Appreciate the help !! That worked

[Star Strider]

As always, my pleasure!