King M.R., Mody N.A. Numerical and Statistical Methods for Bioengineering: Applications in MATLAB

Подождите немного. Документ загружается.

The Lagrangian function for a minimization problem subjected to m equality con-

straints and p inequality constraints is

L x; lðÞ¼f xðÞþ

j¼1

ðxÞþ

k¼1

mþk

xðÞþσ



: (8:16)

Equation (8.16) is the augmented objective function in n + m +2p v ariables with no

constraints. At the optimum point,



ðÞ¼0;



ðÞþσ

¼ 0;

and

L x



; λðÞ¼f x



ðÞ:

A solution is a feasible critical point if it satisﬁes the following conditions:

∂Lðx; λÞ

∂x

¼ 0; i ¼ 1; 2; ...; n;

∂Lðx; λÞ

∂l

¼ 0; j ¼ 1; 2; ...; m þ p;

∂Lðx; λÞ

∂σ

¼ 2l

mþk

¼ 0; k ¼ 1; 2; ...; p :

Note that when σ

6¼ 0, the kth inequality constraint is inactive and l

mþk

¼ 0, since

this inequality constraint does not affect the search for the optimum.

Example 8.4

Minimize the function

; x

ðÞ¼x

þ 1ðÞ

x

; (8:17)

subject to the constraint

; x

ðÞ¼x

 x

 4:

The Lagrangian function for this constrained problem is given by

; x

; l; σ



¼ x

þ 1ðÞ

 x

þ l x

 x

 4 þ σ



A critical point of L satisfies the following conditions:

∂L

∂x

¼ 2 x

þ 1ðÞx

 l ¼ 0;

∂L

∂x

¼x

 l ¼ 0;

∂L

∂l

¼x

 x

 4 þ σ



¼ 0;

∂L

∂σ

¼ 2lσ ¼ 0 :

The last equation states that if the inequality constraint is inactive (σ

0Þthen l ¼ 0, and if the inequality

constraint is active (σ

¼ 0Þthen l 6¼ 0. In other words, both σ and λ cannot be non-zero simultaneously.

When l ¼ 0, only the first two conditions apply. In this case the critical point is given by

¼ 0; x

¼ 2:

527

8.4 Constrained nonlinear optimization

This is a saddle point, as is evident from the contour plot shown in Figure 8.21. The Hessian of the objective

function at point (0, 2) is given by

H ¼

2 1

10



;

and H

¼1

0: Therefore, the Hessian is not positive definite at this critical point, and (0, 2) is

not a local minimum. Note that the Hessian is also not negative definite, indicating that this point is also

not a local maximum.

When l 6¼ 0, the first three conditions reduce to a system of three simultaneous linear equations whose

solution is

¼

; x

¼

; l ¼

Upon inspection of Figure 8.21, we find that this is indeed a local minimum of the constrained function.

Another technique used to solve constrained optimization problems is the penalty

function method, in which a penalty function is constructed. The penalty function is a

parameterized objective function that is not subject to any constraints. The penalty

function can be minimized using any of the techniques discussed in Section 8.3. If the

constrained optimization problem is to minimize fðxÞ, subject to m equality con-

straints h

xðÞ¼0; j ¼ 1; 2; ...; m, then the penalty function is deﬁned as

P xðÞ¼f xðÞþr

j¼1

xðÞ



; (8:18)

Figure 8.21

Contour plot of the objective function with the inequality constraint shown by the dashed line. The feasible region is

located above this line

−6.11538

−4.23077

−2.34615

−0.461538

1.42308

3.30769

5.19231

7.07692

8.96154

10.8462

12.7308

14.6154

16.5

6.5

18.3846

20.2692

22.1538

24.0385

25.9231

27.8077

29.6923

31.5769

33.4615

−4

−2

−4

−3

−2

−1

The penalty function can be constructed in several different ways.

528

Nonlinear model regression and optimization

where the parameter r is a positive number. The penalty function method is an

iterative minimization scheme in which the value of r is gradually increased at each

iteration, and Equation (8.18) is successively minimized to obtain a series of optimal

points. The penalty function is ﬁrst minimized using a small value of r, say 1.0. The

optimum value of x obtained at this iteration is used as the starting guess value for

the next iteration, in which the value of r is increased. The increase in r with each

subsequent iteration shifts the optimum point of Equation (8.18) towards the

optimum x*off xðÞ.Whenr ! ∞, k hxðÞk!0, and the optimum of the penalty

function approaches the optimum of the objective function. Inequality constraints

can also be incorporated into a penalty function through the use of slack variables,

as demonstrated earlier in the Lagrange multiplier method.

Using MATLAB

The fmincon function in Optimization Toolbox performs constrained minimiza-

tion of a function of several variables. The constraints can be of the following types:

(1) linear inequality and equality constraints: Ax  b; A

x ¼ b

(2) nonlinear inequality and equality constraints: Cx  0; C

x ¼ 0,

(3) bounds: l ≤ x ≤ u, where l is the lower bound vector and u is the upper bound vector.

The syntax for fmincon is

x = fmincon(func, x0, A, b, Aeq, beq, lb, ub, confunc, options)

func returns the value of the objective function at any point x; x0 is the starting guess

point, A is the linear inequality coefﬁcient matrix, b is the linear inequality constants

vector, Aeq is the linear equality coefﬁcient matrix, and beq is the linear equality

constants vector. If there are no linear inequality constraints, assign a null vector [] to

A and b. Similarly, if no linear equality constraints exist, set Aeq = [] and beq=[].

Note that lb and ub are vectors that deﬁne the upper and lower bounds of the

variables. If no bounds exist, use a null vector as a place holder for lb and ub.

Also, confunc is a function that deﬁnes the nonlinear inequality and equality con-

straints. It calculates the values of Cx and C

x. The objective function func is

minimized such that Cx  0andC

x ¼ 0. If nonlinear constraints do not exist for

the problem, then the non-existent function should be represented by a null vector [].

There are several different ways in which the functi on can be called. See help

fmincon for more infor mation on its syntax.

Example 8.5

Use fmincon to solve the constrained optimization problem described in Example 8.4.

We create two functions m-files, one to evaluate the objective function, and the other to evaluate the

nonlinear constraint.

MATLAB program 8.20

function f = func(x)

f = (x(1) + 1)^2 - x(1)*x(2);

In MATLAB, we type the following. The output is also displayed below.

A = [-1 -1]; b = 4;

fmincon (‘func’, [0,1], A, b)

529

8.4 Constrained nonlinear optimization

Warning: Trust-region-reﬂective method does not currently solve this

type

of problem, using active-set (line search) instead.

> In fmincon at 439

Local minimum found that satisﬁes the constraints.

Optimization completed because the objective function is non-

decreasing in

feasible directions, to within the default value of the function

tolerance,

and constraints were satisﬁed to within the default value of the

constraint tolerance.

Active inequalities (to within options.Tolcon = le-006) :

lower upper ineqlin ineqnonlin

ans =

-1.5000 -2.5000”

What happens if you provide the initial guess point as (0, 2), which is the saddle point of the objective

function?

8.5 Nonlinear error analysis

Nonlinear regression gives us best-ﬁt values of the parameters of the nonlinear

model. It is very important that we also assess our conﬁdence in the computed

values and future model predictions. The conﬁdence interval for each parameter

conveys how precise our estimate of the true value of the parameter is. If the

conﬁdence intervals are very wide, then the ﬁtted model is of little use. The true

parametric values could lie anywhere within large conﬁdence intervals, giving us

little conﬁdence in the best-ﬁt values of our model. In Chapter 3, we demonstrated

how to obtain the covariance matr ix Σ

of the linear model parameters x:

¼ E x  μ

ðÞx  μ

ðÞ



Along the diagonal of the covariance matrix are the variances associated with each

best-ﬁt parame ter, which are used to estimate the conﬁdence intervals.

Two popular techniques for estimating the error associated with the ﬁtted

nonlinear model parameters are the undersampling method and the bootstrap

method. In both methods, several data sets are generated by reproducing values

from the original data set. The mathematical model is then ﬁtted to each of the data

sets derived from the master data set. This produces several sets of parameters and

therefore a distribution of values for each parameter. The distributions are used to

calculate the variances and covariances for each pair of parameters.

(1) The method of undersampling This method is very useful when a large data set is

available. Suppose N is the number of data points in the data set. If every mth point

in the data set is picked and used to create a new, smaller data set, then this new data

set will contain N/m points, and m distinct data sets can be derived this way. In this

method each data point in the original data set is used only once to create the m

derived sets. The model is ﬁtted to each of the m data sets or “samples” to produce m

estimates of the model parameters. The standard deviation of each parameter can be

530

Nonlinear model regression and optimization

calculated from their respective distributions. Alternatively, the covariance matrix

can be estimated using the following formula:

¼ s

m  1

 x



: (8:19)

Equation (8.19) is the estimated covariance between the ith and jth model parameter,

where i; j ¼ 1; 2; ...; n. Two limitations of this method are (i) N must be large and (ii)

the nonlinear regression must be performed m times, which can be a time-consuming

exercise.

(2) The bootstrap method This method is appropriate for large as well as moderately

sized data sets. The N data points in the data set serve as a representative sample of

the inﬁnite number of observable data points. Another data set of size N can be

created by performing random sampling with replacement (via a Monte Carlo

simulation) from the N data points of the original data set. In this manner, m data

sets are constructed from the original data set. Nonlinear regression is then per-

formed m times (once for each newly constructed data set) to obtain m sets of the

parameter values. The covariance between each pair of parameters is calculated using

Equation (8.19). You may ask, How many N sized data sets should be produced from

the original data? Note that m should be large enough to reduce sampling bias from

the bootstrap procedure, which can inﬂuence the covariance estimate.

Box 8.4B Kidney functioning in human leptin metabolism

Estimate the variance in the model parameters K

and R

max

using the bootstrap method.

The data set has N = 16 points. We generate m = 10 data sets each with 16 data points, by

replicating the data points in the original data set. Program 8.22 performs the bootstrap routine and calls

fminunc to perform a nonlinear regression to obtain best-fit parameter values for all ten data sets.

Note the use of the MATLAB function cov(X) in Program 8.21. This function calculates the

covariance between the n variables that are represented by the n columns in the m × n matrix X,

where each row is a different observation.

A nonlinear optimization algorithm may not always converge to the correct optimal solution. Be

aware that, occasionally, automated optimization procedures can produce strange results. This is why

good starting guesses as well as good choices for the options available in minimizing a function are

critical to the success of the technique. We present here the results of a successful run. The ten sets of

best-fit values for the two parameters of the model were obtained as follows:

pdistribution =

4.4690 29.7062

7.1541 70.2818

4.3558 32.4130

3.2768 18.2124

4.4772 21.9166

8.1250 50.4744

3.2075 26.1971

4.1410 22.5261

7.0033 60.6923

3.5784 21.1609

The mean of each distribution is

max

¼ 4:979;

¼ 35:358:

531

8.5 Nonlinear er ror analysis

The covariance matrix is computed as

Σ ¼

3:147 29:146

29:146 339:49



From the covariance matrix we obtain the standard deviations of the distribution of each parameter as

max

¼ 1:774

and

¼ 18:42:

MATLAB program 8.21

% Bootstrap routine is performed to estimate the covariance matrix

clear all

global S

global R

% Data

% Plasma leptin concentration

Sorig = [0.75; 1.18; 1.47; 1.61; 1.64; 5.26; 5.88; 6.25; 8.33; ...

10.0;11.11; 20.0; 21.74; 25.0; 27.77; 35.71];

% Renal Leptin Uptake

Rorig = [0.11; 0.204; 0.22; 0.143; 0.35; 0.48; 0.37; 0.48; 0.83; ...

1.25; 0.56; 3.33; 2.5; 2.0; 1.81; 1.67];

% Other variables

N = length(Sorig); % size of original data set

m = 10; % number of derived data sets of size N

% Bootstrap to generate m data sets

dataset(1:N,2,m) = 0;

for i = 1:m

for j = 1:N

randomno = ﬂoor(rand()*N + 1);

dataset(j,:,i) = [Sorig(randomno) Rorig(randomno)];

end

% Perform m nonlinear regressions

pdistribution(1:m, 1:2) = 0; % preallocating m sets of parameter

values

for i = 1:m

S = dataset(:,1,i);

R = dataset(:,2,i);

p0 = [1.732, 10.87]; % initial guess value obtained from Box 3.9A

p1 = fminunc(‘SSEleptinuptake’, p0);

pdistribution(i,:) = p1;

end

% Calculate the covariance matrix

sigma = cov(pdistribution);

muRmax = mean(pdistribution(:,1));

muKm = mean(pdistribution(:,2));

532

Nonlinear model regression and optimization

8.6 End of Chapter 8: key points to consider

(1) The optimal value of the parameters of a nonlinear model are obtaine d by perform-

ing nonlinear least-squares regression or minimization of

SSE ¼

i¼1



ðÞ

Least-squares regression provides the “best-ﬁt” for the model parameters if the

following assumptions are true:

(a) variability is dominant in the values of the dependent variable y (i.e. the inde-

pendent variables are precisely known),

(b) the observations are normally distributed about the y-axis,

(d) the data exhibit homoscedasticity.

(2) A point at which the gradient (or ﬁrst derivative) of a function is zero is called a

critical point. This point can either be a minimum, maximum, or saddle point of the

function.

(3) Unconstrained minimization of an objective function in one variable can be carri ed

out using Newton’s method, successive parabolic interpolation, or the golden sec-

tion search method.

(a) Newton’s method has a second-order rate of convergence and is preferred when

the ﬁrst and second derivatives of the object ive function are well-deﬁned and easy

to calculate.

(b) The golden section search method has a ﬁrst-order rate of convergence and

is preferred when the derivative of the function does not exist at one or

more points near or at the minimum, or when the derivative is difﬁcult to

calculate.

compute the function derivative and has superlinear convergence.

(4) For a single-variabl e optimization problem, the nature of the critical point can be

tested using the second derivative test. For a multivariable problem, the nature of the

critical point can be established by evaluating the matrix of second partial deriva-

tives of the objective function, called the Hessian matrix. A critical point x* is a local

minimum if Hx



ðÞis positive deﬁnite, and is a local maximum if Hx



ðÞis negative

deﬁnite.

(5) Unconstrained multidimensional optimization can be carried out using the

method of steepest descent, Newton’s multidimensional method, or the simplex

method.

(a) The steepest descent method evaluates the gradient of the objective function

to determine a search direction. An advantage of this method is that the

starting guess values do not need to be close to the minimum to attain

convergence.

(b) Newton’s multidimensional method evaluates the gradient vector and the Hessian

matrix of the function at each step. This method has a fast rate of convergence,

but convergence is guaranteed only for good initial guess values.

tions) at multiple points in the variable space. It draws a geometrical structure

called a simplex that gradually shifts towards the minimum.

533

8.6 Key points

(6) The variables of an objective function fðxÞ can have constraints placed on them.

Constraints can be of two types:

(a) equality constraints h

xðÞ¼0; j ¼ 1; 2; ...; m,or

(b) inequality constraints g

xðÞ0; k ¼ 1; 2; ...; p.

(7) The Lagrange multiplier method performs constrained minimization by converting

an optimization problem in n variables and m + p constraints into an unconstrained

problem in n + m +2p variables. The inequality constraints are converted to

equality constr aints with the use of slack variables.

(8) Bootstrapping and unders ampling are used to estimate the covariance matrix, and

thereby the conﬁdence intervals, of model parameters whose values are obtained

using nonlinear regression.

8.7 Problems

8.1. Nonlinear regression of a pharmacokinetic model The following equation

describes the plasma drug concentration for a one-compartment pharmacokinetic

model with ﬁrst-order drug absorption:

CðtÞ¼C

expðk

tÞexpðk

tÞ½:

Such a model has been successfully used to describe oral drug delivery. Given a set of

experimental data (tdat,Cdat), write a MATLAB function that can be fed into the

built-in function fminsearch to perform a nonlinear regression for the unknown

model parameters C

, k

,andk

. Try your function on the data set given in Table P8.1.

8.2. A classic test example for multidimensional minimization is the Rosenbrock banana

function:

fðxÞ¼100 x

 x



þ 1  x

ðÞ

The traditional starting point is (−1.2, 1). Start by deﬁning a function m-ﬁle that

takes in x

and x

as a single vector x and outputs the value of fðxÞ. Next, use the

built-in MATLAB function fminsearch to solve for the minimum of the above

equation, starting from the speciﬁed point of (−1.2,1). Use the optional two-output-

argument syntax of fminsearch so that the minimum x and the function value are

reported. Use your resul ts to answer the following “riddle”: when is a root-ﬁnding

problem the same as a minimization problem?

8.3. Leukocyte migration in a capillary Many times in (biomedical) engineering, after

an external stimulus or signal is applied, there is a time delay where nothing

happens ... and then a linear response is observed. One example might be a

leukocyte (white blood cell) that is squeezing through a small capillary. A positive

pressure is applied at one end of the capillary, and, foll owing a short delay, the

leukocyte begins to move slowly through the capillary at a constant velocity.

Table P8.1. Time course of drug concentration in plasma

Time (h) 0123456789101112

Plasma

conc.

(mM)

0.004 0.020 0.028 0.037 0.040 0.036 0.028 0.021 0.014 0.009 0.007 0.004 0.002

534

Nonlinear model regression and optimization

Table P8.2 provides an example of such a data set, which exhibits some random error

inherent in experimental measurements.

If you plot the data in Table 8.3 you will see that the pos ition of the cell is relat ively

constant around x =1μm until t  1 s when it begins to move in a linearly time-

dependent fashion. We want to model this data set as a line after the (unknown) time

delay. Thus, our theoretical curve is a piecewise linear function (see Figure P8.1).

To regress our data correctly to such a model we must use nonlinear regression

techniques. Here’s how it’s done.

(a) Create a function that takes in the three regression parameters (slope,

intercept, t0) as a vector called params, and outputs the sum of the

squared residuals (SSE) between the data and the model. The ﬁrst two lines of

your function should look like:

function SSE = piecewisefunc(params)

global t x

Table P8.2. Time-dependent position of a leukocyte

in a capillary subjected to a pressure gradient

Time, t(s) Position, x(μm)

0.0 0.981

0.3 1.01

0.6 1.12

0.9 1.13

1.2 2.98

1.5 2.58

1.8 3.10

2.1 5.89

Figure P8.1

Intercept

Slope

× t

+ intercept

Slope

= ?

535

8.7 Problems

Note that we are passing the data to the function as two global variables, so the same

global command should appear at the beginning of your main program. Next,

your function must calculate the SSE (sum of the squared error). Build a for loop

that executes once for each element of t. Inside your for loop, test to see whether

t(i) is greater than t

(the time delay). If it is, then the SSE is calculated in the usual

way as:

SSE = SSE + (x(i) – (slope*t(i) + intercept))^2;

If t(i) is lower than t

, then the SSE is deﬁned as the squared difference

between the data and the constant value of the model before t

. This looks like:

SSE = SSE + (x(i) – (slope*t0 + intercept))^2;

where slope*t0 + intercept is just the constant value of x before the linear

regime of behavior. Remember to initialize SSE at zero before the for loop.

(b) In your main program, you must enter the data, and declare the data as global so

that your function can access these values. Then, you may use the fminsearch

function to perform the nonlinear regression and determ ine the optimum values

of the params vector. This is accomplished with the statement:

paramsopt = fminsearch(‘piecewisefunc’, ...

[slopeguess,interceptguess,t0guess])

Note that to obtain the correct solution we must feed a good initial guess to the

program. You can come up with suitable guesses by looking at your plot of the data.

Window. Also, prepare a plot comparing the data and best-ﬁt model.

8.4. Optimizing drug dosage A patient is being treated for a certain disease simulta-

neously with two different drugs: x units of the ﬁrst drug and y units of the second

drug are administered to the patient. The body’s response R to the drugs is measured

and is modeled as a function of x and y as

fx; yðÞ¼x

y 15  x  yðÞ:

Find the values of x and y that maximize R using the method of steepest ascent. Since

you are maximizing an objective function, you will want solve the one-dimensional

problem given by

max

f x

þ αrf x

ðÞðÞ

to yield the iterative algorithm

ðkþ1Þ

¼ x

ðkÞ

þ α

opt

rf x

ðkÞ



Carry out the line search (i.e. ﬁnd α

opt

) using successive parabolic interpolation. To

maximize fðx; yÞ, at each iteration you will need to minimize

h αðÞ¼f x

ðkÞ

þ αrf x

ðkÞ



Use an initial guess of x

=(x

, y

) = (8, 2).

(a) To ensure that the minimization is successful, you will need to choose the

maximum allowable value of the step size α carefully. When the magnitude of

the gradient is large, even a small step size can move the solution very far from

the original starting point. To ﬁnd an acceptable interval for α, compute h αðÞfor

536

Nonlinear model regression and optimization