pfrtempopt_c.m

Matlab在化学工程中的应用

function PFRTempOpt_C
% 求反应管最优温度分布(用最优化方法搜索k，并用函数bvp4c求解ODE-BVP问题)
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/06/04 $

% 已知参数
k0 = 1; 
Z = 1;
tspan = [0.00  0.12  0.24  0.36  0.50  0.62  0.74  0.86  1.00];
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;
A = zeros(9,9);
b = zeros(9,1);
for i = 1:8
    A(i,i) = -1;
    A(i,i+1) = 1;
end

k0 = k;
k = fmincon(@ObjFunc,k0,A,b,[],[],lb,ub,[],[],tspan)    % 优化

% 结果输出
solinit = bvpinit(tspan,[1 0 0.2 0.3]);
sol = bvp4c(@ODEs,@BCfun,solinit,[],tspan,k);
x = deval(sol,tspan);
data = [tspan' x' k'];
fprintf('\n\tResults:\n')
fprintf('\t\tt\t\tx1\t\tx2\tλ_1\t\tλ_2hat\tk\n')
for i = 1:length(tspan)
    fprintf('\t%.3f',data(i,:))
    fprintf('\n')
end


% --------------------------------------------------------------------------
function f = ObjFunc(k,tspan)                       % 目标函数
solinit = bvpinit(tspan,[1 0 0.2 0.3]);
sol = bvp4c(@ODEs,@BCfun,solinit,[],tspan,k);
x = deval(sol,tspan);   % x(3,:)=rambda(1,:), x(4,:)=rambda(2,:)
f = mean(abs(x(1,:).*(x(4,:)-x(3,:)) + 2*x(2,:).*k.*(x(3,:)-3*x(4,:))))

% --------------------------------------------------------------------------
function dxdt = ODEs(t,x,tspan,k)                   % ODE方程
k = spline(tspan,k,t);
rambda(1) = x(3);
rambda(2) = x(4);
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(3) = k*(rambda(1)-rambda(2));
dxdt(4) = k^2*(3*rambda(2)-rambda(1));

dxdt = dxdt';

% --------------------------------------------------------------------------
function bc = BCfun(ya,yb,tspan,k)                  % 边界条件
bc = [ya(1)-1; ya(2); yb(3); yb(4)-1];