Code

1
Code

• load wave1.txt
• wave wave1
• N length(wave)
• time 0N-1/200
• figure(1)clf
• plot(time,wave(,2))
• set(gca,'Fontsize',16)
• grid
• xlabel('time sec')
• ylabel('pressure mmHg')
• print -depsc2 fig1.eps

2
Data
3
Code

• idx1 max(find(timelt175))
• idx2 max(find(timelt240))
• figure(2)clf
• plot(time(idx1idx2),wave(idx1idx2,2))
• grid
• set(gca,'Fontsize',16)
• xlabel('time sec')
• ylabel('pressure mmHg')
• xlim(175 240)
• print -depsc2 fig2.eps

4
(No Transcript)
5
Code

• idx1 max(find(timelt175))
• idx3 max(find(timelt190))
• figure(3)clf
• hplot(time(idx1idx3),wave(idx1idx3,2))
• set(h,'Linewidth',2)
• grid
• set(gca,'Fontsize',16)
• xlabel('time sec')
• ylabel('pressure mmHg')
• xlim(175 190)
• print -depsc2 fig3.eps

6
(No Transcript)
7
Sensitivities

• CV model
• Model output
• Sensitivities

7
NC STATE University
8
Computation of sensitivities

• From the model using the chain rule for
  differentiation
• with initial conditions
  
• Remarks The derivatives can be computed
  by hand
• Automatic Differentiation can be used to compute
  sensitivities

NC STATE University
9
Computation of sensitivities

• Finite differences
• where

NC STATE University
10
CV model

• Model output
• Parameters
• Sensitivities

11
CV model

• Relative sensitivities
• Ranked sensitivities (scaled 2-norm)

12
Model
13
Model
14
Experimental data
15
Sensitivities
Es
16
Ranked Sensitivities
17
Sensitivities.m - code

• x0 load_global
• sols,pas cbf_resp(x0)
• h0.01 (something small sqrt(odetol)10)
• for i 1NP
• soli,pasi cbf_resp(x0dh)
• s(,i) (pasi-pas)/h
• sr(,i) s(,i)x0/pas
• sn(i) norm(sr/N,2)
• end

18
Correlated parameters

• Let x be given by
• The solution of x is
• Let the model output
• Parameter identifiability It is clear that
  parameters and can be identified, while
  and are correlated.
• But both parameters will be sensitive!

19
Subset selection

• Compute Jacobian and SVD
• Determine the numerical rank of
• Partition the matrix of eigenvectors
• Determine a permutation matrix P by constructing
  a QR decomposition with column pivoting for
• That is determine P such that
• Use P to reorder the parameter vector as
• Make the partition
• Fix at the a priori estimates
• Compute the new estimate of the parameter vector
• by solving the minimization problem.

20
Optimization

• Minimize the difference between computed and
  measured values of blood pressure
• Use Matlabs fminsearch

21
Fminsearch

• x fminsearch(fun,x0)
• x0 initial parameters
• fun function to be minimized
• Algorithm fminsearch uses the simplex search
  method. This is a direct search method that does
  not use numerical or analytic gradients.

22
The simplex method

• It uses the concept of a simplex, which is a
  polytope of
• vertices in dimension (a line segment
  in one-parameter space, a triangle in
  two-parameter, a tetrahedron in three-parameter
  space, etc.)
• Consider a 2-parameter estimation
  , simplex is a triangle
• Evaluate and sort them so that
  
• Idea Replace the worst point with a point with
  a lower cost value!

23
The simplex method

• Compute the centroid of the simplex
  (not including the worst point)
• Replace with
• If none of these values are better than the
  previous worst, algorithm shrinks the simplex
  towards the best point

24
The simplex method

• Remarks
• In theory, the simplex algorithm is not
  guaranteed to converge to a minimum and can
  stagnate at a suboptimal point. However, in
  practice, the method performs well and often
  results in an initial rapid decrease of the
  objective function value.

25
Experimental data (note HR not constant)
26
Periods
27
Periods
Period Time (sec) T (sec) HR 60/T (BPM)
1 0 0.64 93.75
2 0.64 0.64 87.59
3 1.325 0.685 95.24
4 1.955 0.63 90.23
5 2.62 0.665 90.91
6 3.28 0.66 86.96
7 3.97 0.69 86.96
8 4.66 0.69 88.89
9 5.335 0.675 94.49
28
Heart rate
29
Optimized blood pressure
30
Optimization code

• x0 load_global
• load data.mat (load td, pd, periods)
• x0 Rp Rs Cap Cvp Ca Cv Ts Tr Edl Esl Edr Esr
  Vdl Vdr

31
Optimization code

• x fv exit out t sol pas optimize
• x0 load_global
• x0 sqrt(x0)
• opt optimset(MaxFunEvals,1e8,TolX,1e-8,MaxI
  ter,5000)
• x fv exit out fminsearch(_at_cbf_resp,x0,opt)
• xx2
• t,solode15s(_at_cvmodel,..)
• Plot(td,pas,td,pd)

32
Optimization code

• J cbf_resp(x)
• global td, pd
• x x2
• t, solode15s(_at_cvmodel,td,INIT,options,x)
• passol(,3)
• J sum(abs(pd-pas)2)

33
To run an optimization

• Type
• x, fv, exit, out, t, sol, pas optimize