Per: /* Plot result */

2012-02-14T09:26:20Z

Plot result

← Older revision		Revision as of 04:26, 14 February 2012
Line 200:		Line 200:

	[[File:tubularReactor_01.png]]		[[File:tubularReactor_01.png]]

			[[Category:PROPT Examples]]

Mbot at 12:13, 9 November 2011

2011-11-09T12:13:32Z

Mbot: Created page with "{{Part Of Manual|title=the PROPT Manual|link=PROPT Manual}} Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seide..."

2011-11-08T13:09:04Z

Created page with "{{Part Of Manual|title=the PROPT Manual|link=PROPT Manual}} Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seide..."

New page

{{Part Of Manual|title=the PROPT Manual|link=[[PROPT|PROPT Manual]]}}

Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seider

Case Study V: Global optimization of a bifunctional catalyst blend in a tubular reactor

==Problem Description==

Luus et al and Luus and Bojkov consider the optimization of a tubular reactor in which methylcyclopentane is converted into benzene. The blend of two catalysts, for hydrogenation and isomerization is described by the mass fraction u of the hydrogenation catalyst. The optimal control problem is to find the catalyst blend along the length of the reactor defined using a characteristic reaction time t in the interval 0 <= t <= t_f where t_f = 2000 gh/mol corresponding to the reactor exit such that the benzene concentration at the exit of the reactor is maximized.

Find u(t) to maximize:

<math> J = x_7(t_f) </math>

Subject to:

<math> \frac{dx_1}{dt} = -k_1*x_1 </math>
<math> \frac{dx_2}{dt} = k_1*x_1-(k_2+k_3)*x_2+k_4*x_5 </math>
<math> \frac{dx_3}{dt} = k_2*x_2 </math>
<math> \frac{dx_4}{dt} = -k_6*x_4+k_5*x_5 </math>
<math> \frac{dx_5}{dt} = k_3*x_2+k_6*x_4-(k_4+k_5+k_8+k_9)*x_5+k_7*x_6+k_{10}*x_7 </math>
<math> \frac{dx_6}{dt} = k_8*x_5-k_7*x_6 </math>
<math> \frac{dx_7}{dt} = k_9*x_5-k_{10}*x_7 </math>

where xi, i = 1,..,7 are the mole fractions of the chemical species (i = 1 for methylcyclopentane, i = 7 for benzene), the rate constants are functions of the catalyst blend u(t):

<math> k_i = c(i,1)+c(i,2)*u+c(i,3)*u^2+c(i,4)*u^3 </math>

and the coefficients cij are given by Luus et al. The upper and lower bounds on the mass fraction of the hydrogenation catalyst are.

<math> 0.6 <= u <= 0.9 </math>

The initial vector of mole fractions is:

<math> x(t_0) = [1 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]' </math>

<source lang="matlab">
% Copyright (c) 2007-2008 by Tomlab Optimization Inc.

% Data
c = [0.2918487e-2, -0.8045787e-2, 0.6749947e-2, -0.1416647e-2;...
0.9509977e1, -0.3500994e2, 0.4283329e2, -0.1733333e2; ...
0.2682093e2, -0.9556079e2, 0.1130398e3, -0.4429997e2; ...
0.2087241e3, -0.7198052e3, 0.8277466e3, -0.3166655e3; ...
0.1350005e1, -0.6850027e1, 0.1216671e2, -0.6666689e1; ...
0.1921995e-1, -0.7945320e-1, 0.1105666, -0.5033333e-1;...
0.1323596, -0.4696255, 0.5539323, -0.2166664; ...
0.7339981e1, -0.2527328e2, 0.2993329e2, -0.1199999e2; ...
-0.3950534, 0.1679353e1, -0.1777829e1, 0.4974987; ...
-0.2504665e-4, 0.1005854e-1, -0.1986696e-1, 0.9833470e-2]';

% Independent variable
toms t

% Initial guess
x1opt = 1; x2opt = 0;
x3opt = 0; x4opt = 0;
x5opt = 0; x6opt = 0;
x7opt = 0; uopt = 0.9;

% Loop over number of collocation points
for n=[10 80]
</source>

==Problem setup==

<source lang="matlab">
p = tomPhase('p', t, 0, 2000, n);
setPhase(p);

tomStates x1 x2 x3 x4 x5 x6 x7
tomControls u

% Collocate initial guess
x0 = {icollocate({x1 == x1opt; x2 == x2opt
x3 == x3opt; x4 == x4opt; x5 == x5opt
x6 == x6opt; x7 == x7opt})
collocate(u == uopt)};

% Box constraints
cbox = {icollocate({
0 <= x1 <= 1; 0 <= x2 <= 1
0 <= x3 <= 1; 0 <= x4 <= 1
0 <= x5 <= 1; 0 <= x6 <= 1
0 <= x7 <= 1})
0.6 <= collocate(u) <= 0.9};

% Boundary constraints
cbnd = initial({x1 == 1; x2 == 0
x3 == 0; x4 == 0; x5 == 0
x6 == 0; x7 == 0});

% ODEs and path constraints
uvec = [1, u, u.^2, u.^3];
k = (c'*uvec')';
ceq = collocate({
dot(x1) == -k(1)*x1
dot(x2) == k(1)*x1-(k(2)+k(3))*x2+k(4)*x5
dot(x3) == k(2)*x2
dot(x4) == -k(6)*x4+k(5)*x5
dot(x5) == k(3)*x2+k(6)*x4-(k(4)+k(5)+k(8)+k(9))*x5+...
k(7)*x6+k(10).*x7
dot(x6) == k(8)*x5-k(7)*x6
dot(x7) == k(9)*x5-k(10)*x7});

% Objective
objective = -final(x7);
options = struct;
options.name = 'Tubular Reactor';
options.d2c = false;
Prob = sym2prob('con',objective,{cbox, cbnd, ceq}, x0, options);
Prob.xInit = 35; % Use 35 random starting points.
Prob.GO.localSolver = 'snopt';
if n<=10
Result = tomRun('multiMin', Prob, 1);
else
Result = tomRun('snopt', Prob, 1);
end
solution = getSolution(Result);

% Store optimum for use in initial guess
x1opt = subs(x1,solution);
x2opt = subs(x2,solution);
x3opt = subs(x3,solution);
x4opt = subs(x4,solution);
x5opt = subs(x5,solution);
x6opt = subs(x6,solution);
x7opt = subs(x7,solution);
uopt = subs(u,solution);
</source>

<pre>
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Tubular Reactor - Trial 10 f_k -0.010093483541388937
sum(|constr|) 0.000000004530098742
f(x_k) + sum(|constr|) -0.010093479011290195

Solver: multiMin with local solver snopt. EXIT=0. INFORM=0.
Find local optima using multistart local search
Did 35 local tries. Found 1 global, 31 minima. TotFuncEv 1488. TotConstrEv 1418

FuncEv 1488 GradEv 1418 ConstrEv 1418 ConJacEv 23 Iter 907
CPU time: 2.480416 sec. Elapsed time: 2.489000 sec.

</pre>

<pre>
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license 999007. Valid to 2011-12-31
=====================================================================================
Problem: --- 1: Tubular Reactor f_k -0.010082374603677158
sum(|constr|) 0.000000000037585022
f(x_k) + sum(|constr|) -0.010082374566092137
f(x_0) -0.010093483541388937

Solver: snopt. EXIT=0. INFORM=1.
SNOPT 7.2-5 NLP code
Optimality conditions satisfied

FuncEv 46 GradEv 44 ConstrEv 44 ConJacEv 44 Iter 31 MinorIter 531
CPU time: 2.059213 sec. Elapsed time: 2.055000 sec.

</pre>

<source lang="matlab">
end
</source>

==Plot result==

<source lang="matlab">
t = collocate(t);
x7 = collocate(x7opt);
u = collocate(uopt);

subplot(2,1,1)
plot(t,x7,'*-');
legend('x7');
title('Tubular Reactor state variable x7');

subplot(2,1,2)
plot(t,u,'+-');
legend('u');
title('Tubular Reactor control');
</source>

[[File:tubularReactor_01.png]]

← Older revision		Revision as of 07:13, 9 November 2011
Line 12:		Line 12:

	<math> J = x_7(t_f) </math>		<math> J = x_7(t_f) </math>


	Subject to:		Subject to:

	<math> \frac{dx_1}{dt} = -k_1*x_1 </math>		<math> \frac{dx_1}{dt} = -k_1*x_1 </math>

	<math> \frac{dx_2}{dt} = k_1x_1-(k_2+k_3)x_2+k_4*x_5 </math>		<math> \frac{dx_2}{dt} = k_1x_1-(k_2+k_3)x_2+k_4*x_5 </math>

	<math> \frac{dx_3}{dt} = k_2*x_2 </math>		<math> \frac{dx_3}{dt} = k_2*x_2 </math>

	<math> \frac{dx_4}{dt} = -k_6x_4+k_5x_5 </math>		<math> \frac{dx_4}{dt} = -k_6x_4+k_5x_5 </math>

	<math> \frac{dx_5}{dt} = k_3x_2+k_6x_4-(k_4+k_5+k_8+k_9)x_5+k_7x_6+k_{10}*x_7 </math>		<math> \frac{dx_5}{dt} = k_3x_2+k_6x_4-(k_4+k_5+k_8+k_9)x_5+k_7x_6+k_{10}*x_7 </math>

	<math> \frac{dx_6}{dt} = k_8x_5-k_7x_6 </math>		<math> \frac{dx_6}{dt} = k_8x_5-k_7x_6 </math>

	<math> \frac{dx_7}{dt} = k_9x_5-k_{10}x_7 </math>		<math> \frac{dx_7}{dt} = k_9x_5-k_{10}x_7 </math>


	where xi, i = 1,..,7 are the mole fractions of the chemical species (i = 1 for methylcyclopentane, i = 7 for benzene), the rate constants are functions of the catalyst blend u(t):		where xi, i = 1,..,7 are the mole fractions of the chemical species (i = 1 for methylcyclopentane, i = 7 for benzene), the rate constants are functions of the catalyst blend u(t):

	<math> k_i = c(i,1)+c(i,2)u+c(i,3)u^2+c(i,4)*u^3 </math>		<math> k_i = c(i,1)+c(i,2)u+c(i,3)u^2+c(i,4)*u^3 </math>


	and the coefficients cij are given by Luus et al. The upper and lower bounds on the mass fraction of the hydrogenation catalyst are.		and the coefficients cij are given by Luus et al. The upper and lower bounds on the mass fraction of the hydrogenation catalyst are.

	<math> 0.6 <= u <= 0.9 </math>		<math> 0.6 <= u <= 0.9 </math>


	The initial vector of mole fractions is:		The initial vector of mole fractions is:

	<math> x(t_0) = [1 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]' </math>		<math> x(t_0) = [1 \ 0 \ 0 \ 0 \ 0 \ 0 \ 0]' </math>


	<source lang="matlab">		<source lang="matlab">

PROPT Tubular Reactor - Revision history

Per: /* Plot result */

Mbot at 12:13, 9 November 2011

Mbot: Created page with "{{Part Of Manual|title=the PROPT Manual|link=PROPT Manual}} Global Optimization of Chemical Processes using Stochastic Algorithms 1996, Julio R. Banga, Warren D. Seide..."