MATLAB: System identification Windkessel – pulse wave

biomedicalmodelmodellingparameterssystem identificationwindkessel

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

Best Answer

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.
Related Question