ode45 with a vector of inputs with respect to time stamp

I previously used lsim for my simulation but I wanted to try out ode45 this time.

I having a problem using ode45 since I have a discrete set of inputs and time stamp i.e. I have an input for each time stamp. I tried looking up on Matlab Answers and tried integrating it into my code but still have errors.

Here is my Code:

% In my system I have an A and B matrix, my state variable vector which is a vector of 3 state variables, my input is just 1 which is delta_r also called Rudder Angle

syms v r psi delta_r
PSV1 = [v; r; psi]; % Yaw Axis State Variables vector
PSV1_0 = [0; 0; 0];

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

%% Time
t_new = transpose(linspace(0,100,length(Input_Data_2019_06_13(768:2616,:)))) ;

%% Integration Function
tResult = [];
yResult = [];
for index = 2:length(t_new)
  % Integrate:
  delta_r = Rudder_Angle(index - 1);
  t_new = t_new(index-1:index);
  [t_new, y_ode45] = ode45(@(v, r, psi)  AUV_ode(PA_A1, PA_B1, PSV1, delta_r), t_new, PSV1_0);
    % Collect the results:
    tResult = cat(1, tResult, t_new);
    yResult = cat(1, yResult, y_ode45);
    % Final value of x is initial value for next step:
    PSV1_0(1) = y_ode45(end, 1);
    PSV1_0(2) = y_ode45(end, 2);
    PSV1_0(3) = y_ode45(end, 3);
end

figure
subplot(2,2,1)
plot(tResult,yResult(1));
subplot(2,2,2)
plot(tResult,yResult(2));
subplot(2,2,3)
plot(tResult,yResult(3));


function Yaw_Axis_Equations = AUV_ode(A_Matrix, B_Matrix, State_Variable, Input_data)
Yaw_Axis_Equations = A_Matrix*State_Variable + B_Matrix*Input_data;                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  );
end

[Japnit Sethi]

The error that I am currently getting is:

Error: File: testing_10_21_2019.m Line: 51 Column: 806
Invalid expression. When calling a function or indexing a variable, use parentheses. Otherwise, check for mismatched delimiters.

[Walter Roberson]

Scroll way to the right on your second last line. You will see you have an extra

);

there.

[Japnit Sethi]

Thanks !! I couldnt see that at all !!

Answer 1

[Mani Teja]

Hello Japnit,

From your question, I am not sure what kind of errors are occuring. After running the program, I got this error: "Inputs must be floats, namely single or double". If thats the case, this could be helpful:

https://in.mathworks.com/matlabcentral/answers/399165-inputs-must-be-floats-namely-single-or-double-ode45-function

Regards

Mani

[Walter Roberson]

  [t_new, y_ode45] = ode45(@(v, r, psi)  AUV_ode(PA_A1, PA_B1, PSV1, delta_r), t_new, PSV1_0);

PSV1 is a symbolic vector of variable names

  PSV1 = [v; r; psi]; % Yaw Axis State Variables vector

and therefore

  Yaw_Axis_Equations = A_Matrix*State_Variable + B_Matrix*Input_data;

is going to be a symbolic result involving unresolved variable names. This is fundamentally incompatible with ode45() use, which can only work with numeric outputs as it is a numeric solver.

[Japnit Sethi]

I made the appropriate changes as suggested by you.

I removed the syms(), and also the line delta_r

%% Removed lines
% syms v r psi delta_r
% PSV1 = [v; r; psi];
% delta_r = Rudder_Angle(index - 1);

The changes are:

Rudder_Angle = Rudder_Angle(index - 1);
[t_new, y_ode45] = ode45(@(v, r, psi)  AUV_ode(PA_A1, PA_B1, Rudder_Angle), t_new, PSV1_0);

function Yaw_Axis_Equations = AUV_ode(A_Matrix, B_Matrix, Input_data)
Yaw_Axis_Equations{1} = @(v,r,psi) A_Matrix(1,1)*v + A_Matrix(1,2)*r + A_Matrix(1,3)*psi + B_Matrix(1)*Input_data;
Yaw_Axis_Equations{2} = @(v,r,psi) A_Matrix(2,1)*v + A_Matrix(2,2)*r + A_Matrix(2,3)*psi + B_Matrix(2)*Input_data;
Yaw_Axis_Equations{3} = @(v,r,psi) A_Matrix(3,1)*v + A_Matrix(3,2)*r + A_Matrix(3,3)*psi + B_Matrix(3)*Input_data;
end

But now the error I am getting I am not sure whats wrong:

Error using odearguments (line 93)
@(V,R,PSI)AUV_ODE(PA_A1,PA_B1,RUDDER_ANGLE) must return a column vector.

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

Error in testing_10_21_2019 (line 35)
  [t_new, y_ode45] = ode45(@(v, r, psi)  AUV_ode(PA_A1, PA_B1, Rudder_Angle), t_new, PSV1_0);

[Walter Roberson]

function Yaw_Axis_Equations = AUV_ode(A_Matrix, B_Matrix, Input_data)
Yaw_Axis_Equations{1,1} = @(v,r,psi) A_Matrix(1,1)*v + A_Matrix(1,2)*r + A_Matrix(1,3)*psi + B_Matrix(1)*Input_data;
Yaw_Axis_Equations{2,1} = @(v,r,psi) A_Matrix(2,1)*v + A_Matrix(2,2)*r + A_Matrix(2,3)*psi + B_Matrix(2)*Input_data;
Yaw_Axis_Equations{3,1} = @(v,r,psi) A_Matrix(3,1)*v + A_Matrix(3,2)*r + A_Matrix(3,3)*psi + B_Matrix(3)*Input_data;
end

However, the function you invoke must return numeric values, not function handles.

You are trying to manipulate families of equations, and that does not work in any of the ode* routines. You can only work with specific equations in which everything is pinned down except for t and the current boundary conditions for the various derivatives.

[Japnit Sethi]

I have put the delta_r back since there is no use to make it a symbolic variable and since my Rudder_Angle was getting messed up

%% Integration Function
tResult = [];
yResult = [];
for index = 2:length(t_new)
  % Integrate:
  delta_r = Rudder_Angle(index - 1);
  t_new = t_new(index-1:index);
  [t_new, y_ode45] = ode45(@(v, r, psi)  AUV_ode(PA_A1, PA_B1, delta_r), t_new, PSV1_0);
    % Collect the results:
    tResult = cat(1, tResult, t_new);
    yResult = cat(1, yResult, y_ode45);
    % Final value of x is initial value for next step:
    PSV1_0(1) = y_ode45(end, 1);
    PSV1_0(2) = y_ode45(end, 2);
    PSV1_0(3) = y_ode45(end, 3);
end
function Yaw_Axis_Equations = AUV_ode(A_Matrix, B_Matrix, Input_data)
Yaw_Axis_Equations{1,1} = @(v,r,psi) A_Matrix(1,1)*v + A_Matrix(1,2)*r + A_Matrix(1,3)*psi + B_Matrix(1)*Input_data;
Yaw_Axis_Equations{2,1} = @(v,r,psi) A_Matrix(2,1)*v + A_Matrix(2,2)*r + A_Matrix(2,3)*psi + B_Matrix(2)*Input_data;
Yaw_Axis_Equations{3,1} = @(v,r,psi) A_Matrix(3,1)*v + A_Matrix(3,2)*r + A_Matrix(3,3)*psi + B_Matrix(3)*Input_data;
end

I would assume its taking a specific value of Rudder_Angle everytime for my input delta_r and since I do not have any more symbolic variables, not sure why am i still getting this error:

Error using odearguments (line 113)
Inputs must be floats, namely single or double.

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

Error in testing_10_21_2019 (line 35)
  [t_new, y_ode45] = ode45(@(v, r, psi)  AUV_ode(PA_A1, PA_B1, delta_r), t_new, PSV1_0);

[Walter Roberson]

@(v, r, psi)  AUV_ode(PA_A1, PA_B1, delta_r)

ode45 passes only two variables to the function you give (unless you use an obsolete syntax that has not been documented for close to 20 years.) The first variable a scalar that is the current time (or analog), and the second is the vector of boundary conditions. You are instead expecting three variables to be passed in, and you are ignoring all three of them. You are also attempting to return a cell array of function handles, which is not permitted.

Now, it is, generally speaking, valid to call a function that constructs a function handle that is to be passed into ode45. It would be valid, for example, to have

F = construct_AUV_ode(PA_A1, PA_B1, delta_r);
[t_new, y_ode45] = ode45(F,  t_new, PSV1_0);

in conjunction with

function Yaw_Axis_Equations = construct_AUV_ode(A_Matrix, B_Matrix, Input_data)
Equations{1} = @(v,r,psi) A_Matrix(1,1)*v + A_Matrix(1,2)*r + A_Matrix(1,3)*psi + B_Matrix(1)*Input_data;
Equations{2} = @(v,r,psi) A_Matrix(2,1)*v + A_Matrix(2,2)*r + A_Matrix(2,3)*psi + B_Matrix(2)*Input_data;
Equations{3} = @(v,r,psi) A_Matrix(3,1)*v + A_Matrix(3,2)*r + A_Matrix(3,3)*psi + B_Matrix(3)*Input_data;
Yaw_Axis_Equations = @(v, r, psi) deal( [Equations{1}(v, r, psi); Equations{2}(v, r, psi); Equations{3}(v, r, psi)], []);
end

however, it is not obvious what psi would be here; for that matter it is not obvious that tnew corresponds to v, or that the boundary conditions correspond to r.

My speculation is that you would want more like

function Yaw_Axis_Equations = construct_AUV_ode(A_Matrix, B_Matrix, Input_data)
Equations{1} = @(v,r,psi) A_Matrix(1,1)*v + A_Matrix(1,2)*r + A_Matrix(1,3)*psi + B_Matrix(1)*Input_data;
Equations{2} = @(v,r,psi) A_Matrix(2,1)*v + A_Matrix(2,2)*r + A_Matrix(2,3)*psi + B_Matrix(2)*Input_data;
Equations{3} = @(v,r,psi) A_Matrix(3,1)*v + A_Matrix(3,2)*r + A_Matrix(3,3)*psi + B_Matrix(3)*Input_data;
Yaw_Axis_Equations = @(t, vrp) deal( [Equations{1}(vrp(1), vrp(2), vrp(3)); Equations{2}(vrp(1),vrp(2),vrp(3)); Equations{3}(vrp(1),vrp(2),vrp(3))], []);
end

Also, be careful:

  t_new = t_new(index-1:index);

replaces t_new with a two-element subset of t_new . You want a separate variable for the current timespan, such as t_cur = t_new(index-1:index); that you would pass in as the time range.