PROPT Oil Shale Pyrolysis

This page is part of the PROPT Manual. See PROPT Manual.

Dynamic Optimization of Batch Reactors Using Adaptive Stochastic Algorithms 1997, Eugenio F. Carrasco, Julio R. Banga

Case Study II: Oil Shale Pyrolysis

Problem description

A very challenging optimal control problem is the oil shale pyrolysis case study, as considered by Luus (1994). The system consists of a series of five chemical reactions:

A1 -> A2

A2 -> A3

A1+A2 -> A2+A2

A1+A2 -> A3+A2

A1+A2 -> A4+A2

This system is described by the differential equations

${\frac {dx_{1}}{dt}}=-k_{1}*x_{1}-(k_{3}+k_{4}+k_{5})*x_{1}*x_{2}$ ${\frac {dx_{2}}{dt}}=k_{1}*x_{1}-k_{2}*x_{2}+k_{3}*x_{1}*x_{2}$ ${\frac {dx_{3}}{dt}}=k_{2}*x_{2}+k_{4}*x_{1}*x_{2}$ ${\frac {dx_{4}}{dt}}=k_{5}*x_{1}*x_{2}$

where the state variables are the concentrations of species, Ai, i = 1, ..., 4. The initial condition is

$x(t_{0})=[1\ 0\ 0\ 0]'$

The rate expressions are given by:

$k_{i}=k_{i0}*exp(-{\frac {Ei}{R*T}}),i=1,2,3,4,5$

where the values of ki0 and Ei are given by Luus (1994). The optimal control problem is to find the residence time t_f and the temperature profile T(t) in the time interval 0 <= t <= t_f so that the production of pyrolytic bitumen, given by x2, is maximized. Therefore, the performance index is

$J=x_{2}(t_{f})$

The constraints on the control variable (temperature) are:

$698.15<=T<=748.15$

% Copyright (c) 2007-2008 by Tomlab Optimization Inc.

Problem setup

toms t
toms t_f

ai = [8.86; 24.25; 23.67; 18.75; 20.70];
bi = [20300; 37400; 33800; 28200; 31000]/1.9872;

for n=[4 10 20 30 35]

    p = tomPhase('p', t, 0, t_f, n);
    setPhase(p);

    tomStates x1 x2 x3 x4
    tomControls T

    % Initial guess
    if n == 4
        x0 = {t_f == 9.3
            collocate(T == 725)};
    else
        x0 = {t_f == tfopt
            icollocate({
            x1 == x1opt; x2 == x2opt
            x3 == x3opt; x4 == x4opt
            })
            collocate(T == Topt)};
    end

    % Box constraints
    cbox = {9.1 <= t_f <= 12
        icollocate({0 <= x1 <= 1; 0 <= x2 <= 1
        0 <= x3 <= 1; 0 <= x4 <= 1})
        698.15 <= collocate(T) <= 748.15};

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

    % ODEs and path constraints

    ki1 = exp(ai(1)-bi(1)./T);
    ki2 = exp(ai(2)-bi(2)./T);
    ki3 = exp(ai(3)-bi(3)./T);
    ki4 = exp(ai(4)-bi(4)./T);
    ki5 = exp(ai(5)-bi(5)./T);

    ceq = collocate({
        dot(x1) == -ki1.*x1-(ki3+ki4+ki5).*x1.*x2
        dot(x2) == ki1.*x1-ki2.*x2+ki3.*x1.*x2
        dot(x3) == ki2.*x2+ki4.*x1.*x2
        dot(x4) == ki5.*x1.*x2});

    % Objective
    objective = -final(x2);

Solve the problem

    options = struct;
    options.name = 'Oil Pyrolysis';
    solution = ezsolve(objective, {cbox, cbnd, ceq}, x0, options);

    x1opt = subs(x1, solution);
    x2opt = subs(x2, solution);
    x3opt = subs(x3, solution);
    x4opt = subs(x4, solution);
    Topt = subs(T, solution);
    tfopt = subs(final(t), solution);

Problem type appears to be: lpcon
Time for symbolic processing: 0.5047 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Oil Pyrolysis                  f_k      -0.357327805323273240
                                       sum(|constr|)      0.000000000957541036
                              f(x_k) + sum(|constr|)     -0.357327804365732190
                                              f(x_0)      0.000000000000000000

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

FuncEv    1 ConstrEv   93 ConJacEv   93 Iter   50 MinorIter  197
CPU time: 0.078000 sec. Elapsed time: 0.070000 sec. 

Problem type appears to be: lpcon
Time for symbolic processing: 0.46746 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Oil Pyrolysis                  f_k      -0.354368283541904860
                                       sum(|constr|)      0.000000002719365042
                              f(x_k) + sum(|constr|)     -0.354368280822539790
                                              f(x_0)     -0.357327805323273070

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

FuncEv    1 ConstrEv  208 ConJacEv  208 Iter  112 MinorIter  305
CPU time: 0.171601 sec. Elapsed time: 0.177000 sec. 

Problem type appears to be: lpcon
Time for symbolic processing: 0.47201 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Oil Pyrolysis                  f_k      -0.351747594492437030
                                       sum(|constr|)      0.000000280199856233
                              f(x_k) + sum(|constr|)     -0.351747314292580830
                                              f(x_0)     -0.354368283541905970

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

FuncEv    1 ConstrEv  124 ConJacEv  124 Iter   71 MinorIter  281
CPU time: 0.156001 sec. Elapsed time: 0.154000 sec. 

Problem type appears to be: lpcon
Time for symbolic processing: 0.46076 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Oil Pyrolysis                  f_k      -0.352833701465704920
                                       sum(|constr|)      0.000000018165177960
                              f(x_k) + sum(|constr|)     -0.352833683300526950
                                              f(x_0)     -0.351747594492436640

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

FuncEv    1 ConstrEv  409 ConJacEv  409 Iter  194 MinorIter  697
CPU time: 0.686404 sec. Elapsed time: 0.685000 sec. 

Problem type appears to be: lpcon
Time for symbolic processing: 0.48836 seconds
Starting numeric solver
===== * * * =================================================================== * * *
TOMLAB - TOMLAB Development license  999007. Valid to 2011-12-31
=====================================================================================
Problem: ---  1: Oil Pyrolysis                  f_k      -0.352618526247765740
                                       sum(|constr|)      0.000016167955157926
                              f(x_k) + sum(|constr|)     -0.352602358292607830
                                              f(x_0)     -0.352833701465704590

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

FuncEv    1 ConstrEv   64 ConJacEv   64 Iter   46 MinorIter  364
CPU time: 0.187201 sec. Elapsed time: 0.191000 sec. 

end

t  = subs(collocate(t),solution);
x1 = subs(collocate(x1opt),solution);
x2 = subs(collocate(x2opt),solution);
x3 = subs(collocate(x3opt),solution);
x4 = subs(collocate(x4opt),solution);
T  = subs(collocate(Topt),solution);

Plot result

subplot(2,1,1)
plot(t,x1,'*-',t,x2,'*-',t,x3,'*-',t,x4,'*-');
legend('x1','x2','x3','x4');
title('Oil Pyrolysis state variables');

subplot(2,1,2)
plot(t,T,'+-');
legend('T');
title('Oil Pyrolysis control');