System identification Windkessel - pulse wave
Show older comments
Hello, I would like to know if there is any way how to identify parameters if I have specified input to model and desired output. To be more specific I am modelling Windkessel model. It´s electrical model from two resistors one inductor and one capacitor. Wanted output is pulse wave.
Model:
And this is desired output of model (voltage at resistor R):
Can anyone help me please? Thank you
1 Comment
Rodrigo Vialle
on 30 Jul 2019
Hello, could you share your latest model?
Accepted Answer
Star Strider
on 4 Mar 2017
Edited: Star Strider
on 4 Mar 2017
If you provide the input and output as a ‘.mat’ file, and your model, it should be possible to derive the component values. (I can derive a circuit model for it if you want.)
You aren’t doing system identification since you already have a model for your system. You’re doing parameter estimation of your model.
EDIT —
The circuit model I derived for your model (I let the Symbolic Math Toolbox do the algebra) is:
syms C L R r Pi Po P2 s t
Z1 = 1/r;
Z2 = s*C + 1/R;
Z3 = 1/(s*L);
Zt = Z1 + Z2 + Z3;
I = Pi/Zt;
Po = Z2*I;
Out = ilaplace(partfrac(Po, s));
Out = simplify(Out, 'Steps',20);
Out_fcn = matlabFunction(Out, 'Vars',{t,Pi,r,C,R,L})
Out_fcn = @(t,Pi,r,C,R,L) Pi.*dirac(t)-(Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);
The ‘Z’ values are the three lumped impedances (and ‘Zt’, the total impedance), where the output is taken across ‘Z2’. The idea is straightforward: the current ‘I’ (or blood flow here) through the series of impedances is simply ‘Pi/Zt’. It is then straightforward to calculate the voltage (or pressure) drop across ‘Z2’ as ‘Z2*I’.
In order to fit your function and estimate the parameters, we need three data vectors: ‘t’, ‘Pi’, and ‘Po’, the time, input pressure with respect to time, and output pressure with respect to time. The function fits the output pressure it estimates with specific parameters to the ‘Po’ vector.
The procedure for estimating the parameters would then be:
% % % MAPPING: b(1) = r, b(2) = C, b(3) = R, b(4) = L
Out_fcn = @(t,Pi,r,C,R,L) Pi.*dirac(t)-(Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);
t = t(:); % Column Vector
Pi = Pi(:); % Column Vector
Po = Po(:); % Column Vector
b0 = ['Initial Parameter Estimates For r,C,R,L — IN THAT ORDER — As A Vector'];
nrmrsd = @(b) norm(Po - Out_fcn(t,Pi,b(1),b(2),b(3),b(4))); % Norm Of Residuals (Errors)
bvals = fminsearch(nrmrsd, b0);
The ‘bvals’ vector will be the estimated parameters.
NOTE — This is UNTESTED CODE since I do not have your data to test it with.
24 Comments
David Hatle
on 4 Mar 2017
Star Strider
on 4 Mar 2017
My pleasure!
I’ll put ‘r’ in parallel with ‘L’ and let the Symbolic Math Toolbox re-derive the model. I’ll post it here.
I can’t guarantee success with the parameter estimation, because actual physiological systems are notoriously difficult to identify.
I also need a time vector, so if you didn’t provide one (I’ve not loaded the 'mat files yet), I need to know the sampling interval (time in seconds between samples, assuming a constant sampling interval). I assume the lengths are the same for both ‘input’ and ‘output’, so the duration should be straightforward once I have the sampling interval.
I have done modeling with various physiological systems including Windkessel models (I’m a Board Certified Internal Medicine Physician and M.S. Biomedical Engineer). However I need for you to provide approximate initial values for the parameters ‘r,C,R,L’.
I will re-derive the model, and wait for you to provide the approximate initial parameter estimates. Without those, there could be any number of sets of parameter estimates that could fit, all but one of which would be completely meaningless to you, or that could result in a complete failure to fit your data. That’s simply the nature of nonlinear parameter estimation. So I very definitely need to know what the initial estimates are.
EDIT —
The re-derived model is:
Out_fcn = @(t,Pi,r,C,R,L) (Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);
All I now need are your initial parameter estimates and we’re cleared to go!
David Hatle
on 4 Mar 2017
Star Strider
on 4 Mar 2017
The time vectors both have to be the same. I can use the resample function (preferably) or interp1 (last resort for signal processing) to put them on the same time reference. They won’t match (there are phase and other considerations), but they have to have the same time base and length or MATLAB will throw a deservedly unrecoverable error.
If you know that the sampling interval is 1 sec, then I would truncate the longer sequence to the length of the shorter sequence. I will do so in my code.
Please note that 1.0s is a very long sampling interval in circulatory dynamics, since the pulse interval of humans is about 0.8s or less. I would choose a sampling interval of 0.5s, and preferably 0.1s, for circulatory dynamics studies.
I know it consumes memory and disk space. It also results in publishable data. Err on the side of publishable data. That’s why whatever Deity you believe in (or don’t) created DVD disks and 64 GB micro SD cards.
The initial parameter estimates are important, since I don’t know the system you’re studying. Values differ with respect to age and species. Yours will provide an acceptable initial estimate for the system you’re studying.
David Hatle
on 4 Mar 2017
Edited: David Hatle
on 4 Mar 2017
Star Strider
on 4 Mar 2017
I may use the Global Optimization Toolbox functions to find the best parameters. That could take a while.
I derived the ‘Out_fcn’ (in R2016b) from the information in your revised circuit diagram descriptions:
syms C L R r Pi Po P2 s t
Z1 = 1/(1/r + 1/(s*L)); % Parallel ‘rL’ Network
Z2 = 1/(s*C + 1/R); % Parallel ‘RC’ Network
Zt = Z1 + Z2; % Total Impedance
I = Pi/Zt; % Network Current
Po = Z2*I; % Output Pressure
Out = ilaplace(partfrac(Po, s)); % Inverse Laplace Transform
Out = simplify(Out, 'Steps',20);
Out_fcn = matlabFunction(Out, 'Vars',{t,Pi,r,C,R,L})
The rest of my code remains the same.
The radial artery information helps. I may go online with PubMed to see if anyone else has identified relevant parameters.
Star Strider
on 5 Mar 2017
I opted to interpolate the longer sequence to the shorter 201-element sequence. That way, I avoid creating data that aren’t actually present in the shorter sequence.
That aside, we have a problem with the data. Specifically, when I load your data, I notice some significant differences in the magnitudes of the source and output pressures:
inp = load('David Hatle input.mat');
otp = load('David Hatle output.mat');
Po = interp1(otp.t_output, otp.x_output, inp.time);
q1 = [inp.src(1:10)' Po(1:10)']
q1 =
0 8.9981e+03
2.1565e-02 8.9660e+03
4.2998e-02 8.9357e+03
6.4166e-02 8.9064e+03
8.4938e-02 8.8760e+03
1.0519e-01 8.8446e+03
1.2479e-01 8.8144e+03
1.4362e-01 8.7866e+03
1.6157e-01 8.7596e+03
1.7852e-01 8.7314e+03
The input pressures are in the first column, and the output pressures are in the second column. The network (that I believe I have derived correctly) is purely passive. It does not amplify.
From a purely technical perspective, another problem is that the source data begins at zero. That simply does not occur in normal physiology. The systemic arterial blood pressure limits are characteristically between about 80 mm Hg (Torr) or 10.5 kPa to 120 mm Hg or 16 kPa in a normal, healthy human. (The diastolic can be as low as 40 mm Hg or 5.3 kPa with aortic valve regurgitation, since the pressures are significantly lower in the pulmonary circulation, and left ventricular pressure has to drop below the pulmonary vein pressure in order to fill.) Aortic root pressure values in the range you have them are incompatible with life.
What are the units of the output pressures? Pascals?
I will rely on you to make the necessary corrections and scale the pressures appropriately. It is necessary to get the pressure units correct in order for the identified component values to have any meaning. Then we can try this again.
David Hatle
on 5 Mar 2017
Star Strider
on 5 Mar 2017
My code now estimates the parameters without returning a NaN value because of the initial zero value for the pressure. It still doesn’t give anything close to an appropriate fit. I even used patternsearch to see if it could find a better parameter set, but the result was essentially the same as with fminsearch.
I also tried shifting the aortic pressure waveform to approximately coincide with the radial artery waveform (see the Wikipedia article on the Cardiac cycle), again without success.
I used this code to do the shift:
idx = find(t <= 0.1, 1, 'last');
Pi = circshift(Pi, idx+1);
The problem may be your model. I went online to PubMed (using this search string) and found 18 articles that looked at Windkessel models in radial artery circulatory dynamics. I found one slightly informative (and free) article Noninvasive Pulse Wave Analysis for the Early Detection of Vascular Disease. Hypertension 1995;26:503-508. It’s not recent, but you can use its PMID: 7649589, as well as the search string, to search for others. (Most of my Windkessel experience involved ‘designing’ various pulse waveforms to be programmed into a pulsatile perfusion pump for cardiopulmonary bypass during my brief but spectacular foray into private industry long, long ago in a galaxy far, far away. My model belonged to my employer, so I wasn’t able to take it with me. It wasn’t published outside of our successful FDA Class III PMA application for the pump, and would probably be redacted anyway should anyone request it.) None of my Windkessel experience is recent.
I’ll be glad to work with you on this as my time permits. It’s certainly an interesting problem.
David Hatle
on 5 Mar 2017
Edited: David Hatle
on 5 Mar 2017
Star Strider
on 5 Mar 2017
My pleasure.
I generally do not give help other than here on MATLAB Answers. (The exceptions are only to resolve minor ambiguities that do not affect the substance of the Answers threads.) To do otherwise destroys the evolution of problem solutions, and creates discontinuities that others searching for help on the same (or similar) topics would need to solve their own problems.
My M.S. involved attempting to model the regulatory dynamics of carbohydrate metabolism. (I am Certified by the American Board of Internal Medicine, did a fellowship in Endocrinology.) That problem turned out to be much more complicated than I ever imagined, and is a problem that remains unsolved because of the complexity of the physiological and biochemical processes involved. I got ‘close enough’ to satisfy my committee. I am convinced that they awarded me my degree at least in part because they no longer wanted to hear my screams of frustration.
Physiological regulatory processes (one of which you are modeling) evolved over millions of years to be adaptive to myriad external disturbances in order to keep us alive, providing endless M.S. and Ph.D. research opportunities in the process (and endless frustration for clinicians like me when things go awry), since there is more to them than will ever be dreamt of in any of our philosophies.
That is the reason I find your research personally interesting. There is much to be discovered, and it is quite possible that you could find something new.
Physiological processes are in general nonlinear, nonstationary, and nontrivial. Enjoy the challenge!
David Hatle
on 5 Mar 2017
Star Strider
on 5 Mar 2017
I have the Toolbox, but I have no experience with it.
Simply reading the documentation, note that it first requires an accurate model, and that the model already be present as a Simulink model, so it seems you have to do that first. It then will allow you import the data and to tweak the parameters. You quite probably need to explore alternate models, preferably one that comes close to reproducing your radial pulse waveform.
Some Windkessel models that I have seen incorporate a rectifier or gated rectifier. It would likely be best for you first to explore the literature to see what others have done. That is where I always begin.
David Hatle
on 6 Mar 2017
Star Strider
on 6 Mar 2017
My pleasure.
As I remember, modeling the dicrotic notch requires a rectifier of some sort in order to simulate the closing of the aortic valve. I do not recall the details at present. (A rectifier with a bias can be a challenge to model. I do not remember now how I did it before.)
You can get some information from a 4-component model, as did the authors of the paper I referenced. It depends on what you want from the model.
David Hatle
on 6 Mar 2017
With some model adjustments I was able to find parameters for all components. And these are the normalised results:
Now I will focus on modelling dicrotic notch.
Star Strider
on 6 Mar 2017
The ‘model adjustments’ must make your new model significantly different from the model you first posted?
Your posted model could not produce this latest result.
David Hatle
on 6 Mar 2017
Star Strider
on 6 Mar 2017
Reconfiguring the model, adding 4 parameters, certainly are as well.
I assume you no longer need my help.
David Hatle
on 7 Mar 2017
Jan
on 7 Mar 2017
@David: Please press the "accept" button to accept the answer.
Star Strider
on 7 Mar 2017
@David — If you still need my help, I will provide it as much as I am able. I prefer to respond on MATLAB Answers to provide continuity. I simply felt ‘out of the loop’ with your new model (that you have not shared here).
@Jan — Thank you!
David Hatle
on 13 Mar 2017
Star Strider
on 14 Mar 2017
I am not aware of any way to convert units of liquid pressure, flow, and resistance to electrical units. (The exception being calibrated transducers where a transducer would be calibrated in terms of volt/torr or some such.) They are completely different physical properties. You would probably have to go back to the original definitions of the relevant units to see if you could determine direct relationships between them. (Blood has a critical Reynolds number of 2000 if I remember correctly, and while it generally stays below that value, blood is non-Newtonian and has nonlinear characteristics.)
I have never seen any conversions, but they may exist. I am not certain how to search the literature or the Internet for such information, although a physiology reference text could have them. A brief search found Units and Conversion Factors and List of equations in fluid mechanics. A more extensive search could discover what you want.
More Answers (1)
Hi David
1.
I have simulated the circuit of your question in ADS and I get the following transient
2.
tuning C1 and L1 both min values
3.
C2 and L1 both max values
4.
R2 value is quite irrelevant, the one that counts is R1/R2 it has to be kept small
5.
In MATLAB you can either build a model from ground zero, or could start with Simulink. Do you have Simscape?
So, David
if you find this answer useful would you please be so kind to mark my answer as Accepted Answer?
To any other reader, please if you find this answer of any help solving your question,
please click on the thumbs-up vote link,
thanks in advance
John BG
3 Comments
David Hatle
on 4 Mar 2017
John BG
on 5 Mar 2017
1.
I use MATLAB as often as possible, but for the task you asked to solve and for the time I allocated to your question, ADS save time, even if using Simscape.
2.
So, Simscapte it is, would you please attach to your question your Simscape circuit?
3.
Also, do you want to generate an ECG signal, or the wavelet transform of an ECG?
4.
Would the following in anyway help?
the Illustrated wavelet transform, by Paul Addison, CRC press, page 297
David Hatle
on 5 Mar 2017
Edited: David Hatle
on 5 Mar 2017
Categories
Find more on Transform Data in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
