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

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 !

[J. Alex Lee]

It is clearer than in your last question that you expect "a 4x1 double solution as output" of ode45 for any given value of u in u_opt. But it's still (maybe more so) unclear what you really want. Are you looking for just the terminal value of each variable, at time = tf?

By the way it doesn't help you that you happen to have 50 time steps and also u_opt happens to be 50.

[Ivan]

(#1) is a place holder for "new" trajectory, so that You know, which one I'm considering, but it should be trivial

Hi Alex.

As I already wrote in the 1st post, u_opt is a 50x1 double vector I obtained from a previous calculation as a solution of an optimal control problem, which represents the force which has to be applied on the cart in order to obtain the best (optimal) trajectory.

What I'm looking for is to check the plausibilty of this input vector. I want to hand u_opt over to my system of equations (described by the function calc_only_f) and determine by using ode45 a "new" (#1) solution for these values of the force. That "new" solution can be seen as a new feasible trajectory, which I'll compare to the optimal trajectory calculated (using fmincon) before.

Are you looking for just the terminal value of each variable, at time = tf (tf = 5 ... the final time is just a number)

No, I'm not. Each of the values in the input vector u_opt were determined at different time instant in tspan, as follows:

t=t1 ... u_opt(t1) = u_opt1

t=t2 ... u_opt(t2) = u_opt2

...

t=N ... u_opt(tN) = u_optN

Maybe You can help me further.

(if You don't know how to implement it or if it's not possible at all ... please just read this text below)

Since this request seems to be not realisable, I could determine the "new" (#1) trajectory by selecting just ONE arbitrarily time instant, let us call it t_c and for this time t_c let the ode45 determine the "new" trajectory for that single element of u_opt (i.e. a 1x4 double )and after N=50 interations we should be able to find a "new" trajectory calculated at time t_c (50x4 double) and check if at the observed time t_c the "new" state (also a single value) of the trajectory is equal/ similar (surely considering numerical approximations) or different to the "optimal" one calculated by fmincon.

Maybe You can give me a hint in how to implement this. I've tried by myself as You seen, now I'm stuck since a couple of days.

Thx

[J. Alex Lee]

As far as I understand, it's a non-starter as soon as you assert you want ode45 to return a [4x1] matrix upon giving it a 1x50 sized time (domain) vector.

ode45 has to return something at least 2 elements long and M-elements tall, where M is the number of equations you specify in your odefun (or the transpose, i'm continually confused about how matlab likes to return columns and rows).

If you can revise this part, maybe the community can invest a bit more in understanding what you actaully want (which I am afraid I still cannot grasp).

[Ivan]

Actually I want ode45 to return a 50x4 double upon giving it a 50x1 sized vector.

Well I reformulated the problem twice, I tried to explain it in an easier way, that's why I opened a second thread, furthermore I gived to the community additional informations about the goal of my task (as last part of my final bachelor thesis in control engineering/ optimisation at university) for some background informations.

What I want to acheive is >> I want ode45 to return a 50x4 double upon giving it a 50x1 sized vector <<

ode45 has to return something at least 2 elements long and M-elements tall, where M is the number of equations you specify in your odefun (or the transpose, i'm continually confused about how matlab likes to return columns and rows).

> Yes, you're right. That's why I set the start point x0 to be a 1x4 double sized vector. (see code) In this case 4 represents the nr. of colums which are equal to the 4 state of my inverse pendulum's system (I've added several comments in my code, so it should be clear enough).

Maybe You can share it with the community, since I do not have such a visibilty as You have.

I really hope, that someone can help me. I did all the rest by myself (it was not that trivial, but I could managed it), now I'm a bit stuck :(

Looking forward to hearing Your thoughs.

Bye

[J. Alex Lee]

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]

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.

Answer 1

[J. Alex Lee]

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.

[Ivan]

Thank Your for Yours solution.

Why did you invert columns amd rows in the definitions of tAll and xAll?

I ran both versions, Yours and mine, unfortunately we get the same mistake

% My version
tAll = zeros(length(u_opt),1); % initializing the vectors using zeros(...) and NaN(...)
xAll = zeros(length(u_opt),4);
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);    
    u_opt(k) % Added to see if we jump correctly in the loop --> OK
    u_opt(idx) % Same --> OK
    tAll(k) = t(idx);
    xAll(:,k) = x(:,idx);
end

I get the following error from the command window

Index in position 2 exceeds array bounds (must not exceed 4).
Error in main (line 203)  xAll(:,k) = x(:,idx);

Only the 1st cycle has been executed

u_opt(k) % k=1 .. uopt(k) has a certain value, which corresponds
u_opt(idx) % Same 

For k=2 the process interrupts and we get the error I mentioned above.

[Walter Roberson]

x(idx, :)

[J. Alex Lee]

I mistakenly thought ode45 returns rows rather than columns (I think because I've started using the struct output, and therein the solutions are returned as rows...thanks TMW for making things easy). Your transposition and Walter's fix partially makes sense, I guess you want

xAll(k,:) = x(idx,:)

By the way, my case-in-point about the confusion of the equality of lengths of your time vector and u_opt vector is evident here...there is no reason to want (and I think you have admitted as much) idx==k.

[Ivan]

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]

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.