Why does this ode23 simulation not run/plot

I am trying to do something which is slightly above my paygrade as a MATLAB noob, but my code is as follows:

first function:

function v=pwm(t,T,d)
v=zeros(1,numel(t)); 
for i=1:numel(t)
k = floor(t(i)./T);
if ((k*T<=t(i)) & (t(i)<(k+d)*T))
    v(i)=1;
end
if (((k+d)*T<=t(i)) & (t(i)<(k+1)*T)) %#ok<*AND2>
    v(i)=0;
end
end
end 

second function (differential equation):

function dx=boost(t, x, v_in, R, C, L, z) %#ok<INUSL>
dx=zeros(2,1);
for i=1:numel(z)
if z(i)==1
    dx(1) = v_in./L ;
    dx(2) = -x(2)./R*C ;
end
if z(i)==0 & x(1)>=0
    dx(1) = (v_in - x(2))./L ;
    dx(2) = x(1)./C - x(2)./R*C ;
end
if z(i)==0 & x(1)<0           %#ok<*AND2>
    dx(1) = 0 ;
    dx(2) = -x(2)./R*C ;
end
end

differential equation solving/simulation :

v_in = 1.5;
R = 680;
C = 200e-6;
L = 180e-6;
T = 1e-3; 
d = 0.1;
t = linspace(0,0.2,2e6);
z = pwm(t,T,d); 

tspan = (0:1e-7:0.2);
x0 = [0 0];
options=odeset('MaxStep', 1e-5);
[tv,xm] = ode23(@(t,x) boost(t, x, v_in, R, C, L, z), tspan, x0, options);
figure(1)
plot(tv,xm)
grid
legend('i_L', 'V_O', 'PWM')

My goal with this code is for the function 'pwm' to have a time vector input, and output which is a vector containing a series of 1's and 0's. This is my pwm regulated control signal, but all it needs to be thought of is as a vector of 1's and 0's. We then have the boost function which defines 3 sets of differential equations as seen by the three if statements. Which if statement, and by extension which set of differential equations is used it based on z(i), (and in some cases x(1)). Moving onto the simulation code, we can see that I have assigned z to the output vector produced by the pwm function which has an input time vector defined by t = linspace(0,0.2,2e6). Finally, I attempt to actually solve the system of differential equations using ode23, and then plot the result.

when it comes to running the simulation it pretty much never concludes, and I never end up seeing a result - any help as to why would be great

Answer 1

[Star Strider]

Give it time!

I shortened ‘tspan’ to 150 elements to see if there was a problem with the code (there isn’t):

tspan = linspace(min(tspan), max(tspan), 150);

and it took 585 seconds, about 4 seconds per element. At that rate, with ‘tspan’ of elements, the integration should finish in 90.28 days.

The problem may be that you are passing ‘z’ and evaluating every element of it (it is the same length as ‘t’, and ‘tspan’ however it is evaluated at different time points than ‘tspan’) in a for loop, rather than calling ‘pwm’ as a function as in the earlier code, and with the time points corresponding to ‘tspan’.

So, it might be best to shorten ‘tspan’ withoiut affecting it limits. One way is to re-define it just after the original ‘tspan’ assignment (as I did) using linspace, choosing an appropriate number of elements for the new ‘tspan’.

Also, it is llikely best to use yyaxis for the plots, since ‘xm(:,1)’ and ‘xm(:,2)’ have significantly different magnitudes:

figure(1)
yyaxis left
plot(tv,xm(:,1))
yyaxis right
plot(tv,xm(:,2))
grid
legend('i_L', 'V_O', 'PWM')

You will only be pltting two values, so the third (’PWM’) will not appear. Plotting it will require a separate plot call, then plotting it as a function of whatever time argument you choose to evaluate it with.

[Jonah Foley]

ok, glad it was coded correctly, but not useful if its gonna take 90 days haha! Honestly, before I tried this method I had the following code, which did not use this additional variable z :

differential equation:

function dx=boost(t, x, v_in, R, C, L, T, d)
dx=zeros(2,1);
if pwm(t,T,d)==1 %evalutates single element of pwm or whole output vector?
    dx(1) = v_in./L ;
    dx(2) = -x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)>=0
    dx(1) = (v_in - x(2))./L ;
    dx(2) = x(1)./C - x(2)./R*C ;
end
if pwm(t,T,d)==0 & x(1)<0           %#ok<*AND2>
    dx(1) = 0 ;
    dx(2) = -x(2)./R*C ;
end
end

simulation:

v_in = 1.5;
R = 680;
C = 200e-6;
L = 180e-6;
T = 1e-3; 
d = 0.5;
t = linspace(0,0.2,2e6);
tspan = linspace(0,0.2,2e6); 
x0 = [0 0];

z = pwm(t,T,d); %necessary?

options=odeset('MaxStep', 1e-5);
[tv,xm] = ode23(@(t,x) boost(t, x, v_in, R, C, L, T, d), tspan, x0, options);
figure(1)
plot(tv, xm(:,1))
hold on
plot(tv, xm(:,2))
hold on
plot(t,z)
grid
legend('i_L', 'V_O', 'PWM')

the pwm function remained the same

Now, this code worked, and produced graphs for x(1) and x(2), but i wasnt sure if they were the correct graphs I was expecting, which lead me to the second method. the reason I thought this wasnt working correctly was because in the 'boost' function I use

if pwm(t,T,d) == 1

and by my understanding this would evaulate the whole vector poduced by the pwm function and check if it had a value of 1, whereas what I really wanted to do was have the pwm signal of 1's and 0's running along the same time scale as the boost function, then evalutate at each division in the time scale whether the element in the pwm vector was a 1 or 0, and then from there chose the correct if statement for ode23 to integrate and simulate. I do hope I'm making sense with this.

First, could you confirm whether the above code would do as I originally intended, if its clear what that was. And second, considering the boost function features those if statements which utilise the pwm function, is it necessary to call the pwm function in the simulation in order to generate that vector of outputs, or when those if statements are evaluate is the function pwm going to be executed within the boost function anyway.

Thank you for the help so far.

[Star Strider]

My pleasure.

That looks like the previous code, other than the value for ‘d’. (If it isn’t, please mention what changed.)

A couple things I just now noticed that could be problems:

First, this assignment:

-x(2)./R*C

is equivalent to:

-x(2)*C/R

so if you want to evaluate the product of ‘R’ and ‘C’ instead, do this:

-x(2)./(R*C)

and so for the others.

Second, this test:

if pwm(t,T,d) == 1

will be true (or 1) only if the match is exact. It may be appropriate to use ranges (such as ‘((pwm(t,T,d) >= 0.99) & (pwm(t,T,d) <= 1.01))’) if there is reason to suspect that is a problem.

I will help as I can!

[Jonah Foley]

ok so the difference is that in the initial code i showed you, in the simulation part I did

z = pwm(t,T,d)

which could produce the vector z with a range of values of 0 and 1

then in the boost function the if statements were evaluated as so

if z(i) == 1 .... // do something

which means I was evaluating each element in the z vector each time. However in the latest version I showed you notice that in the if statements I use

if pwm(t,T,d) == 1 ... // do something

first things first, the output of the pwm function is only ever exactly 1 or exactly 0, so I do not suspect that is the problem but I take you point. what I was mainly wondering is what w ould the difference be between the two methods described above. i.e with the statement ' if pwm(t,T,d) == 1 ' what exactly am I doing with this ... am I evaluating the whole vector (somehow) or am I evaluating the element of the vector which corresponds to the current time in the simulation, so that as the time progresses in the simulation I work through the full vector of 1's and 0's (which is what I want) ?

[Star Strider]

In this version, since ‘t’ and ‘tspan’ now appear to be identical, there’s no difference. It’s likely best to evaluate ‘pwm’ for various values of ‘t’ in the ‘boost’ function. I doubt that any efficiency is gained by pre-calculating it and then comparing it to the instant value of ‘t’ in a specific iteration. Evaluating in each iteration (rather than the vector comparison) is much more efficient, and is likely to produce the correct result of the logical comparison.

‘... am I evaluating the whole vector ...’

In the previous (originally posted version here), yes, however without the benefit of using either the any or all functions, so it is likely not possible to determine exactly what the result is.

‘... or am I evaluating the element of the vector which corresponds to the current time in the simulation ...’

If call ‘pwm’ with the value of ‘t’ of the instant iteration, yes, otherwise, no.

The previous (previous Question) code is correct with respect to the output and use of ‘pwm’, at least as I understand what you’re doing (no promises, there).

[Jonah Foley]

ok, this was really useful, thank you

[Star Strider]

As always, my pleasure!