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

16 views (last 30 days)
Ivan on 29 Jul 2020
Edited: Ivan on 10 Sep 2020
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)];
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!
Ivan on 29 Jul 2020
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

Sign in to comment.

Answers (1)

Ayush Gupta
Ayush Gupta on 10 Sep 2020
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,:);
  1 Comment
Ivan on 10 Sep 2020
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.

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!