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

6 次查看（过去 30 天）
Ivan
Ivan 2020-7-31
评论： J. Alex Lee 2020-8-3
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)];
end
% +++++++++++++++++++++++++++++++++++++++++++++++++++
% -------------------- 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
end
I appreciate any help!
Cheers !
  6 个评论
显示 4更早的评论
J. Alex Lee
J. Alex Lee 2020-8-1
Your original post:
"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."
It sounds to me that you do want each of 50 calls to ode45 (using a single value of u drawn from some vector, which in your words are irrelevant details to the problem) to produce a 4x1 solution. It's going to give you solutions at all the time steps you requested in tSpan, which in this case is 50 long. Do you want to choose some specific time step to pick out of the 4x50 solution that it will give you, as you hinted above with t_c? Which one? What does it mean to let ode45 "choose"?
Ivan
Ivan 2020-8-2
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.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（1 个）

J. Alex Lee
J. Alex Lee 2020-8-2
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);
end
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.
  5 个评论
显示 3更早的评论
Ivan
Ivan 2020-8-3
Ok nice! I post the solution for the community HERE :)
GOAL: get a 50x4 double solution from a 50x1 double array
At least we now get for at time instant idx a 50x4 double solution (i.e. a trajectory in control engineering for the input u_opt of the dynamic system) as follows:
N=50; tf=5; % set your own values
u_opt = rand(N,1,'double'); % Nx1 double vector I obtained from my previous calculation.
% I defined here a random vector so that you can run your code :)
tspan = linspace(0,tf,N); % N .. nr of partitions of the time interval, tf .. final time
tAll = zeros(length(u_opt),1);
xAll = zeros(length(u_opt),4);
idx = 25; % or whatever index of tspan you want
for k = 1:numel(u_opt)
[tode,xode] = ode45(@(t,x)calc_only_f(t,x,u_opt(k),param),tspan,x0);
tAll(k) = tode(idx);
xAll(k,:) = xode(idx,:);
end
That's partially what I was looking for, if I find a more dynamic approch, I'll post it below.
Thx guys!
J. Alex Lee
J. Alex Lee 2020-8-3
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.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by