Chapter Contents 
Previous 
Next 
Nonlinear Optimization Examples 
In this class of problems, function evaluations, as well as gradient evaluations, are not done in full precision. Evaluating a function involves the numerical solution of a differential equation with some prescribed precision. Therefore, two choices exist for evaluating first and secondorder derivatives:
With the finite difference approach, a very delicate balance of all the precision requirements of every routine must exist. In the examples that follow, notice the relative levels of precision that are imposed on different modules. Since finite differences are used to compute the first and secondorder derivatives, it is incorrect to set the precision of the ODE solver at a coarse level because that would render the numerical estimation finite difference worthless.
A coarse computation of the solution of the differential equation cannot be accompanied by very fine computation of the finite difference estimates of the gradient and the Hessian. That is, you cannot set the precision of the differential equation solver to be 1E4 and perform the finite difference estimation with a precision of 1E10. In addition, this precision must be wellbalanced with the termination criteria imposed on the optimization process.
In general, if the precision of the function evaluation is , the gradient should be computed by finite differences , and the Hessian should be computed with finite differences .^{*}
The rates of flow of the drug are described by the following pair of ordinary differential equations:
Also, a nonnegativity constraint is imposed on the parameters , , and , although for numerical purposes, you may need to use a small value instead of zero for these bounds (such as 1E7).
Suppose y_{i} is the observed serum concentration at time t_{i}. The parameters are estimated by minimizing the sum of squares of the differences between the observed and predicted serum concentrations:
data tetra; input t c @@; datalines; 1 0.7 2 1.2 3 1.4 4 1.4 6 1.1 8 0.8 10 0.6 12 0.5 16 0.3 ; proc iml; use tetra; read all into tetra; start f(theta) global(thmtrx,t,h,tetra,eps); thmtrx = ( theta[1]  0 ) // ( theta[1]  theta[2] ); c = theta[3]//0 ; t = 0 // tetra[,1]; call ode( r1, "der",c , t, h) j="jac" eps=eps; f = ssq((r1[2,])`tetra[,2]); return(f); finish; start der(t,x) global(thmtrx); y = thmtrx*x; return(y); finish; start jac(t,x) global(thmtrx); y = thmtrx; return(y); finish; h = {1.e14 1. 1.e5}; opt = {0 2 0 1 }; tc = repeat(.,1,12); tc[1] = 100; tc[7] = 1.e8; par = { 1.e13 . 1.e10 . . . . }; con = j(1,3,0.); itheta = { .1 .3 10}; eps = 1.e11; call nlpdd(rc,rx,"f",itheta) blc=con opt=opt tc=tc par=par;
The output from the optimization process is shown in Output 11.3.1.
Output 11.3.1: Printed Output for Tetracycline Diffusion Problem
Output 11.3.1: (continued)

The closed form of the solution requires a change to the function f(·). The functions needed as arguments of the ODE routine, namely the der and jac modules, can be removed.
start f(th) global(theta,tetra); theta = th; vv = v(tetra[,1]); error = ssq(vvtetra[,2]); return(error); finish; start v(t) global(theta); v = theta[3]*theta[1]/(theta[2]theta[1])* (exp(theta[1]*t)exp(theta[2]*t)); return(v); finish; call nlpdd(rc,rx,"f",itheta) blc=con opt=opt tc=tc par=par;
The parameter estimates, which are shown in Output 11.3.2, are close to those obtained by the first solution.
Output 11.3.2: Second Set of Parameter Estimates for Tetracycline Diffusion

data oil ( drop=temp); input temp time bitumen oil; datalines; 673 5 0. 0. 673 7 2.2 0. 673 10 11.5 0.7 673 15 13.7 7.2 673 20 15.1 11.5 673 25 17.3 15.8 673 30 17.3 20.9 673 40 20.1 26.6 673 50 20.1 32.4 673 60 22.3 38.1 673 80 20.9 43.2 673 100 11.5 49.6 673 120 6.5 51.8 673 150 3.6 54.7 ; proc iml; use oil; read all into a; /****************************************************************/ /* The INS function inserts a value given by FROM into a vector */ /* given by INTO, sorts the result, and posts the global matrix */ /* that can be used to delete the effects of the point FROM. */ /****************************************************************/ start ins(from,into) global(permm); in = into // from; x = i(nrow(in)); permm = inv(x[rank(in),]); return(permm*in); finish; start der(t,x) global(thmtrx,thet); y = thmtrx*x; if ( t <= thet[5] ) then y = 0*y; return(y); finish; start jac(t,x) global(thmtrx,thet); y = thmtrx; if ( t <= thet[5] ) then y = 0*y; return(y); finish; start f(theta) global(thmtrx,thet,time,h,a,eps,permm); thet = theta; thmtrx = ((theta[1]+theta[4])  0  0 )// (theta[1]  (theta[2]+theta[3])  0 )// (theta[4]  theta[2]  0 ); t = ins( theta[5],time); c = { 100, 0, 0}; call ode( r1, "der",c , t , h) j="jac" eps=eps; /* send the intermediate value to the last column */ r = (c r1) * permm; m = r[2:3,(2:nrow(time))]; mm = m` a[,2:3]; call qr(q,r,piv,lindep,mm); v = det(r); return(abs(v)); finish; opt = {0 2 0 1 }; tc = repeat(.,1,12); tc[1] = 100; tc[7] = 1.e7; par = { 1.e13 . 1.e10 . . . .}; con = j(1,5,0.); h = {1.e14 1. 1.e5}; time = (0 // a[,1]); eps = 1.e5; itheta = { 1.e3 1.e3 1.e3 1.e3 1.}; call nlpqn(rc,rx,"f",itheta) blc=con opt=opt tc=tc par=par;
The parameter estimates are shown in Output 11.3.3.
Output 11.3.3: Parameter Estimates for Oil Shale Pyrolysis

start f(thet) global(time,a); do i = 1 to nrow(time); t = time[i]; v1 = 100; if ( t >= thet[5] ) then v1 = 100*ev(t,thet[1],thet[4],thet[5]); v2 = 0; if ( t >= thet[5] ) then v2 = 100*thet[1]/(thet[2]+thet[3]thet[1]thet[4])* (ev(t,thet[1],thet[4],thet[5]) ev(t,thet[2],thet[3],thet[5])); v3 = 0; if ( t >= thet[5] ) then v3 = 100*thet[4]/(thet[1]+thet[4])* (1.  ev(t,thet[1],thet[4],thet[5])) + 100*thet[1]*thet[2]/(thet[2]+thet[3]thet[1]thet[4])*( (1.ev(t,thet[1],thet[4],thet[5]))/(thet[1]+thet[4])  (1.ev(t,thet[2],thet[3],thet[5]))/(thet[2]+thet[3]) ); y = y // (v1  v2  v3); end; mm = y[,2:3]a[,2:3]; call qr(q,r,piv,lindep,mm); v = det(r); return(abs(v)); finish; start ev(t,a,b,c); return(exp((a+b)*(tc))); finish; con = { 0. 0. 0. 0. . . . , . . . . . . . , 1 1 1 1 . 1 1.e7 }; time = a[,1]; par = { 1.e13 . 1.e10 . . . .}; itheta = { 1.e3 1.e3 1.e2 1.e3 1.}; call nlpqn(rc,rx,"f",itheta) blc=con opt=opt tc=tc par=par;
The parameter estimates are shown in Output 11.3.4.
Output 11.3.4: Second Set of Parameter Estimates for Oil Shale Pyrolysis

Chapter Contents 
Previous 
Next 
Top 
Copyright © 1999 by SAS Institute Inc., Cary, NC, USA. All rights reserved.