Can anyone tell me how to implement these timevarying state space equations in matlab

Answer 1

[Ced]

I would use one of the ode functions, e.g. ode45. You can then define an arbitrarily complex, nonlinear, time-varying, etc. system.

In your case, something like that (parameters are rubbish, but you have those)

function [Tout, Xout] = run_ode
% Some Parameters
Rm = 1000;
Ra = 1000;

% Time-varying parts:
C = @(t)(sqrt(t+1)^3);
dC = @(t)(3*sqrt(t+1)/2);
r = @(x)x*(-x<0);   % Ramp function
p = @(x)[r(x(2)-x(1))/Rm ; r(x(1)-x(4))/Ra];

% Set up simulation
dx = @(t,x)ode_sys(t,x,C(t),dC(t),p(x));
Ts = 1e-12;
t_end = 1e-10;
t_sim = 0:Ts:t_end;

x0 = randn(5,1);

[ Tout, Xout ] = ode45(dx,t_sim,x0);

end

% Define system equations
function dy = ode_sys(t,x,C,dC,p)
% Some Parameters
Rs = 100;
Cr = 1e-6;
Cs = 2*Cr;
Ca = 0.1*Cr;
Ls = 6e-9;
Rc = 1e6;

b1 = 1/(Rs*Cr);
b2 = 1/(Rs*Cs);

% Defined time-varying system matrices elementwise
Ac = [ -dC/C 0 0 0 0 ; 
    0 -b1 b1 0 0 ; 
    0 b2 -b2 0 1/Cs;
    0 0 0 0 -1/Ca;
    0 0 -1/Ls 1/Ls -Rc/Ls ];

Pc = [ 1/C -1/C; -1/Cr 0; 0 0 ; 0 1/Ca ; 0 0 ];

dy = Ac*x + Pc*p;

end

Cheers

[Bilal]

Thanks alot for all this, so here you have selected the C = @(t)(sqrt(t+1)^3); randomly ? secondly i need to extract x(1) x(2) x(3) x(4) x(5)

[Bilal]

this is how i have calculated the time varying C(t) and dC(t), can you please tell me how to write C(t) which is the reciprocal of E(t) in ode45 syntax. C(t)=1/E(t)

clc;

close all;

clear all;

%%

Emax = 2;

Emin = 0.06;

HR = 60;

ts = 0.0025;

t=[0:ts:1];

tc = 60/HR;

Tmax = 0.2 + 0.15 *tc;

tn=t/Tmax;

En= 1.55*(((tn./0.7).^1.9)./(1+(tn./0.7).^1.9)).*((1./(1+(tn./1.17).^21.9)));

E= ((Emax - Emin) * En)+Emin;

plot(t,E);

grid on

%% for dC/dt

C=[1./E];

dcdt=diff(C)./ts;

figure

C=C(1:length(C)-1);

CC=dcdt./C;

plot(t(1:(length(t)-1)),-CC);

title('\Delta{C}/\Delta{t}')

grid on

[Ced]

Yes, I chose C(t) and dC(t) randomly.

One important thing is not quite clear to me. Do your functions C,dC, and E have to be discrete? or can you define them as continuous functions of time t?

[Ced]

Assuming you can define them as continuous functions: You don't need any particular syntax for C(t), dC(t) with ode45. The only function which needs a particular syntax is the function "dx" which I passed to ode45 itself. "dx" needs to have the variables (t,x). For the rest, you just need to define "standard" matlab functions.

Basically, there are two ways to define functions:

a) anonymous functions: "inline" functions as above with the @ sign

b) functions defined outside of the rest of your code. The syntax there would be "function output = function_name(input)". You can have them as separate .m files, or locally/nested in your main script file. All details are in the link, but this is not so important for the problem at hand.

You can find a good introduction with examples etc here: http://www.mathworks.com/help/matlab/functions.html

In your example, I would say that anonymous functions are enough, but I don't know the full complexity of your program. If there are a lot of interdependent and/or complex functions, I would define them as "local functions". See the link above for examples.

If you want to define E(t) as an anonymous function:

En = @(tn)1.55*(((tn./0.7).^1.9)./(1+(tn./0.7).^1.9)).*((1./(1+(tn./1.17).^21.9))); % note that this is a function, and not a vector of numbers
E = @(t) ((Emax - Emin) * En(t/Tmax))+Emin; % again, E(t) as a function

Let me know if something is unclear.

[Bilal]

E --> Elastance , C ----> Compliance , dC ----> Change in compliance all are continuous

[Bilal]

so do i need to give the range of tn ??

and what if i want to take the derivative of E here as i need C and dC which is C=1/E so from above equation E do i need to manually calculate the derivative can't i use the diff function any more as its a function

[Bilal]

from where it will get the value of tn ?

[Ced]

no, you don't need to give the range of tn for ode45. You just define E the same way you do for C, dC and the rest. When you then call e.g. E(3), then the function E(t) will be evaluated at t = 3;

The integration range is then passed to ode45. In my first answer, I had used

Ts = 1e-12;
t_end = 1e-10;
t_sim = 0:Ts:t_end;

x0 = randn(5,1);

[ Tout, Xout ] = ode45(dx,t_sim,x0);

t_sim is only passed to ode45, and does not appear anywhere else. Instead of defining t_sim as a time vector with all time steps, I could also only pass [ t_start t_end ], and then ode45 output it's own time steps. The reason I prefer passing t_sim as "full vector" is that it makes things easier to plot later, but that is a personal preference.

For the derivative: Yes, calculate the derivative by hand, and then hardcode the function. There are ways to compute derivatives symbolically using the symbolic math toolbox, but that would just make things unnecessarily complicated here. I think this derivative is still ok. If you want to check that your derivative is correct, you can always use some calculator or mathematica to check your answer. Also, note that "diff" does not compute derivatives as such, but only takes the difference between the elements in your vector, e.g. diff([ 1 2 3]) would return [ 1 1 ].

Answer 2

[Bilal]

yes you are right but the thing is E(t) is the re scaled version of En(tn) so it would definitely effect E(t) if i don't define tn anywhere if you run my code, you will see i m getting the write results as given in this paper. check the attachment

And yes if i accept the answer can we still continue communicating here. ?

[Ced]

Yes, I'm pretty sure we can keep adding to the answer once you accept it. Worst case, email me afterwards.

Thanks for the paper. But you do add Emin, Emax, and tc as known constants, correct? So tn is nothing else than t/Tmax, meaning that instead of writing E(t) = (Emax-Emin)*En(tn)+Emin, write it as a function of t, i.e.

E = @(t) ((Emax - Emin) * En(t/Tmax))+Emin;

Of course, if you really want, you could define

tn =@(t)t/Tmax

as well an use it like this, but I see no reason to do this. The only reason they introduce tn is that it makes the formula of En(tn) more compact for the paper.

*EDIT* Just realized this may not have come across clearly: We are only talking about ode45 here. If you want to plot the functions, you of course have to tell 'plot' at which points you want your functions evaluated. But let's maybe stick to solving the ode first, and then we can tackle the plotting part.

[Bilal]

i don't have any problem in writing t/Tmax or tn both are same all i wanted to know if i write express like this E = @(t) ((Emax - Emin) * En(t/Tmax))+Emin; in ode45 it means i need to have the values of all the contants except (t) as the function is dependent on t. Am i right, yes everything is known Emin Emax and tc.

[Ced]

yes, exactly.

[Bilal]

I have modified the code but now i am getting a very strange error in ode45. please have a look

[Bilal]

everything is as it is except the En(tn) and E(t) to find out C(t) and dC(t)

[Ced]

Yes, the error message is not so intuitive. The problem is that in your definitions of En(t), C(t), etc., which are functions of "t", you are using "tn". You did define "tn" as a function, which is fine, but this means that you have to pass an argument when you use it. In short: just replace "tn" by "tn(t)" in your definitions of En, E, C, etc., and you're good to go.

I'll have to pass for today, but I'll be back tomorrow.

Cheers

[Bilal]

Thanks alot for your effort and time it worked but u know what now Xout matrix is undefined having (NaN) in it

[Ced]

Yes, that's unfortunate :D. My best guess at the moment: Your initial conditions are still set as "randn" as far as I can tell. 'randn' takes random values from the normal distribution. If nothing is known, this is usually as good a guess as any other, but since you have a biological system, there may be some limitations to what makes sense.

Things you could do:

- Check all equations (in particular, signs)

- Check if the initial conditions are given somewhere in the paper.

- check if there are any upper/lower bounds on the values of x (e.g. a negative pressure makes little sense to me), make sure the initial conditions fulfil them

- There is most likely a division by zero somewhere. Try to find out where. To do this, you can e.g. write 'keyboard' in your code. The run will then pause there, and you can evaluate the code one line at a time by pressing f10. There are lots of other debugging tools for matlab, you can check the documentation for a list.

- your time scale seems to be similar to the ones in the paper, so that's good, and since the curve is very smooth and does not have high frequency noise or so, I don't think it's the solver tolerance. You could however try another solver (look at 'doc ode45' for what other versions are available), but since your solution is already NaN after 1 timestep, I don't think it's the solver.

[Bilal]

yes there was something wrong in my dC equation (the derivative of c equation) i choose dc = C just to check the NaN problem it worked but now all the values are same means nothing is updating, and there is no info in the paper regarding initial condition but if you see the plots all the values are positive so i made randn to rand but the values are small is that okay as its just the initial condition. one thing i wanna ask here is there any way to plot E(t) here, so that i can verify that i m getting the same plot as i was getting before and then i can plot C and dC if i m getting the same plots then I would shift to next step. Hope you can understand what i mean to say

[Ced]

Glad it works! Not sure I understand, no. What do you mean by "plot here"? Since E, C, and dC are only functions of time, they are not dependent on your solution from the ode, so I'm not sure what exactly you want to check.

But you can plot them whenever you like using "plot", e.g.

plot(t_sim,E(t_sim))

[Bilal]

yup i did it and i m getting the same plot

function [Tout, Xout] = run_ode

% Some Parameters

Rm = 0.0050;

Ra = 0.0010;

Emax = 2;

Emin = 0.06;

HR = 60;

tc = 60/HR;

Tmax = 0.2 + 0.15 *tc;

ts=0.0025

t=[0:ts:1];

% Time-varying parts:

tn = @(t) t./Tmax;

En = @(t) 1.55*(((tn(t)./0.7).^1.9)./(1+(tn(t)./0.7).^1.9)).*((1./(1+(tn(t)./1.17).^21.9)));

E = @(t) ((Emax - Emin) * En(tn(t)))+Emin;

plot(t,E(t))

[Bilal]

but i don't know how to take the derivative of this complex function E, i tried online solvers for derivative i am getting very wired expression from there.. will solve it then will get back to you.

[Ced]

Sounds good. Yes, solutions can look strange sometimes because they apply weird "simplifications". You could also simply derive it by hand, the first derivative is still fairly simple, and then compare the plot to the one from your online solution. Let me know if you need help with that.

[Bilal]

can you please check the problem in the code, getting some sort of error in ode_sys. Please have a look at the attached file.

[Ced]

I will have some time in a few hours, I'll let you know!

[Ced]

ok, had a look. A new version is in the attachments, but a couple of (hopefully) useful comments:

The error you are (were) getting is:

Error using vertcat Dimensions of matrices being concatenated are not consistent.

Error in tmp>ode_sys (line 60) Ac = [ -dC/C 0 0 0 0 ;

Error in tmp>@(t,x)ode_sys(t,x,C(t),dC(t),p(x))

[ ... ]

Error in tmp (line 45) [ Tout, Xout ] = ode45(dx,t_sim,x0);

Generally speaking, the first line tells you what kind of error it is, and then goes on telling you where this error was produced. In this case, the error is that "dimensions of matrices being concatenated are not consistent", and it happens in (line 60).

This is quite a common error which means that you are trying to concatenate two matrices/vectors (meaning "insert" them into the same array), but that there dimensions don't match. This usually happens when you messed up indices somewhere in a for loop or so. Unfortunately, in your case, there's a bigger problem. What you are doing with "diff" is that you are approximating the derivative by finite differences. But ode45 evaluates your function at a single timestep, and with only one value of C given, computing a finite difference is impossible.

I somewhat fixed your code with a hack which I think should work, but I would really encourage you to use the "real" derivative. I'll see if I find some time tomorrow and set it to you.

A couple of useful side notes:

1) The last index can be reach through "end", e.g.

x(1:end-1) % x without last element

which is a bit easier than your t(1:length(t)-1) construct :).

2) I wouldn't implement a function C referring to a previous version of C. I'm actually surprised this even worked. I simply omitted dC and plugged C directly into CC. If for some reason, you ever had to do that, I would use two different names, e.g. C0 and C1

Cheers

[Bilal]

Oh thanks alot it worked, now coming back to the state space, if you see in the paper, there are five state variables, x1 to x5. they r in dy in our code if im not getting it wrong, but i m getting only single value of x1 x2 x3 x4 x5. I need a vector so that i can plot these values.

[Ced]

Hi

Sorry it took me so long. I wanted to check our equations, but am running into big problems, nothing seems to match. Probably some stupid typo somewhere. Anyway... to your question: Tout is a vector of time, Xout is a matrix containing your states. So x1, x2, ... are in Xout. dy in our code is x_dot.

if your run the function from the matlab command line (you need to comment out "clear all" first in the function, otherwise it won't work), each xi will be in one column, i.e.:

[ t_out, x_out ] = run_ode;
x1 = x_out(:,1);
x2 = x_out(:,2);
etc.

[Bilal]

I know but x1 x2 ... n x5 never gets updated their values remain same all the time.

[Bilal]

although if you see C(t) and dC(t) is changing but the state variable never changes.

[Ced]

Ah, yes. Two things:

1. When computing the derivative with finite differences, you need to take the same "step size" on the top and bottom, e.g. for f([t t+dt1])/dt2, both dt1 and dt2 have to be the same (my mistake, sorry).

Replace the following lines in your code, this should solve the problem:

CC = @(t,ts) diff(C(t))./(C(t(1:end-1))*ts); % dC./C

[...]

plot(t(1:end-1),-CC(t,ts));

[...]

Ts = 0.0025;
dx = @(t,x)ode_sys(t,x,C(t),CC([t t+Ts],Ts),p(x)); % Reverse order!

2. Then, concerning the non-changing data: Have a look at your time step Ts and the end time. You are only simulating for 1e-10 seconds! During this time interval, nothing changes. You could e.g. set Ts = 0.0025 and t_end = 1 for one cycle, but that depends on what you want.

[Bilal]

wanted to plot the output of state variables for 3 cycles

[Bilal]

r = @(x)x*(-x<0); % Ramp function p = @(x)[r(x(2)-x(1))/Rm ; r(x(1)-x(4))/Ra];

how this statement will get the values of x(1) x(2) and x(3) it is written above the ode_sys call ... i m not getting it.

[Bilal]

I have tried to solve the problem, now its the update version of code pls run it and check that it only calculates states for 1 cycle. although now i have 4 or 5 cycles of time varying signal. so the plot in Xout (:,1) ... Xout(:,5) should repeact periodically like E(t) and C(t) ...

And attached herewith are the plots from paper that i have already sent you.

[Ced]

I'm sorry, I don't understand your question(s). Why don't you just plot it? That way, you'll see how many cycles appear? So... just to e.g.

plot(Tout,Xout(3,:))

And if it's not enough/too many cycles, adapt your end time accordingly (or simply only plot a particular section of the results).

[Bilal]

i have plotted them but i am getting only one cycle let me show u

Answer 3

[Bilal]

<<

>>

[Bilal]

If u see there the time range is okay i.e 2secs in Xout plots the problem is that output doesn't repeat itself after 0.4 secs

Answer 4

[Bilal]

[Bilal]

is there any way to check the values of r = @(x)x*(-x<0); % Ramp function

p = @(x)[r(x(2)-x(1))/Rm ; r(x(1)-x(4))/Ra];

p(x) -->

[Ced]

Yes, of course. Since you know Xout, you can also compute p, r, and everything else. If you want to have them as "standard outputs" like Xout, you could do something like :

Nt = length(Tout);
Pout = zeros(Nt,2);
Rout = r(Xout);
for i = 1:Nt
    Pout(i,:) = p(Xout(i,:));
end

after Tout and Xout have been computed and add them to the output, i.e.

function [Tout, Xout, Pout, Rout] = run_odebq

IMPORTANT: You will need to change your ramp function to an element-wise evaluation, otherwise you get some errors. Meaning it should look like this (note the '.'):

r = @(x)x.*(-x<0);   % Ramp function

Once you have Pout, Rout, you can plot them, do some analysis or whatever (you can of course also plot them directly in the function without passing them as outputs)

[Bilal]

okay, i got it thanks alot.

Did you find anything why i am getting only one cycle in Xout and then output gradually goes to zero

[Ced]

No, I'm afraid I'm a bit caught up in work myself at the moment. I'll have a look at it if I find the time. I haven't read the paper in detail, but maybe they have something like mod(t,HR) or so? Maybe you can try to find out what function is not periodic, and if not, why.

Answer 5

[Bilal]

This is what exactly i have done in the generation of C(t) i have use mod(t,HR) but do i need to implement it for the state space as well ? okay please get back to me when u have time. thanks

Answer 6

[Ced]

I mean, that would certainly work. simply replace every instance of "t" by "mod(t,60/HR)" when calling ode_sys, so

dx = @(t,x)ode_sys(mod(t,60/HR),x,C(mod(t,60/HR)), ...

The question is if this makes sense. That's something you could ask your professor or discuss with some phd student. Skimming the paper, I cannot find any reference that they artificially periodicized (that's a word now ^^) their solutions. But then again, little details are often omitted (or forgotten) when describing simulations in papers.

[Bilal]

Thanks alot by doing this the periodicity issue has been resolved, I am getting the shape of all the states exactly like the results on paper (Except X3).

[Bilal]

r = @(x)x*(-x<0); % Ramp function

Can you please tell me how this is interpreted.

we need output equal x as long as x>0 so how the above statement works

[Bilal]

Secondly if you see my plots they change on every cycle, whereas in paper there is no change, the plot is periodic.

[Bilal]

I have tried with fixed initial condition but it doesn't work either.

[Ced]

Great!

What you want the ramp function r(x) to do is the following: If the value x is greater or equal than zero, set r(x) = x, else set r(x) = 0. The

(-x<0)

part returns a logical value, which is 1 if (0 <= x) and zero else (if used in an algebraic expression, the logical is automatically interpreted as a double). I wrote it as -x<0 to have an inequality instead of a <=.

So, x*(-x<0) is x*1 if 0 <= x and 0 else, which is just what the ramp function does.

Now, to the changing period: This is a little strange. To be fair, the first period in the paper is also slightly different than the others. Check if the mod(...) arguments are the same everywhere in your code. Unfortunately, it could be that the equations are very sensitive to changes in the derivative, and since our dC is approximated, that could explain the differences. You could test this by changing the magnitude of the step in the finite difference computation (the second argument of dC). If that does not change anything, let me know, and I'll think a bit harder :)

[Bilal]

thanks for the explanation let me try this

[Bilal]

yes mod argument is consistent throughout the code. this is what i got after increasing the magnitude by the factor of 2.

[Bilal]

if true
  % code
end

[Ced]

hm... looks fishy to me, but I can't really put my finger on the reason. One thing I would certainly change is that you should not use randn for your initial conditions, since your pressure should not be negative.

And then, although I don't think that the derivative is really the problem, there is a strong dependency on it, as the solution changes when I change the discretization. That should not happen. If you could find the algebraic derivative, that would certainly help a lot. If you can't figure out the correct derivative by hand, you could try using the symbolic math toolbox, e.g.

syms t_smb real              % define symbolic variables

f = t_smb.^2 + 3;            % define symbolic expression (would be your function)

df = diff(f,t_smb);   % compute derivative wrt t_smb

df_fun = matlabFunction(df); % transform into matlab function

Note that here, "diff" is overloaded and computes the algebraic derivative and not the finite difference, since your input is symbolic.

Also, I actually find it a bit weird to have mod(t,60/HR). It would make more sense to have mod(t,Tmax) since t is "normalized" by Tmax. However, this diverges. Do you have anyone you can ask (supervisor or so) about this periodicity and normalization who is working in that field? While I know how ode works, I don't really know anything about heart simulations.

[Bilal]

yes i will ask my supervisor. getting too wired expression through this method pls check did i write the expression correctly ?

if true
  % code
end

[Bilal]

there?

[Bilal]

I am unable to plot derivative of C that you told me to generate using syms. its getting mpower error

Answer 7

[Bilal]

Please reply me back as soon as you get free. Thanks