Why ode45 is not working??

2021 7 19
1 个回答
更新时间：2021 7 22
11 次查看（30 天）

0 个投票

Hello,
I've uploaded the two following files:
1. epas_steering_system.m is the main file, and it is the file where ode45 gives me an error,
2. epasFun.m is the function file.
My question is, does somebody know why ode45 won't work?
Thanks.

35 个评论
显示 33更早的评论 隐藏 33更早的评论

Torsten
Torsten 2021-7-19
编辑：Torsten 2021-7-19
If the two files are not too big in size, maybe you could include them here as plain text. I have no possibility to download files at the moment.
Also including the complete error message would be of great help.
Frane
Frane 2021-7-19
I've uploaded the codes in text files and the error message.
Torsten
Torsten 2021-7-19
编辑：Torsten 2021-7-19
1 The definition of I from Ms must be done elementwise:
for i = 1:numel(time)
if Ms(i) >= 3.5
I(i) = ...
....
end
2 Ms(x), I(x) and Fr(x) in the list of parameters are incorrect.
Make a loop when calling the ode solver over the number of elements in the array time:
for i = 1:numel(time)
[T{i},X{i}] = ode45(@(t,x) epasFun(x,...,Ms(i),I(i),Fr(i)),...
end
Frane
Frane 2021-7-19
编辑：Frane 2021-7-19
I've changed the things you wrote (I hope I changed them properly). The code is uploaded in a text file so you can check it. But it still won't work. Maybe I've done it wrong?
Can you explain to me what does numel mean/represent in the for loop?
And one more thing, why can't I see the "I" for the given "Ms"? I can write out all the "Ms" but can't write out "I"?
Ankit
Ankit 2021-7-19
numel: number of array elements (link)
Torsten
Torsten 2021-7-19
Again look at the necessary changes to be made I indicated in my response.
You didn't incorporate them properly in your code.
Frane
Frane 2021-7-19
I got the part for Ms and I okey now as seen in the following (now it gives me the appropiate I for every Ms):
for i = 1:numel(time)
if Ms(i) >=3.5
I(i) = 21.29 * Ms(i) - 69.4;
elseif Ms(i) < 3.5 & Ms(i) >= 2
I(i) = 2.73 * Ms(i) - 4.47;
elseif Ms(i) < 2 & Ms(i) >= 0
I(i) = 0.5 * Ms(i);
elseif Ms(i) < 0 & Ms(i) >= (-2)
I(i) = 0.5 * Ms(i);
elseif Ms(i) < (-2) & Ms(i) >= (-3.5)
I(i) = 2.73 * Ms(i) + 4.47;
elseif Ms(i) < (-3.5)
I(i) = 21.29 * Ms(i) + 69.4;
end
end
But the part for the ode45 is written wrong. I just don't know what? Can you write it in the right way? :/
for i = 1:numel(time)
tspan = 0:0.1:3;
x0 = [0; 0; 0; 0; 0; 0; 0; 0; 0];
[t(i),x(i)] = ode45(@(t,x)epasFun(x,Juc,Jlc,kuc,klc,duc,dlc,rpin,mr,dm,Jm,K,Ms(i),I(i),Fr(i)), tspan, x0);
end
Torsten
Torsten 2021-7-19
编辑：Torsten 2021-7-19
Again check the call to ode45 with my response. Note that {} and () are different in Matlab.
Frane
Frane 2021-7-19
Still won't work. Now I changed the () to {}.
for i = 1:numel(time)
[T{i},X{i}] = ode45(@(t,x)epasFun(x,Juc,Jlc,kuc,klc,duc,dlc,rpin,mr,dm,Jm,K,Ms(i),I(i),Fr(i)), tspan, x0);
end
Torsten
Torsten 2021-7-19
What is the error message ?
Torsten
Torsten 2021-7-19
编辑：Torsten 2021-7-19
One error is that you must replace Fr(i) by Fr in the list of parameters handed to epasFun because Fr=1 is a scalar.
If the error persists, check the size of all variables handed to epasFun in this function. One of them must have size 3x9 instead of 1. So insert the lines
size(x)
size(Juc)
size(Jlc)
etc
at the beginning of function epasFun.
Frane
Frane 2021-7-19
I have founded the variable. It's "K".
But I've set the variable "K" in the begining to be 0.01. How is it 3x9?
Torsten
Torsten 2021-7-20
You overwrote the scalar value for K by K = lqr(A,B,Q,R). Use another name here.
Frane
Frane 2021-7-20
Yeah it works now. If I want to plot the T and X, I can't use plot. What can I use to see the variables? Because in an example I've seen online they used plot. Here is the following code and the plot after it.
figure
plot(t,x,'LineWidth',2);
legend('x', 'v', 'theta', 'omega');
xlabel('Vrijeme');
ylabel('Stanje');
set(gcf,'color','w');
title('Odaziv svih varijabli sustava sa LQR');
Torsten
Torsten 2021-7-20
I'm not completely sure how to deal with the cell arrays T and X. So to keep control, change the last part of the code to
tspan = 0:0.1:3;
x0 = zeros(9,1);
X = zeros(numel(time),numel(tspan),9);
for i = 1:numel(time)
[t,x] = ode45(...);
X(i,:,:) = x(:,:);
end
Now you can plot, e.g. all nine components of the solution for Ms(1) and I(1) :
plot (t,X(1,:,:))
Frane
Frane 2021-7-20
I get the following error message.
Torsten
Torsten 2021-7-20
plot (t,squeeze(X(1,:,:)))
Frane
Frane 2021-7-20
So the code is working. I've uploaded the new code here in text file. Can you explain to me what am I looking at, and by that I mean to these lines represent Ms, I and Fr or something else?
Torsten
Torsten 2021-7-20
Plot shows the 8th solution variable over time for (Ms(1),I(1)), (Ms(2),I(2)),...,Ms(numel(time)),I(numel(time))) (thus there should be 20 curves because numel(time) = 20, I guess).
Frane
Frane 2021-7-20
编辑：Frane 2021-7-20
Yeah the X is 20x3x9. Can I show all the results at the same time or one by one like this? There are different result (example for X(:,:6), x(:,:,7)) and some are 0.
Even by changin tspan I can only see 12 of them.
Frane
Frane 2021-7-20
Got this with this time. How can I zoom it up, because when I change the time frame I either get the pictures above or something like this?
Torsten
Torsten 2021-7-20
help xlim
help ylim
Frane
Frane 2021-7-20
编辑：Frane 2021-7-20
You wrote that there shoud be 20 curves because numel(time) = 20 and the time is:
time = 0.01: 0.01: 0.2;
So why 20 curves? I don't understand? Can you explain to me?
Thanks.
P.S. Do these curve respresent the outputs given in the left side of the equation in the picture below (even though there are 9 outputs)?
Torsten
Torsten 2021-7-20
If you now tell me that Ms and I are interpolation curves dependent on t and that the values Ms(t) and I(t) have to be inserted in the equations when the solver has reached time t, we can start anew modifying your code.
Frane
Frane 2021-7-20
编辑：Frane 2021-7-20
Just to explain Ms and I:
Ms - represents the momentum given to the steering wheel
I - is the current given by the electromotor to help steer, depending on the momentum Ms
The part in the beggining is from the picture below it (hope this answers you question above?).
P.S. Ms, I and Fr are inputs (where Fr is given by the CarMaker; that's why it's value is 1).
for i = 1:numel(time)
if Ms(i) >=3.5
I(i) = 21.29 * Ms(i) - 69.4;
elseif Ms(i) < 3.5 & Ms(i) >= 2
I(i) = 2.73 * Ms(i) - 4.47;
elseif Ms(i) < 2 & Ms(i) >= 0
I(i) = 0.5 * Ms(i);
elseif Ms(i) < 0 & Ms(i) >= (-2)
I(i) = 0.5 * Ms(i);
elseif Ms(i) < (-2) & Ms(i) >= (-3.5)
I(i) = 2.73 * Ms(i) + 4.47;
elseif Ms(i) < (-3.5)
I(i) = 21.29 * Ms(i) + 69.4;
end
end
Torsten
Torsten 2021-7-20
And why do these quantities depend on time ?
Frane
Frane 2021-7-20
编辑：Frane 2021-7-20
Because the driver rotates the wheel differently during a period of time, and depending on the momentum given by the driver (Ms), a current (I) is sent to help the driver to turn the wheel (so he doesn't use too much "force").
Could that be the reason otherwise I don't know. I was given by the proffesor like that.
P.S. When I change the "time" from
time = 0.01: 0.01: 0.2;
to
time = 0.01: 0.01: 0.1;
I get less curves (and the other way around).
Torsten
Torsten 2021-7-20
At the moment, the equations are solved for Ms and I being fixed parameter values over the integration time. Thus for each combination for Ms and I (and there exist as many combinations as the vector "time" has elements), you get 9 solution curves over the time period defined by "tspan". I suspect that the time instants specified in "time" and "tspan " must be identical and that the values of Ms and I must be interpolated to the integration time provided by ODE45 in function epasFun. But this is not what is simulated in the actual code.
Frane
Frane 2021-7-20
Yes when "time" and "tspan" are the same I get 9 curves representing different "Ms" and "I".
I'm doing this based on an example of inverted pendulum on cart (I've uploaded the files). In lqr_pendcart the guy made the controler "u" that in this case controled the force by which the cart was being pulled or pushed. There is only one input (the force F represented by u) and four outputs given in a plot when you run the code (also you get two stepplots one before using LQR and one after). So my task was to make a controler in a similar way that they did with the inverted pendulum on cart, and I needed to get the response of the ouput variables similar like in the given example.
Frane
Frane 2021-7-21
Can I ask you one more thing? In my case i did a stepplot (code and picture is given below). I was wondering am I looking at the picture in the red rounded direction or the blue direction? I know it represents my 3 inputs but I don't know how to look at the plot?
figure
sys = ss(A, B, C, D); % ss = prostor stanja
t1 = 0:0.1:5;
stepplot(sys,t1);
set(gcf,'color','w');
xlabel('Vrijeme');
ylabel('Odziv');
hold on;
Torsten
Torsten 2021-7-21
Sorry, but I have no experience with the output of function "stepplot".
You should open a new question.
Frane
Frane 2021-7-21
So I went to my professor, he changed my code a little bit (some of the code that was in epas_steering_system is now in epasFun). He told me that I shoud do ode45 without for loop (given in the files below). But as you can see all of my variables divaricate, and they shoud converge. Do you maybe know what could be the problem? Maybe I've written some of the equations wrong in epasFun, or maybe the formula for Ms and I is wrong, or some of the parameters are wrong? But expect of the given problems could there be something else?
Torsten
Torsten 2021-7-21
According to the equations you posted, they seem to be implemented correctly in epasFun.
Of course, I can't tell whether they are correct or whether the parameters you chose make sense.
Frane
Frane 2021-7-22
I checked the parameters, they are okey; checked the equations, they are okey.
It could be something with the code itself (probably with the ode45 part). To be honest I don't know what, and I'm stuck at the moment.

请先登录，再进行评论。

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

回答（1 个）

Tesfaye Girma
Tesfaye Girma 2021-7-21

0 个投票

you can try this approach if it worked for you
f = (t,y) [ 4.86*y(3) - 4.86*10^14*y(1)*y(2);
4.86*y(3) - 4.86*10^14*y(1)*y(2);
-4.86*y(3) + 4.86*10^14*y(1)*y(2) ];
opt = odeset('maxstep',1e-13);
tspan = [0 1e-11];
y0 = [1.48e-8; 6.7608e-3; 1];
[t,y] = ode45(f,tspan,y0,opt);
subplot(311)
plot(t,y(:,1))
subplot(312)
plot(t,y(:,2))
subplot(313)
plot(t,y(:,3))

3 个评论
显示 1更早的评论 隐藏 1更早的评论

The @ character is missing when you define .
f = @(t,y) [...
Tesfaye Girma
Tesfaye Girma 2021-7-22
you are right thank you brother!!
Frane
Frane 2021-7-22
编辑：Frane 2021-7-22
Where did you get this numbers that are multiplying y(1), y(2) and y(3) in f=(t,y), and the numbers from y0?
Did you try my code and added these numbers and it worked or?
And what should y(1), y(2) and y(3) represent?

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息

标签

提问:

Frane
Frane
2021-7-19

编辑：

Frane
Frane
2021-7-22

Community Treasure Hunt

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

Start Hunting!

Translated by