Section 10 Chemical engineering -- a chemical kinetics problem

This example may strike you as being fairly complicated. If so, you should compare it with the original solution using the so-called sensitivity equations. The reference is Bard, Nonlinear Parameter Estimation, chapter 8. We consider the chemical kinetics problem introduced on page 233. This is a model defined by a first order system of two ordinary differential equations.

The differential equations describe the evolution over time of the concentrations of the two reactants, , and . There are ten initial parameters in the model, . is the temperature at which the reaction takes place. To integrate the system of differential equations we require the initial concentrations of the reactants, and at time 0.

The reaction was carried out three times at temperatures of 200, 400, and 600 degrees. For the first run there were initially equal concentrations of the two reactants. The second run initially consisted of only the first reactant, and the third run initially consisted of only the second reactant. The initial concentrations of the reactants are known only approximately. They are

The unknown initial concentrations are treated as parameters to be estimated with Bayesian prior distributions on them reflecting the level of certainty of their true values which we have. The concentrations of the reactants were not measured directly. Rather the mixture was analyzed by a ``densitometer'' whose response to the concentrations of the reactants is

where

and

. The differences between the predicted and observed responses of the densitometer are assumed to be normally distributed so that least squares is used to fit the model. Bard employs an ``explicit'' method for integrating these differential equations, that is, the equations are approximated by a finite difference scheme like

over the time period

of length

. Equations 2 are called explicit because the values of

and

at time

are given explicitly in terms of the values of

and

at time

The advantage of using an explicit scheme for integrating the model differential equations is that the derivatives of the model functions with respect to the model parameters also satisfy differential equations -- called sensitivity equations (Bard pg 227-229). It is possible to integrate these equations as well as the model equations to get values for the derivatives. However this involves generating a lot of extra code as well as carrying out a lot of extra calculations. Since with AD Model Builder it is not necessary to produce any code for derivative calculations it is possible to employ alternate schemes for integrating the differential equations.

Let , , and In terms of and we can replace explicit finite difference scheme by the semi-implicit scheme

Now let

and solve equations (3) for

and

to obtain

Implicit and semi-implicit schemes tend to be more stable than explicit schemes over large time steps and large values of some of the model parameters. This stability is especially important when fitting nonlinear models because the algorithms for function minimization will pick very large or ``bad'' values of the parameters from time to time and the minimization procedure will generally perform better when a more stable scheme is employed.


DATA_SECTION
  init_matrix Data(1,10,1,3)  
  init_vector T(1,3) // the initial temperatures for the three runs
  init_vector stepsize(1,3)  // the stepsize to use for numerical integration
  matrix data(1,3,1,10)
  matrix sample_times(1,3,1,10) // times at which reaction was sampled
  vector x0(1,3)           // the beginning time for each of the three
                           // runs
  vector x1(1,3)           // the ending time for each of the three runs
  // for each of the three runs

PARAMETER_SECTION
  init_vector theta(1,10)   // the model parameters
  matrix init_conc(1,3,1,2) // the initial concentrations of the two
                            // reactants over three time periods
  vector instrument(1,2)    // determines the response of the densitometer
  matrix y_samples(1,10,1,2)// the predicted concentrations of the two
                            // reactants at the ten sampling periods
                            // obtained by integrating the differential
                            // equations
  vector diff(1,10)         // the difference between the observed and
                            // readings of the densitometer
  objective_function_value f  // the log_likelihood function
  number bayes_part         // the Bayesian contribution
  number y2
  number x_n
  vector y_n(1,2)
  vector y_n1(1,2)
  number A  // A B C D hold some common subexpressions
  number B
  number C
  number D
PRELIMINARY_CALCS_SECTION
  data=trans(Data);   // it is more convenient to work with the transformed
                   // matrix
PROCEDURE_SECTION

    // set up the begining and ending times for the three runs
    x0(1)=0;
    x1(1)=90;
    x0(2)=0;
    x1(2)=18;
    x0(3)=0;
    x1(3)=4.5;
    // set up the sample times for each of the three runs
    sample_times(1).fill_seqadd(0,10);  // fill with 0,10,20,...,90
    sample_times(2).fill_seqadd(0,2);   // fill with 0,2,4,...,18
    sample_times(3).fill_seqadd(0,0.5); // fill with 0,0.5,1.0,...,4.5

    // set up the initial concentrations of the two reactants for
    // each of the three runs
    init_conc(1,1)=theta(5);
    init_conc(1,2)=theta(6);
    init_conc(2,1)=theta(7);
    init_conc(2,2)=0.0;     // the initial concentrations is known to be 0
    init_conc(3,1)=0.0;     // the initial concentrations is known to be 0
    init_conc(3,2)=theta(8);

    // coefficients which determine the response of the densitometer
    instrument(1)=theta(9);
    instrument(2)=theta(10);
    f=0.0;
    for (int run=1;run<=3;run++)
    {
       // integrate the differential equations to get the predicted
       // values for the y_samples
      int nstep=(x1(run)-x0(run))/stepsize(run);
      nstep++;
      double h=(x1(run)-x0(run))/nstep; // h is the stepsize for integration

      int is=1;
      // get the initial conditions for this run
      x_n=x0(run);
      y_n=init_conc(run);
      for (int i=1;i<=nstep+1;i++)
      {
        // gather common subexpressions
        y2=y_n(2)*y_n(2);
        A=theta(1)*exp(-theta(2)/T(run));
        B=exp(-1000/T(run));
        C=(1+theta(3)*exp(-theta(4)/T(run))*y_n(1));
        C=C*C;
        D=h*A/C;
        // get the y vector for the next time step
        y_n1(1)=(y_n(1)+D*B*y2)/(1.+D);
        y_n1(2)=(y_n(2)+2.*D*y_n(1))/(1.+(2*D*B*y_n(2)));

        // if an observation occurred during this time period save
        // the predicted value
        if (is <=10)
        {
          if (x_n<=sample_times(run,is) && x_n+h >= sample_times(run,is))
          {
            y_samples(is++)=y_n;
          }
        }
        x_n+=h;  // increment the time step
        y_n=y_n1;  // update the y vector for the next step

      }
      diff=(1.0+y_samples*instrument)-data(run); //differences between the
                 // predicted and observed values of the densitometer 
      f+=diff*diff;  // sum of squared differences
    }
    // take the log of f and multiply by nobs/2 to get log-likelihood
    f=15.*log(f); // This is (number of obs)/2. It is wrong in Bard (pg 236).

    // Add the Bayesian stuff
    bayes_part=0.0;
    for  (int i=5;i<=9;i++)
    {
      bayes_part+=(theta(i)-1)*(theta(i)-1); 
    }
    bayes_part+=(theta(10)-2)*(theta(10)-2);
    f+=1./(2.*.05*.05)*bayes_part;

AD Model Builder produces a report containing values, standard deviations, and correlation matrix of the parameter estimates. As discussed below any parameter or group of parameters can easily be included in this report. For models with a large number of parameters this report can be a bit unwieldly so options are provided to exclude parameters from the report if desired.

NEXT