function PFRTempOpt_B
% 求反應管最優溫度分布(用最優化方法搜索k，其它同PFRTempOpt_A.m)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/06/03 $

% 已知參數
k0 = 1; 
Z = 1;
tspan1 = [0.00  0.12  0.24  0.36  0.50  0.62  0.74  0.86  1.00];
tspan2 = inverse(tspan1);
k = [1  0.5  0.5  0.5  0.5  0.5  0.5  0.5  0.5];
lb = zeros(size(k));
ub = ones(size(k));
lb(1) = 1; 
ub(1) = 1;
k0 = k;
A = zeros(9,9);
b = zeros(9,1);
for i = 1:8
    A(i,i) = -1;
    A(i,i+1) = 1;
end
k = fmincon(@ObjFunc,k0,A,b,[],[],lb,ub,[],[],tspan1,tspan2)
[t,x] = ode45(@CEquation,tspan1,[1.0 0.0],[],k,tspan1); 

% 結果輸出
data = [t x k'];
fprintf('\n\tResults:\n')
fprintf('\t\tt\t\tx1\t\tx2\t\tk\n')
for i = 1:length(t)
    fprintf('\t%.3f',data(i,:))
    fprintf('\n')
end


% ------------------------------------------------------------------
function f = ObjFunc(k,tspan1,tspan2)           % 目標函數
[t,x] = ode45(@CEquation,tspan1,[1.0 0.0],[],k,tspan1);    
k2 = inverse(k);
[t,rambda] = ode45(@rambdaEquation,tspan2,[0 1],[],k2,tspan2);
rambda = inverse(rambda);
f = mean(abs(x(:,1).*(rambda(:,2)-rambda(:,1))+2*x(:,2).*k'.*(rambda(:,1)-3*rambda(:,2))))

% ------------------------------------------------------------------
function dxdt = CEquation(t,x,ki,tspan)                 % 濃度方程
k = spline(tspan,ki,t);
k1 = k;
k2 = k*k;
k3 = 2*k*k;
dxdt(1) = k2*x(2)-k1*x(1);
dxdt(2) = k1*x(1)-(k2+k3)*x(2);
dxdt = dxdt';

% ------------------------------------------------------------------
function f = rambdaEquation(t,rambda,k2,tspan)          % 伴隨方程
k = spline(tspan,k2,t);
f(1) = k*(rambda(1)-rambda(2));
f(2) = k^2*(3*rambda(2)-rambda(1));
f = f';

% ------------------------------------------------------------------
function y = inverse(x)                                 % 逆向排序
[row,col] = size(x);
if row==1
    i = col:-1:1;
    j = 1:col;
    y(j) = x(i);
elseif col==1
    i = row:-1:1;
    j = 1:row;
    y(j) = x(i);
else
    i = row:-1:1;
    j = 1:row;
    y(j,:) = x(i,:);
end