How to get the solution to a GIVEN input VECTOR using ode45 solver

2 views (last 30 days)
Ivan on 31 Jul 2020
Commented: J. Alex Lee on 3 Aug 2020
Hi everyone!
I already posted my problem here, but unfortunately we did not come to a proper solution, so I've decided to reformulate my problem in an easier way, hoping it to be clear enough.
In the code below, you'll find the variable u_opt which represents the imput u of my system of equations.
u in function calc_only_f is a single scalar number, not a vector.
u_opt is a 50x1 double vector I obtained from a previous calculation as a solution of an optimal control problem (THIS INFO IS NOT RELEVANT). What I want to acheive is using u_opt as the imput of the ode45 solver and get for EACH component of this vector (50x1 i.e. for 50 different values) a 4x1 double solution as output.
After N=50 interation I would like to have a 50x4 double solution (aka trajectory) for those 50x1 double inputs of vector u_opt
As You can see below (represented by an arrow <-- in second to last line) I obtain for EACH component (i.e. value) of the input vector u_opt a 50x4 double solution. So for N=50 interations, I get 50 of these 50x4 double trajectoies, which are given as 1x50 double cell array.
Is there possibility to acheive exactly what I'm looking for, so getting a 50x4 double solution for those 50x1 double inputs? I spent a lot of hours in finding a solution, but unfortunately I'm still stuck on this problem.
........... calcf .........
% System of equations x' = f(t,x,u)
function f = calc_only_f(t,x,u,param)
% Parameter of the systems given as a struct
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
K = N*L; % L*(M + m*(1-Cx^2))
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)];
% +++++++++++++++++++++++++++++++++++++++++++++++++++
% -------------------- MAIN -------------------------
% +++++++++++++++++++++++++++++++++++++++++++++++++++
N=50; sizeu=1; tf= 5;
u_opt = rand(N*sizeu,1,'double'); % only for the implementation ... dimension of u_opt should be a 50x1 double
x0 = [0;0;pi;0]; % start point 4x1 double ... x0_1 position x0_2 velocity x0_3 angle x0_4 angular velocity
param.m = 1.5;
param.M = 27;
param.L = 2;
param.g = -9.80665;
param.d = 1.2;
param.N = 50;
param.timePart = linspace(t0,tf,N);
tspan = param.timePart;
for k = 1:numel(u_opt)
[t_ode{k},x_ode{k}] = ode45(@(t,x)calc_only_f(t,x,u_opt(k),param),tspan,x0); % <-- I get a 1x50 double cell array
I appreciate any help!
Cheers !
Ivan on 2 Aug 2020
Your answer:
"It's going to give you solutions at all the time steps you requested in tSpan, which in this case is 50 long." --> Yes and as I wrote previoulsly, I got a 1x50 cell array: each cell contains a 50x4 solution for each element of t_opt. This result does not help me because so I do not have anything to compare with my reference.
Yes, you got it right!
For each single value of u drawn from u_opt I want ode45 to produce a 1x4 solution (make sure not to confuse rows and columns ;) ). So for 50 calls, as You say, having a 50x4 double solution.
That's why I'm looking for a 50x4 solution for a given 50x1 input.
Considering only one time instant (for example at a time t_c) was only a hint for the community. This would help me at least to get a trajectory at that time instant t_c, let consider t_c to be at the half of tspan, for instance.
Maybe You or someothers know how to code it.

Sign in to comment.

Answers (1)

J. Alex Lee
J. Alex Lee on 2 Aug 2020
If you're just looking for code to extract some element out of your ode45 solution:
tAll = nan(1,length(u_opt));
xAll = nan(4,length(u_opt));
idx = 50; % or whatever index of tspan you want
for k = 1:numel(u_opt)
[t,x] = ode45(@(t,x)calc_only_f(t,x,u_opt(k),param),tspan,x0);
tAll(k) = t(idx);
xAll(:,k) = x(:,idx);
and just change idx to whatever you want. If you want t_c to be "half of tspan", and you have defined tspan with linspace, then let idx=25, for example.
J. Alex Lee
J. Alex Lee on 3 Aug 2020
Ok, glad it gets you closer to what you want.
But for would-be consumers of this "solution", be warned that in my understanding this was essentially a long protracted discussion about how to extract elements from a matrix, so that is all this answer purports to provide.

Sign in to comment.




Community Treasure Hunt

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

Start Hunting!