1 Introduction

The paper [1] by Jackson presented the analytical solution to the plug flow reactor catalyst blend problem, a classic dynamic optimization example proposed by Gunn and Thomas [2], and thereafter the problem was studied by many researchers [3–6]. Among the studied cases of this problem, the most common ones were the reactor length equals 1 and 12. However, neither the explicit expression nor any numerical value of the optimal objective function was given in [1], which has caused confusion and inaccuracies in the subsequent literatures [7–11]. For example, Dadebo and McAuley [7] supposed that the optimal value of this problem with the reactor length equals 12 is 0.476946 according to [1], and, for the case that reactor length equals 1, Bell and Sargent [12] achieved the best literature result of 0.048080. However, using the time-scaling method [13], the authors derived different results from these reported ones with errors larger than their respective optimization tolerances. The obtained results, such as a better result of 0.477712 for the reactor length of 1, and 0.048055 for the reactor length of 12, are reported in detail in Sect. 4.

Since the exact optimal value of this problem has not been reported, these inaccuracies in the solutions can hardly be recognized and judged by the researchers. Thus, we present the exact optimal solutions and optimal objective values analytically in this study as a supplement and correction for the cases when the reactor length equals 1 and 12. Moreover, these results are verified numerically by using the time-scaling method.

2 Problem Formulation

Chemical reactions of economic importance often take place in several steps through a number of intermediate products, and each stage may be catalyzed by a different catalytic substance [14]. However, under certain circumstances, there are advantages in mixing the catalysts in a single reaction vessel [1]. Gunn and Thomas quoted the example of the reaction scheme

$$\begin{aligned} A \underset{2}{\overset{1}{\rightleftharpoons }} B\buildrel 3 \over \longrightarrow C. \end{aligned}$$
(1)

They remarked that there is an optimum catalyst blend, which gives the highest yield of \(C\) for a given reactor length, and further improvement can be obtained by varying the catalyst blend along the reactor [2].

Referring to reaction scheme (1), we denote the unknown mole fractions of the substances \(A\) and \(B\) in the mixture by \(x_A \) and \(x_B \), respectively, and assume that all the reactions are carried out in an isothermal tubular reactor. Then, the optimization problem of maximizing the concentration of the final product \(C\) can be formulated as follows:

$$\begin{aligned} \begin{array}{ll} \hbox {max}&{}[J=1-x_A (t_f )-x_B (t_f )] \\ \hbox {s.t. }&{}\frac{\hbox {d}x_A }{\hbox {d}t}=u(k_2 x_B -k_1 x_A ), \\ &{}\frac{\hbox {d}x_B }{\hbox {d}t}=u(k_1 x_A -k_2 x_B )-(1-u)k_3 x_B ,\\ &{}x_A (0)=1, x_B (0)=0, \\ &{}0\le t\le t_f , 0\le u\le 1, \end{array} \end{aligned}$$
(2)

where \(t_f \) is the fixed reactor length and \(t\) denotes the distance from the entry of the reactor. The symbols \(k_1 \), \(k_2 \), and \(k_3\) are, respectively, the velocity constants of reactions 1, 2, and 3 in the reactor. The catalyst blend \(u\) denotes the fraction of the catalyst, which catalyzes the reactions \(A\rightleftharpoons B\). This fraction can be varied as required along the reactor by suitable mixing of the two catalysts of the reaction scheme (1). When the blend has the value \(u\), the effective velocity constants are \(uk_1 \), \(uk_2 \), and \((1-u)k_3 \), as indicated in the formulation (2).

3 Analytical Solution

3.1 Structure of the Optimal Solution

Since the objective function \(J\) and the right-hand side of the state differential equations in (2) are Lipschitz continuously differentiable functions of each of their arguments, Pontryagin’s maximum principle can be employed to solve the problem stated above. According to the maximum principle, if the control variable \(u^*(t)\) and its associated state trajectories are optimal, then there exist adjoint variables \(\lambda _A \) and \(\lambda _B \), such that the adjoint equations

$$\begin{aligned} \frac{\hbox {d}\lambda _A }{\hbox {d}t}=-\frac{\partial H}{\partial x_A}, \frac{\hbox {d}\lambda _B }{\hbox {d}t}=-\frac{\partial H}{\partial x_B } \end{aligned}$$
(3)

hold on \([0,t_f ]\) with terminal condition

$$\begin{aligned} \lambda _A (t_f )=\lambda _B (t_f )=-1, \end{aligned}$$
(4)

where

$$\begin{aligned} H=\lambda _A u(k_2 x_B -k_1 x_A )+\lambda _B [u(k_1 x_A -k_2 x_B )-(1-u)k_3 x_B ] \end{aligned}$$
(5)

is the Hamiltonian, and there holds

$$\begin{aligned} H(\lambda _A (t),\lambda _B (t),x_A (t),x_B (t),u^*(t))=\mathop {\max }\limits _{0\le u\le 1} H(\lambda _A (t),\lambda _B (t),x_A (t),x_B (t),u(t))\nonumber \\ \end{aligned}$$
(6)

on \([0,t_f ]\). Additionally, for the Hamiltonian of this problem, there exists the following result.

Proposition 3.1

The maximum of Hamiltonian

$$\begin{aligned} H_{\max } =\mathop {\max }\limits _{0\le u\le 1} H(\lambda _A (t),\lambda _B (t),x_A (t),x_B (t),u(t)) \end{aligned}$$

with respect to the control variable \(u\) is a constant over \([0,t_f ]\) in this problem.

Proof

The derivative of \(H_{\max } \) with respect to \(t\) is

$$\begin{aligned} \frac{\hbox {d}H_{\max } }{\hbox {d}t}&= \frac{\partial H_{\max } }{\partial t}+\frac{\partial H_{\max } }{\partial \lambda _A }\frac{\partial \lambda _A }{\partial t}+\frac{\partial H_{\max } }{\partial \lambda _B }\frac{\partial \lambda _B }{\partial t}+\frac{\partial H_{\max } }{\partial x_A }\frac{\partial x_A }{\partial t}\\&+\,\frac{\partial H_{\max } }{\partial x_B }\frac{\partial x_B }{\partial t}+\frac{\partial H_{\max } }{\partial u}\frac{\partial u}{\partial t}. \end{aligned}$$

When \(H \) is maximized with respect to the control variable, it holds either \(\frac{\partial H_{\max } }{\partial u}=0\), if \(u\) is in the internal of the control region, or \(\frac{\partial u}{\partial t}=0\), if \(u\) achieves the bound. Thus

$$\begin{aligned} \frac{\partial H_{\max } }{\partial u}\frac{\partial u}{\partial t}=0. \end{aligned}$$

Moreover, based on (3) and (5), the followings can be obtained:

$$\begin{aligned} \frac{\partial H_{\max } }{\partial \lambda _A }\frac{\partial \lambda _A }{\partial t}+\frac{\partial H_{\max } }{\partial x_A }\frac{\partial x_A }{\partial t}=0 \end{aligned}$$

and

$$\begin{aligned} \frac{\partial H_{\max } }{\partial \lambda _B }\frac{\partial \lambda _B }{\partial t}+\frac{\partial H_{\max } }{\partial x_B }\frac{\partial x_B }{\partial t}=0. \end{aligned}$$

So, the expression of \(\frac{\hbox {d}H_{\max } }{\hbox {d}t}\) can be simplified as

$$\begin{aligned} \frac{\hbox {d}H_{\max } }{\hbox {d}t}=\frac{\partial H_{\max } }{\partial t}. \end{aligned}$$

Because the Hamiltonian \(H\) does not depend explicitly on time, \(\frac{\partial H_{\max } }{\partial t}=0\); then, we have \(\frac{\hbox {d}H_{\max } }{\hbox {d}t}=0\). Thus, the value of \(H_{\max } \) is a constant in this problem.

The adjoint equations are

$$\begin{aligned} \frac{\hbox {d}\lambda _A}{\hbox {d}t}&= -uk_1 (\lambda _B -\lambda _A),\end{aligned}$$
(7)
$$\begin{aligned} \frac{\hbox {d}\lambda _B}{\hbox {d}t}&= uk_2 (\lambda _B -\lambda _A)+(1-u)k_3 \lambda _B. \end{aligned}$$
(8)

Then, according to (6), a necessary condition for the optimality of \(u(t)\) is that the Hamiltonian

$$\begin{aligned} H=\rho u-\lambda _B k_3 x_B , \rho =(\lambda _B -\lambda _A )(k_1 x_A -k_2 x_B )+\lambda _B k_3 x_B \end{aligned}$$
(9)

has its greatest value with respect to \(u\) for each \(t\). Since \(H\) is linear with \(u\), this implies that \(u=0\) or \(u=1\), depending on the sign of \(\rho \) at the point in question. The catalyst blend may take a value between these bounds only if \(\rho \) vanishes. When this occurs over a finite interval of \(t\), then the corresponding part of the solution is called a singular segment.

Taking into account of (4), the Hamiltonian at \(t=t_f \) is given by \(H=-uk_3 x_B +k_3 x_B \) and is maximized by taking \(u=0\). Thus, the optimal catalyst blend starts back from \(t=t_f \) with \(u=0\) and retains this value until \(\rho \) changes its sign. Moreover, Jackson pointed out that the condition \(\rho = 0\) cannot be satisfied again after switched to \(u=\hbox {1}\) [1]. At this point, we obtain the following result about the structure of the optimal solution.

Proposition 3.2

There are four possible control structures for the optimal solution:

  1. 1.

    \(u=0\) and \(\rho <0\) hold over the whole interval \([0,t_f ]\).

  2. 2.

    The control starts from \(t=0\) with a singular segment (\(\rho =0)\) and then switches to \(u=0\) with \(\rho <0\).

  3. 3.

    The control starts from \(t=0\) with \(u=1\) and \(\rho >0\) and then switches to \(u=0\) with \(\rho <0\).

  4. 4.

    The control starts from \(t=0\) with \(u=1\) and \(\rho >0\) and then, after going through a singular segment (\(\rho =0)\), switches to \(u=0\) with \(\rho <0\).

For the first situation, the control variable \(u=0\) over the whole interval \([0,t_f ]\). It will make the effective velocity constants of reactions 1 and 2 as well as the concentration of the final product \(C\) to be 0, which is unacceptable in the real production process. Thus, we consider the last three situations.

3.2 Properties for a Singular Segment

Suppose that \(u\) takes values between its bounds for the finite time interval \(t_1 <t<t_2 \), it is clearly necessary that

$$\begin{aligned} \rho =0,\quad t_1 <t<t_2. \end{aligned}$$
(10)

In turn, this implies that \(\frac{\hbox {d}\rho }{\hbox {d}t}=0\) for all \(t_1 <t<t_2 \). Using the state and adjoint equations, this condition can be reduced to the form

$$\begin{aligned} \frac{\hbox {d}\rho }{\hbox {d}t}=k_3 (\lambda _B k_1 x_A -\lambda _A k_2 x_B )=0,\quad t_1 <t<t_2. \end{aligned}$$
(11)

Solving (11) for \(\lambda _A \) and substituting \(\lambda _A =\lambda _B k_1 x_A /k_2 x_B \) into (10), the following necessary condition for the singular segment is obtained:

$$\begin{aligned} (\lambda _B /k_2 x_B )[k_2 k_3 x_B^2 -(k_1 x_A -k_2 x_B )^2]=0,\quad t_1 <t<t_2 , \end{aligned}$$
(12)

hence

$$\begin{aligned} \lambda _B =0 \end{aligned}$$
(13)

or

$$\begin{aligned} k_2 k_3 x_B^2 =(k_1 x_A -k_2 x_B )^2. \end{aligned}$$
(14)

From (12), it follows that (13) implies \(\lambda _A =\lambda _B k_1 x_A /k_2 x_B =0\), and hence \(H=0\), at all the points of the singular segment. However, we have already found the form of \(H\) at \(t=t_f \) and seen that it is maximized by \(u=0\). The corresponding maximum value of \(H\) is \(H_{\max } (t_f )=k_3 x_B >0\). Thus, a singular segment on which \(H=0\) everywhere cannot belong to the optimal solution, since it is known from Proposition 3.1 that \(H_{\max } \) must remain constant throughout the entire solution. Therefore, the possibility (13) can be discarded, and we only need to consider (14), which reduces to

$$\begin{aligned} k_1 x_A =k_2 x_B (1\pm \sqrt{k_3 /k_2 } ). \end{aligned}$$
(15)

Then, from (11), the corresponding relationship between the adjoint variables is

$$\begin{aligned} \lambda _A /\lambda _B =1\pm \sqrt{k_3 /k_2 }. \end{aligned}$$
(16)

The singular segment must also be a solution of state equations for a suitable choice of \(u\). To determine the necessary form of \(u(t)\), we obtain the following differential equation from (15)

$$\begin{aligned} \frac{\hbox {d}x_B }{\hbox {d}x_A }=\frac{k_1 }{k_2 (1\pm \sqrt{k_3 /k_2 } )} \end{aligned}$$
(17)

and compare it with the corresponding differential equation obtained from the state equations, namely,

$$\begin{aligned} \frac{\hbox {d}x_B }{\hbox {d}x_A }=\frac{u(k_1 x_A -k_2 x_B )-(1-u)k_3 x_B }{u(k_2 x_B -k_1 x_A )}. \end{aligned}$$
(18)

It can be seen that (18) reduces to (17) if and only if,

$$\begin{aligned} \frac{k_1 }{k_2 (1\pm \sqrt{k_3 /k_2 } )}=\frac{\pm u\sqrt{k_3 /k_2 } -(1-u)k_3 /k_2 }{\mp u\sqrt{k_3 /k_2 } }. \end{aligned}$$
(19)

Define the quantities

$$\begin{aligned} \alpha :=\sqrt{k_3 /k_2 } , \beta :=k_1 /k_2 , \end{aligned}$$
(20)

in order to get a simple expression for \(u\). Taking the upper signs in (19), we then obtain the relation

$$\begin{aligned} u=\alpha (1+\alpha )/[\beta +(1+\alpha )^2] \end{aligned}$$
(21)

and, correspondingly, the followings are obtained from (15) and (16):

$$\begin{aligned} \lambda&:= \lambda _A /\lambda _B =1+\alpha ,\end{aligned}$$
(22)
$$\begin{aligned} z&:= x_B /x_A =\beta /(1+\alpha ). \end{aligned}$$
(23)

Similarly, after taking the lower signs in (19), we deduce that

$$\begin{aligned} u=-\alpha (1-\alpha )/[\beta +(1-\alpha )^2], \end{aligned}$$
(24)

together with

$$\begin{aligned} \lambda :=\lambda _A /\lambda _B =1-\alpha \end{aligned}$$
(25)

and

$$\begin{aligned} z:=x_B /x_A =\beta /(1-\alpha ). \end{aligned}$$
(26)

If \(\alpha <1\), (24) gives a negative value for \(u\), which is physically unacceptable, while, if \(\alpha >1\), (26) gives a negative value for \(z\), which is also physically unacceptable. Thus, Eqs. (21)–(23) describe the properties of the only acceptable singular segment for this problem. From (21), it is easy to see that \(u<1\) for all \(\alpha >0\) and \(\beta >0\); therefore, this value of \(u\) lies in the permitted interval, whatever the reaction kinetics.

3.3 Solution with a Singular Segment

Based on the properties of the singular segment, we consider the second and fourth situations of Proposition 3.2. Taking into account of the second situation, since \(x_A (0)=1\) and \(x_B (0)=0\), there exists \(\rho (0)=k_1 (\lambda _B -\lambda _A )=0\). Then, \(\lambda _A =\lambda _B \) can be obtained, which is in contradiction with the Eq. (22). Thus, the second situation can be discarded, leaving the fourth situation.

From (7) and (8), we obtain a differential equation for \(\lambda \), namely,

$$\begin{aligned} \frac{\hbox {d}\lambda }{\hbox {d}t}=-u(1-\lambda )(k_1 +\lambda k_2)-(1-u)\lambda k_3 , \end{aligned}$$
(27)

which, on the segment \(u=0\), reduces to

$$\begin{aligned} \frac{\hbox {d}\lambda }{\hbox {d}t}=-\lambda k_3. \end{aligned}$$
(28)

The segment \(u=0\) has the terminal abscissas \(t=t_2 \) and \(t=t_f \), and we know that \(\lambda (t_f )=1\) from (4). Thus, the integration of (28) leads to

$$\begin{aligned} \lambda (t_2 )=\exp [k_3 (t_f -t_2 )]. \end{aligned}$$
(29)

The above relation and (22) permit one to express \(t_f -t_2\) as follows:

$$\begin{aligned} k_3 (t_f -t_2 )=\log (1+\alpha ). \end{aligned}$$
(30)

From the differential state equations, it can be seen that

$$\begin{aligned} \frac{\hbox {d}z}{\hbox {d}t}=u(1+z)(k_1 -k_2 z)-(1-u)k_3 z, \end{aligned}$$
(31)

which, on the segment \(u=1\), reduces to

$$\begin{aligned} \frac{\hbox {d}z}{\hbox {d}t}=(1+z)(k_1 -k_2 z). \end{aligned}$$
(32)

Since \(z=0\) at \(t=0\), integrating (32) over the range \([0,t_1 ]\) gives

$$\begin{aligned} k_3 t_1 =[\alpha ^2/(1+\beta )]\log \{\beta [1+z(t_1 )]/[\beta -z(t_1 )]\}. \end{aligned}$$
(33)

Because there exists (23) in the singular segment, we obtain

$$\begin{aligned} k_3 t_1 =[\alpha ^2/(1+\beta )]\log [(1+\alpha +\beta )/\alpha ]. \end{aligned}$$
(34)

On the singular segment, the expression of \(\frac{\hbox {d}x_B}{\hbox {d}t}\) in (2) can be reduced to

$$\begin{aligned} \frac{\hbox {d}x_B }{\hbox {d}t}=\frac{-k_2 \beta x_B }{\beta +(1+\alpha )^2}, \end{aligned}$$
(35)

and this can be integrated to obtain the following relation:

$$\begin{aligned} \log [x_B (t_1 )/x_B (t_2 )]=k_2 \beta (t_2 -t_1 )/[\beta +(1+\alpha )^2]. \end{aligned}$$
(36)

On the segment \(u=1\), the effective velocity of reaction 3 is zero, so \(x_A +x_B =1\) at any point in the range of \(t\le t_1 \). Moreover, from (23), the value of \(x_B (t_1 )\) can be obtained as follows:

$$\begin{aligned} x_B (t_1 )=\beta /(1+\alpha +\beta ), \end{aligned}$$
(37)

so that (36) is reduced to

$$\begin{aligned} k_3 (t_2 -t_1 )=[1+(1+\alpha )^2/\beta ]\log [\beta /x_B (t_2 )(1+\alpha +\beta )]. \end{aligned}$$
(38)

Then, adding (30), (33), and (38) gives the following expression of \(t_f \):

$$\begin{aligned} k_3 t_f&= \log (1+\alpha )+[\alpha ^2/(1+\beta )]\log [(1+\alpha +\beta )/\alpha ] \nonumber \\&+\,[1+(1+\alpha )^2/\beta ]\log [\beta /x_B (t_2 )(1+\alpha +\beta )], \end{aligned}$$
(39)

and the value of \(x_B (t_2 )\) can be obtained as

(40)

From the equation of the singular segment, the following can be deduced:

$$\begin{aligned} x_A (t_f )=x_A (t_2 )=x_B (t_2 )(1+\alpha )/\beta . \end{aligned}$$
(41)

By integrating the differential state equation of \(x_B \) along the final segment \(u=0\), the following can be obtained:

$$\begin{aligned} x_B (t_f )=x_B (t_2 )\exp [-k_3 (t_f -t_2 )]=x_B (t_2 )/(1+\alpha ). \end{aligned}$$
(42)

Finally, from (41) and (42), \(J\) can be expressed in terms of \(x_B (t_2 )\):

$$\begin{aligned} J=1-x_A (t_f )-x_B (t_f )=1-x_B (t_2 )[1/(1+\alpha )+(1+\alpha )/\beta ]. \end{aligned}$$
(43)

Note that \(t_1 <t_2 \) for the singular segment. Thus, according to (30) and (34), the derived solution with a singular segment holds when

$$\begin{aligned} t_f >[\log (1+\alpha )]/k_3 +[\alpha ^2/k_3 (1+\beta )]\log [(1+\alpha +\beta )/\alpha ] \end{aligned}$$
(44)

3.4 Bang–Bang Solution

Now, consider the third situation of Proposition 3.2. Let \(t_s \) be the value of \(t\) at which the switch from \(u=1\) to \(u=0\) is made. Then, \(\rho (t_s )=0\), and from (9), this implies that

$$\begin{aligned} \lambda (t_s )=1+\alpha ^2z(t_s )/[\beta -z(t_s )]. \end{aligned}$$
(45)

Meanwhile, similar to (29), the following equation can be obtained:

$$\begin{aligned} \lambda (t_s )=\exp [k_3 (t_f -t_s )]. \end{aligned}$$
(46)

Based on (45) and (46), \(t_f -t_s \) can be expressed in the term of \(z(t_s )\) as

$$\begin{aligned} k_3 (t_f -t_s )=\log \{1+\alpha ^2z(t_s )/[\beta -z(t_s )]\}. \end{aligned}$$
(47)

And, according to the derivation of (33), a second expression relating \(t_s \) and \(z(t_s )\) can be obtained:

$$\begin{aligned} k_3 t_s =[\alpha ^2/(1+\beta )]\log \{\beta [1+z(t_s )]/[\beta -z(t_s )]\}. \end{aligned}$$
(48)

By adding (47) and (48), \(t_f \) can be expressed as follows:

$$\begin{aligned} k_3 t_f&= [\alpha ^2/(1+\beta )]\log \{\beta [1+z(t_s )]/[\beta -z(t_s )]\}\nonumber \\&+\,\log \{1+\alpha ^2z(t_s )/[\beta -z(t_s )]\}. \end{aligned}$$
(49)

Finally, we calculate \(x_A (t_f )\), \(x_B (t_f )\), and hence \(J\), as functions of \(z(t_s )\). At \(t=t_s \),

$$\begin{aligned} x_A (t_s )+x_B (t_s )=1, \end{aligned}$$

whence

$$\begin{aligned} x_A (t_s )=1/[1+z(t_s )], x_B (t_s )=z(t_s )/[1+z(t_s )]. \end{aligned}$$
(50)

Since \(x_A \) does not change on the terminal segment \(u=0\),

$$\begin{aligned} x_A (t_f )=x_A (t_s )=1/[1+z(t_s )]. \end{aligned}$$
(51)

Then, relate \(x_B (t_f )\) to \(x_B (t_s )\) by integrating the differential state equation of \(x_B \) after setting \(u=0\). It gives that

$$\begin{aligned} x_B (t_f )=x_B (t_s )\exp [-k_3 (t_f -t_s )]. \end{aligned}$$

Because of (47) and (50), this reduces to

$$\begin{aligned} x_B (t_f )=\{z(t_s )/[1+z(t_s )]\}/\{1+\alpha ^2z(t_s )/[\beta -z(t_s )]\}. \end{aligned}$$
(52)

Thus, the objective function becomes

$$\begin{aligned} J&= 1-x_A (t_f )-x_B (t_f )=1-1/[1+z(t_s )]\nonumber \\&-\,\{z(t_s )/[1+z(t_s )]\}/\{1+\alpha ^2z(t_s )/[\beta -z(t_s )]\}. \end{aligned}$$
(53)

Equations (49) and (53) provide the relationship between \(J\) and \(t_f \) in parametric form, with \(z(t_s )\) as the parameter. Similarly, Eqs. (48) and (49) provide a parametric relationship between \(t_s \) and \(t_f \).

3.5 Global Optimum

The four possible control structures for the optimal solution have been discussed above. We can find that, if the reactor is sufficiently small, corresponding to a small value of \(t_f \) such that

$$\begin{aligned} t_f \le [\log (1+\alpha )]/k_3 +[\alpha ^2/k_3 (1+\beta )]\log [(1+\alpha +\beta )/\alpha ], \end{aligned}$$
(54)

there will be no solution involving the singular segment. Therefore, for the value of \(t_f \) that satisfies (54), the optimal blending policy is the bang–bang policy described in Sect. 3.4. On the other hand, for larger values of \(t_f \), there are two alternative policies: the policy with a singular segment and the bang–bang policy. Among these policies, the one giving greater value of the objective function can be found only by direct comparison of the results, since both policies satisfy the optimality necessary condition.

Based on the obtained conclusions, we are already to solve the cases when the reactor length equals 1 and 12. Considering the values of the kinetic constants assumed by Gunn and Thomas, namely, \(k_1 =k_3 =1,k_2 =10\), which are also the same as those used in [7], the high-accuracy global optimum of reactor length \(t_f =12\) is calculated first.

From (20), the value of \(\alpha \) and \(\beta \) can be obtained easily:

$$\begin{aligned} \alpha ={\sqrt{10} } / {10}, \quad \beta =1/10. \end{aligned}$$

As \(t_f =12\) satisfies Eq. (44), both policies should be computed for comparison. Using (30) and (34), we find that the policy with a singular segment requires \(u=1\) for an interval of length

$$\begin{aligned} t_1 =\frac{\log (1+1.1\sqrt{10} )}{11}\approx \hbox {0.136299034594555} \end{aligned}$$

at the entry to the reactor and \(u=0\) from the length

$$\begin{aligned} t_2 =12-\log (1+{\sqrt{10} } / {10})\approx \hbox {11.725230107591655} \end{aligned}$$

to the reactor exit. Throughout the remainder of the reactor length \(t_1 \le t\le t_2 \), the blend takes the intermediate value

$$\begin{aligned} u=\frac{\hbox {5}\sqrt{10} -4}{\hbox {52}}\approx 0.227142082708498, \end{aligned}$$

computed from (21). Then, from (37), the value of \(x_B (t_2 )\) can be obtained as

$$\begin{aligned} x_B (t_2 )&= \exp \left[ \frac{\log (14.552022\sqrt{10} +45.979671)-132}{22\sqrt{10} +132}\right] \Big /\sqrt{10} +11 \\&\approx 0.037515231832340 \\ \end{aligned}$$

From (43), the following is obtained:

$$\begin{aligned} J(12)\approx 0.477712020050041. \end{aligned}$$

For the bang–bang policy, the objective function value can be likewise calculated according to the derived equations in Sect. 3.4. From (49), the value of \(z(t_s )\) can be obtained as

$$\begin{aligned} z(t_s )\approx 0.0999999989. \end{aligned}$$

Then, from (53), the objective function value of the bang–bang policy is obtained

$$\begin{aligned} J(12)\approx 0.09090908, \end{aligned}$$

which is less than the value of the policy with a singular segment. Thus, the optimum blend policy for \(t_f =12\) is the policy with a singular segment, which is illustrated in Fig. 1. This is consistent with the result obtained by England [15]. The optimal objective function value is \(J^*(12)\approx 0.477712020050041\), rather than the value 0.476946 obtained in [7].

Fig. 1
figure 1

Optimal blend policy for \(t_f =12\)

Full size image

Similarly, the optimal solution and optimal objective function value can be calculated for the case when the reactor length \(t_f =1\). Again, the optimal policy is the policy with a singular segment, where

$$\begin{aligned} t_1&= \frac{\log (1+1.1\sqrt{10} )}{11}\approx 0.136299034594555,\\ t_2&= 1-\log (1+{\sqrt{10} } / {10})\approx 0.725230107591655, \end{aligned}$$

and on the singular segment

$$\begin{aligned} u=\frac{\hbox {5}\sqrt{10} -4}{52}\approx 0.227142082708498. \end{aligned}$$

Correspondingly, \(J^*(1)\approx 0.048055685860877\), rather than the best literature result of 0.048080 under the optimization tolerance \(10^{-6}\). The optimum blend policy for \(t_f =1\) is shown in Fig. 2 as an illustration.

Fig. 2
figure 2

Optimal blend policy for \(t_f =1\)

Full size image

4 Numerical Verification

To verify the derived results, the problem with kinetic constants \(k_1 =k_3 =1\) and \(k_2 =10\), which were assumed by Gunn and Thomas [2], and also used in [3–13], is solved numerically under different tolerances using the time-scaling method [13], since it can optimize both the control parameters and the switching times effectively.

The time-scaling method was proposed by Teo et al. [16]. It maps the variable switching times to a new time scale as pre-fixed time nodes and treats the lengths of partitions as variables to be optimized. Denote the control and length parameters of the \(i\)-th segment by \(\sigma _i \) and \(\theta _i \); then, the gradients of the objective function \(J\) with respect to \(\sigma _i \) and \(\theta _i\), respectively, can be obtained by

$$\begin{aligned} \frac{\partial J}{\partial \sigma _i }&= \int _0^{t_f } {\frac{\partial H}{\partial \sigma _i }} dt,\\ \frac{\partial J}{\partial \theta _i }&= \int _0^{t_f } {\frac{\partial H}{\partial \theta _i }} dt. \end{aligned}$$

On this basis, both the control and length parameters can be optimized easily as a nonlinear programming problem by gradient-based optimization technique. For more details of the time-scaling method, please refer to [13, 17].

Here, the ODE45 function of MATLAB is used to solve the relevant differential equations, and the SQP method implemented in the FMINCON function is used for optimization. The optimization is started with the whole reactor length divided into 3 segments with equal size, and the initial guess of the control is set to 0.5. All these initial values of the parameters are set without using any prior knowledge about the optimal solution. The obtained results listed in Tables 1 and 2 indicate that, as the tolerances of ODE and NLP solver increase, the control parameters and switching times get closer to their corresponding analytical values, and the objective function value \(J\) gradually converges to the value of \(J^*(t_f )\).

Table 1 The solution of time scaling for \(t_f =12\)
Full size table
Table 2 The solution of time scaling for \(t_f =1\)
Full size table

5 Conclusions

The analytical solutions of the plug flow reactor catalyst blend problem with the reactor length \(t_f =1\) and 12 are presented in this study. Correspondingly, the correct high-accuracy global optimums are given, which are verified numerically using the time-scaling method. It can be found that the best literature result of 0.048080 for \(t_f =1\) and the supposed optimal value 0.476946 for \(t_f =12\) are inaccurate. The exact optimal values for \(t_f =1\) and 12 should be, respectively, \(0.048055685860877\) and \(0.477712020050041\).