function DOE
% D-優化試驗設計
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/07/12 $

clear all
clc

% 各因素的水平
p = [470, 300,120
     285, 190, 65
     100, 80, 10];

% D-優化試驗設計：先生成settings，再將settings轉變為相應實驗條件expCond
settings = cordexch(3,13,'q');      % settings: 因子設置矩陣
mr = p(2,:);                        % mr: middle row, i.e., middle level
mr = mr(ones(13,1),:);
hr = (p(1,:) - p(3,:))/2;
hr = hr(ones(13,1),:);
expCond = settings.*hr + mr;        % expCond: settings的相應實驗條件

% 由反應模擬器生成實驗數據data
p1 = expCond(:,1);
p2 = expCond(:,2);
p3 = expCond(:,3);
y = zeros(13,1);
for k = 1:13
    y(k)  = 1.25*(p2(k) - p3(k)/1.5183)./(1 + 0.064*p1(k)     ...
            + 0.0378*p2(k) + 0.1326*p3(k))*normrnd(1,0.02);
end
data = [expCond y]     % data為實驗數據矩陣（實驗條件及其對應的反應速率）

% 數據分析: 由非線性模型Nonlinear Model估計參數
x = data(:,1:3); 
y = data(:,4);
xname = str2mat('Hydrogen','n-Pentane','Isopentane');
yname = 'Reaction Rate';
beta0 = [1.2 0.1 0.01 0.1 1.5];             % 參數初值
nlintool(x,y,@hougen,beta0,[],xname,yname);
% rstool(x,y,[],[],xname,yname);
rstool(x,y,'quadratic',[],xname,yname);

[beta,resid,j] = nlinfit(x,y,@hougen,beta0)
ci = nlparci(beta,resid,j)

% 參數辨識結果:β1、β2、... β5
fprintf('Estimated Parameters:\n')
fprintf('\tβ1 = %.2f ± %.2f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tβ2 = %.3f ± %.3f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tβ3 = %.4f ± %.4f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tβ4 = %.4f ± %.4f\n',beta(4),ci(4,2)-beta(4))
fprintf('\tβ5 = %.4f ± %.4f\n',beta(5),ci(5,2)-beta(5))