How to get the state trajectory with a GIVEN input VECTOR using ode45 solver

Hi everyone!

I would like to check whether my previous trajectory is plausible or not.

From the previous calculations I obtained a 50x1 double VECTOR I called u_opt , which represents the input (u) of the system of equations You see below.

function f = calcf(t,x,u,param)

% Physical parameter in the struct "param"
M = param.M;
g = param.g;
L = param.L;
m = param.m;
d = param.d;

Sx = sin(x(3)); % sin(x)
Cx = cos(x(3)); % cos(x)
N = M + m*(1-Cx^2); % denominator

f = [x(2); (1/N)*(-m*g*Cx*Sx + m*L*x(4)^2*Sx - d*x(2)+ u); x(4); (1/(N*L))*((m+M)*g*Sx -Cx*(m*L*x(4)^2*Sx - d*x(2)) - Cx*u)];

end

Now I would like to hand this vector u_opt over to the ode45 solver to obtain a 50x4 double x_ode that I can compare with my previous trajectory and see if it's plausible or not as follows:

[t_ode,x_ode] = ode45(@(t,x)calcf(t,x,u_opt,param),tspan,x0);

with tspan and x0 being just

tspan = linspace(0,tf,N); % N=50 - tf=5s ... final time
x0 = [0;0;pi;0]; % start state - x0_1 position , x0_2 velocity , x0_3 angle , x0_4 angular velocity

The Command window shows this error

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

Error using odearguments (line 95)
@(T,X)CALCF(T,X,U_OPT,PARAM)returns a vector of length 102, but the length of initial conditions vector is 4. The vector returned by
@(T,X)CALCF(T,X,U_OPT,PARAM) and the initial conditions vector must have the same number of elements.

The ode45 solver seems to work with a given input being a scalar, but I would like to calculate the solution of the diff. equation for every state at a time 't' being element of tspan, so that I ca compare my results / plots. That's why I want the vector u_opt to be the input.

Can anyone help me?

Thanks and cheers!

[Star Strider]

‘Can anyone help me?’

Only if you tell us what ‘param’ and ‘u’ are and what they contain. My guess is that one or more of them are vectors, and that is causing the problem.

[Ivan]

Thank You for the answer.

param - as already written in the comment - is a struct containing the scalar (numbers ... i.e. M=27, m=1.5, g=9.81, etc.) parameters of the system of equations and should not be a problem. param is only an elegant way to hand those parameters over to the function calcf.

Here one more example, maybe I explained myself not as clear as I thought.

u0 = 5; % with u being a scalar
[te,xe] = ode45(@(t,x)calcf(t,x,u0,param),tspan,x0);
% ----
% no problem at all ... ode45 returns a solution xe being a 50x4 double (4 because of the dimension of x0)

u - as you can see in the 3rd line at the top - represents the input of the system and it's a single number and NOT a vector.

It make sense to me that the ode45 solver does not accept u_opt because of the different data type.

As I said before, u_opt is a 50x1 vector. Maybe I have to give as an input every single component of this vector to the ode45 solver.

What I want to acheive is for every single component of the vector u_opt a single 1x4 double solution x_ode (see example in my 1st thread) for a single instant of tspan, so at the end for 50x1 single inputs 50x4 solutions.

i.e.: u_opt=[u1, u2, u3, .... , u49, u50]^T --> ode45 --> x_ode = 50x4 double

I hope that I was clear enough :)

Maybe You can help me to solve this problem

Answer 1

[Ayush Gupta]

ode45() is not suitable for finding multiple solutions, except through the mechanism of running multiple times with different u until the number of distinct solutions found has accumulated to the number desired. A loop can be run over the ode and the solution can be concatenated to get the desired solution. Refer to the following code: 

N=50; tf=5; 

% number  of partitions of the time interval, tf .. final time 
tspan = linspace(0,tf,N); 

tAll = zeros(length(u_opt),1); 

xAll = zeros(length(u_opt),4); 
% Could be whatever the index of tspan you want 
idx = 25; 

for k = 1:numel(u_opt) 

[tode,xode] = ode45(@(t,x)calc_f(t,x,u_opt(k),param),tspan,x0); 

tAll(k) = tode(idx); 

xAll(k,:) = xode(idx,:); 

end 

[Ivan]

Hi, thank You for the response.

Well, you're right in saying that ode45 is not suitable for finding multiple solutions for this type of settings.

The actual reason why I had troubles in getting my desired output, was the fact that I was trying to solve a continuos-time problem represented by the system dynamics calc_f by feeding the ode45 solver with the discrete-time vector u_opt.

In order to solve this mistake I interpolated the values in the input vector numerically and I "reconstruced" the required continuos trend.