Why ode45 is not working??

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.

[Torsten]

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]

I've uploaded the codes in text files and the error message.

[Torsten]

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]

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]

numel: number of array elements (link)

[Torsten]

Again look at the necessary changes to be made I indicated in my response.

You didn't incorporate them properly in your code.

[Frane]

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]

Again check the call to ode45 with my response. Note that {} and () are different in Matlab.

[Frane]

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]

What is the error message ?

[Frane]

[Torsten]

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]

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]

You overwrote the scalar value for K by K = lqr(A,B,Q,R). Use another name here.

[Frane]

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]

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]

I get the following error message.

[Torsten]

plot (t,squeeze(X(1,:,:)))

[Frane]

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]

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]

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]

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]

help xlim

help ylim

[Frane]

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]

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]

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]

And why do these quantities depend on time ?

[Frane]

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]

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]

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]

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]

Sorry, but I have no experience with the output of function "stepplot".

You should open a new question.

[Frane]

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]

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]

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.

Answer 1

[Tesfaye Girma]

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))

[Ildeberto de los Santos Ruiz]

The @ character is missing when you define .

f = @(t,y) [...

[Tesfaye Girma]

you are right thank you brother!!

[Frane]

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?