Contents
- Demonstration of Profiled Estimation of the Nylon Differential Equations
- Set up function paths
- Load the data
- Reformat the data into cell arrays
- Plot observations and order the plots by temperature
- Profiling optimisation control
- Setting up functional data objects and perform a model based smooth.
- Perform the Profiled Estimation
- plot the results
- Calculate a Sample Information Matrix
Demonstration of Profiled Estimation of the Nylon Differential Equations
This page demonstrates the generalized profiling estimation for the nylon differential equations. The nylon data is from Zheng et al 2005 using a 4 parameter version of their differential equations and rescaled parameters as described in Campbell et al. 2007. These equations include calculations of $W_{eq}$ and $C_T$ from the input water pressure using equations from Schaffer et al. 2003.
Wei Zheng, Kim B. McAuley, E. Keith Marchildon and Yao K. Zhen, 2005, "Effects of End-Group Balance on Melt-Phase Nylon 612 Polycondensation: Experimental Study and Mathematical Model", Ind. Eng. Chem. Res. vol 44, p 2675 - 2686
M.A. Schaffer, Kim B. McAuley, M.F. Cunningham, and E. Keith Marchildon, 2003 "Experimental Study and Modeling of Nylon Polycondensation in the Melt Phase", Ind. Eng. Chem. Res. 42, p2946
%The nylon differential equation system: % %$$ \frac{dA}{dt} = -k_p_0*10^{-3}*(AC-LW/K_a)$$ %$$ \frac{dC}{dt} = -k_p_0*10^{-3}*(AC-LW/K_a)$$ %$$ \frac{dW}{dt} = k_p_0*10^{-3}*(AC-LW/K_a) - 24.3(W-W_{eq})$$ % % Where $K_p_0$ is unknown and $K_a$ is a function of the three % other parameters $[\gamma, K_{a0}, \Delta H]$, % experimentally controlled temperature $T$, reference % temperature $T_0=549.15$, ideal gas constant $R=8.3145*10^{-3}$ and % $K_{T}=20.97\exp[-9.624+3613/T]$ % % $$ K_a = (1+W_{eq}\gamma 10^{-3})K_T K_{a0}*10^{-3}* \exp\left[-\frac{\Delta H}{R}\left(1/T-1/T_0\right)\right] $$ %
Set up function paths
clear clc addpath(strcat(matlabroot,'\Profiling\nylon\')) addpath(strcat(matlabroot,'\Profiling')) odefn = @nylonfn1; % Function for ODE solver (exact), used to simulate data fn.fn = @nylonfn; % RHS function fn.dfdx = @nylondfdy; % Derivative wrt inputs (Jacobian) fn.dfdp = @nylondfdpar; % Derviative wrt parameters fn.d2fdx2 = @nylond2fdy2; % Hessian wrt inputs fn.d2fdxdp = @nylond2fdydpar; % Hessian wrt inputs and parameters fn.d2fdp2 = @nylond2fdpar2; % Hessian wrt parameters. fn.d3fdx3 = @nylond3fdy3; % Third derivative wrt inputs. fn.d3fdx2dp = @nylond3fdy2dp; % Third derivative wrt intputs, inputs and pars. fn.d3fdxdp2 = @nylond3fdydp2; % Third derivative wrt inputs, pars and pars. alg = [1,1,1]; % None of the 3 DEs are algebraic equations
Load the data
the nylon_data.txt file has columns: 1 - A number indentifying the experimental run 2 - The temperature, constant within an experimental run 3 - The observation time 4 - Input $P_w$ which must be converted to $W_eq$ using temperature 5 - Observed values of A 6 - Observed values of C
The changepoint_times.txt are in a seperate file. There are two times where the input $P_w$ undergoes a step change in each experimental run.
nylon_data = load('nylon_data.txt'); run = nylon_data(:,1); Nruns = length(unique(run)); T = nylon_data(:,2); Time = nylon_data(:,3); Pwvec = nylon_data(:,4); A = nylon_data(:,5); C = nylon_data(:,6); runs=unique(run); changes = load('changepoint_times.txt'); changepoints=cell(Nruns,1); for i=1:Nruns changepoints{i}=changes(i,2:3); end
Reformat the data into cell arrays
The observations and times must be in cell arrays to be used by the profile functions. Since $P_w$ and temperature $T$ are used to compute the step function $W_{eq}$, a functional data object is created to pass $W_{eq}$ to the profile nylon functions. This is passed along with temperature to the nylon functions through fn_extras.
Also extract the time range of the observed values, set the initial parameter estimates and set up a fine grid of points, used for plotting purposes later.
Pw=cell(Nruns,1);
Tcell = cell(Nruns,3);
Ycell = cell(Nruns,3);
tfine = cell(Nruns,3);
Temps = cell(Nruns,1);
range = cell(Nruns,1);
C_time_index = cell(Nruns,1);
index={};
for i = 1:Nruns
index{i}=find(run==i);
Tcell{i,1} = Time(index{i}) ; %observation times for A
Tcell{i,2} = Tcell{i,1}(~isnan(C(index{i}))); %Cobservation times for C
Tcell{i,3}=[]; % component W is unobserved
Ycell{i,1} = A(index{i}(~isnan(A(index{i})))); %observations for A
C_time_index{i} = ~isnan(C(index{i}));
Ycell{i,2} = C(index{i}(C_time_index{i})); %observations for C
Pw{i}=Pwvec(index{i}); %Experimental conditions
Temps{i}=T(index{i}(1));
range{i}=[0,max(Tcell{i,1})];
tfine(i,:) = repmat({linspace(0,Tcell{i,1}(end),2000)},1,3);
end
time_temp=cell(Nruns,1);
fn_extras=cell(Nruns,1);
W_eq = nylon_W_eq_fd(Temps, Pw,range, cell2mat(changepoints));
for i = 1:Nruns
time_temp{i}=[min([Tcell{i,1}(1),Tcell{i,2}(1)]),changepoints{i},max([Tcell{i,1}(end),Tcell{i,2}(end)])];
end
W_eq_temp = eval_fdcell(time_temp,W_eq,0);
fn_extras=cell(6,1);
for i = 1:Nruns
fn_extras{i}={Temps{i},W_eq{i}};
end
% Initial Parameter Estimates
startpars=[17.7,17.0,15,-78.1];
parind=repmat(1:length(startpars),Nruns,1);
Plot observations and order the plots by temperature
Observations are grouped within experiments, each experiment is seperated by lines on the plot. Plots are sorted by temperature, $T$, which is listed above the plot. The subscript on $T$ in the plot is the run # identifyer.
Temps=cell2mat(Temps); [ii,sort_index]=sort(Temps); dex=[12,10,11,3,2,1]; for k = 1:Nruns i=sort_index(k); subplot(Nruns,3,dex(k)) plot(Tcell{i,1},Ycell{i,1},'.') xlim([0,Tcell{i,1}(end)]) if(dex(k)==1 ||dex(k)==10) ylabel('\fontsize{18}A') end title(strcat('T_',num2str(runs(i)),'=',num2str(Temps(i))),'FontSize',18) subplot(Nruns,3,dex(k)+3) plot(Tcell{i,2},Ycell{i,2},'.') xlim([0,Tcell{i,1}(end)]) if(dex(k)==1 ||dex(k)==10) ylabel('\fontsize{18}C') end subplot(Nruns,3,dex(k)+6) plot([0,changepoints{i}(1),changepoints{i},changepoints{i}(2),Tcell{i,1}(end)],... [W_eq_temp{i}(1);W_eq_temp{i}(1:2);W_eq_temp{i}(2:3);W_eq_temp{i}(3)],'-') xlim([0,Tcell{i,1}(end)]) if(dex(k)==1 ||dex(k)==10) ylabel('\fontsize{18}W_e_q') end end annotation(gcf,'line',[0.36 0.36],[0.95 0.05]); annotation(gcf,'line',[0.05 0.95],[0.51 0.51]); annotation(gcf,'line',[0.65 0.65],[0.95 0.05]);
Profiling optimisation control
Set the convergence criteria for the inner and outer optimization loops.
lsopts_out = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','iter','MaxIter',500,'TolFun',1e-8,'TolX',1e-8); lsopts_in = optimset('DerivativeCheck','off','Jacobian','on',... 'Display','off','MaxIter',500,'TolFun',1e-8,'TolX',1e-8,... 'JacobMult',@SparseJMfun);
Setting up functional data objects and perform a model based smooth.
Set up the B-spline knots, quadrature points and bases.
We use a fifth order b-spline basis with knots at each observation point and additional knots at a rate of 5 per hour of experimental time, keeping only unique knots from this process. Furthermore additional knots are used to allow a break in the first derivative at the step changes in $P_w$.
We also get an initial estimate of the spline coefficients using a penalty on the second derivatives.
In Zheng et al. (2005), the approximate measurement error was estimated from repeated observations in an additional experiment. The inverse of this expected measurement variance is included in the term $wts$.
The break in the derivative at times of step function changes in input $P_w$ can cause confusion for the software in computing the integrated PEN term in the smoothing step. The software may a mixture of left and right hand derivatives to estimate the value at these discontinuities in the derivative. Consequently, we avoid this problem by omitting a small $\delta$ neighbouhood from the integral in the spline fitting criteria. The boundaries of this neighbourhood are shown in the plot of the smooth.
norder = 5; % Order of B-spline approximation % Set up the knot locations: knots_cell = cell(size(Tcell)); for i=1:Nruns knots_cell(i,:) = repmat({unique([Tcell{i,1}',linspace(0,Tcell{i,1}(end),5*floor(Tcell{i,1}(end)))])},1,3); for j=1:(size(Tcell,2)) break_knots=sum(knots_cell{i,j}==changepoints{i}(1)); if(norder-1-break_knots>0) knots_cell{i,j}=sort([knots_cell{i,j},repmat(changepoints{i},1,(norder-1-break_knots))]); end if(knots_cell{i,j}(end)~=range{i}(2)) knots_cell{i,j}=[knots_cell{i,j};range{i}(2)]; end end end % Set up bases LFD_pen_order=2; lambda = 1e2; lambda_first = 1e-1; rough_ord = 2; delta = 1e-6; % half the open neighbourhood around the changepoints that we omit from the integrated penalty. wts=repmat([1/0.6^2,1/2.4^2,0],Nruns,1); % this is 1/sigma_A, 1/sigma_C and 0 for [W] since it's unobserved % Make the basis and obtain an initial estimate for the b-spline % coefficients DEfd_cell=cell(size(Ycell)); for i=1:Nruns DEfd_cell(i,:) = SplineEst_with_jumps(fn,Tcell(i,:),Ycell(i,:),startpars,... knots_cell(i,:),wts,lambda,lambda_first ,rough_ord,alg,lsopts_in,... [],fn_extras{i},norder,changepoints{i},delta); end
Perform the Profiled Estimation
Start with small $\lambda$ and increase up to $\lambda=10^3$
lambda=10.^[1,2,3,4];
[newpars(:,1),newDEfd_cell{1}] = Profile_GausNewt_rep([10,10,10,10]',lsopts_out,parind,...
DEfd_cell,fn,lambda(1),Ycell,Tcell,wts,[],lsopts_in,fn_extras);
for i=2:length(lambda)
[newpars(:,i),newDEfd_cell{i}] = Profile_GausNewt_rep(newpars(:,i-1),lsopts_out,parind,...
newDEfd_cell{i-1},fn,lambda(i),Ycell,Tcell,wts,[],lsopts_in,fn_extras);
end
disp('Starting values and the parameter estimates as \lambda increases')
disp(num2str([startpars',newpars]))
Iteration steps Residual Improvement Grad-norm parameters
1 1 127177 0.809031 2.45e+003 13.53 14.801 20.9233 -12.3014
2 1 16250.3 0.872223 411 17.2116 20.1524 35.7092 -29.6118
3 1 2055.45 0.873513 58.6 19.4854 23.5516 47.8743 -37.2451
4 1 1419.36 0.309466 3.62 19.7494 24.1559 52.3389 -38.6598
5 1 1416.11 0.0022876 0.0356 19.6293 24.038 52.8638 -38.7617
6 1 1416.11 1.17032e-006 0.000108 19.6149 24.025 52.8817 -38.7511
7 18 1416.11 -9.88625e-008 0.000117 19.6149 24.025 52.8817 -38.7511
Iteration steps Residual Improvement Grad-norm parameters
1 1 1663.42 0.00398545 1.13 20.481 26.4105 50.5208 -36.549
2 1 1663.13 0.000176332 0.00442 20.5243 26.496 50.5466 -36.7287
3 1 1663.13 1.67892e-008 3.49e-005 20.5265 26.5041 50.5394 -36.7296
4 1 1663.13 8.68112e-009 1.22e-006 20.5266 26.5045 50.539 -36.7299
Iteration steps Residual Improvement Grad-norm parameters
1 1 1702.9 3.81976e-005 0.00884 20.5032 26.7082 50.3551 -36.5566
2 1 1702.9 2.54081e-008 8.13e-005 20.5027 26.7129 50.3516 -36.559
3 1 1702.9 1.10333e-006 2.48e-005 20.5028 26.7131 50.3514 -36.559
4 14 1702.9 -3.47375e-009 2.49e-005 20.5028 26.7131 50.3514 -36.559
Iteration steps Residual Improvement Grad-norm parameters
1 1 1871.1 0.0519612 8.71 17.2297 20.7774 55.6318 -40.9925
2 1 1856.26 0.00792966 0.286 18.0242 21.8718 55.1675 -39.9295
3 1 1856.06 0.000109537 0.175 17.864 21.4308 55.6563 -40.4519
4 1 1856.04 8.79626e-006 0.0287 17.9088 21.5232 55.563 -40.3343
5 1 1856.04 6.61218e-007 0.0104 17.8985 21.4977 55.5905 -40.3665
6 21 1856.24 -0.000108414 0.00542 17.8985 21.4977 55.5905 -40.3665
Starting values and the parameter estimates as \lambda increases
17.7 19.6149 20.5266 20.5028 17.8985
17 24.025 26.5045 26.7131 21.4977
15 52.8817 50.539 50.3514 55.5905
-78.1 -38.7511 -36.7299 -36.559 -40.3665
plot the results
i=find(lambda==1e3);
devals = eval_fdcell(tfine,newDEfd_cell{i},0);
[ii,sort_index]=sort(Temps);
dex=[12,10,11,3,2,1];
for k = 1:Nruns
i=sort_index(k);
subplot(Nruns,3,dex(k))
plot(tfine{i,1},devals{i,1},'-b',Tcell{i,1},Ycell{i,1},'.')
xlim([0,Tcell{i,1}(end)])
if(dex(k)==1 ||dex(k)==10)
ylabel('\fontsize{18}A')
end
title(strcat('T_',num2str(runs(i)),'=',num2str(Temps(i))),'FontSize',18)
subplot(Nruns,3,dex(k)+3)
plot(tfine{i,2},devals{i,2},'-b',Tcell{i,2},Ycell{i,2},'.')
xlim([0,Tcell{i,1}(end)])
if(dex(k)==1 ||dex(k)==10)
ylabel('\fontsize{18}C')
end
subplot(Nruns,3,dex(k)+6)
plot(tfine{i,3},devals{i,3},'-b')
xlim([0,Tcell{i,1}(end)])
if(dex(k)==1 ||dex(k)==10)
ylabel('\fontsize{18}W')
end
end
annotation(gcf,'line',[0.36 0.36],[0.95 0.05]);
annotation(gcf,'line',[0.05 0.95],[0.51 0.51]);
annotation(gcf,'line',[0.65 0.65],[0.95 0.05]);
Calculate a Sample Information Matrix
Use
i=find(lambda==1e3); % Hessian of squared error with respect to parameters d2Jdp2 = make_d2jdp2_rep(newDEfd_cell{i},fn,Ycell,Tcell,lambda(i),newpars(:,i),parind,alg,wts,fn_extras); % Second derivatives with respect to parameters and observations d2JdpdY = make_d2jdpdy_rep(newDEfd_cell{i},fn,Ycell,Tcell,lambda(i),newpars(:,i),parind,[],wts,fn_extras); % Resulting derivative of parameters with respect to observations dpdY = -d2Jdp2\d2JdpdY; % Variance of observations: S = make_sigma(newDEfd_cell{i},Tcell,Ycell,0); % Resulting parameter covariance matrix: Cov = dpdY*S*dpdY'; upper=newpars(:,i)+2*sqrt(diag(Cov)); lower=newpars(:,i)-2*sqrt(diag(Cov)); CI=[lower,newpars(:,i),upper]; disp('Approximate 95% Confidence Interval for Nylon Reaction Parameters') disp('lower, point, and upper estimates') disp(num2str(CI));
Approximate 95% Confidence Interval for Nylon Reaction Parameters lower, point, and upper estimates 17.3419 20.5028 23.6636 20.0205 26.7131 33.4058 44.0804 50.3514 56.6224 -44.1034 -36.559 -29.0146