Black-Box Optimization of Operating Condition in a Chemical Reactor¶

As a practical example of black-box optimization (we may call it FMQA), we will work on optimization using numerical simulations of nonlinear physical phenomena as the objective function.

In this tutorial, the goal is to optimize the operating condition of a plant reactor to maximize production. The chemical reaction and transport phenomena of reactants and products in the reactor are predicted by numerical simulations using the finite difference method. Based on the predicted results, optimization is performed to maximize the amount of material production in the reactor.

Note that the simulation code used in this example is implemented in Python and can be called directly from the FMQA class which is introduced later. However, even if this is not the case (e.g. if the machines for FMQA execution and simulation execution are different or the objective function is based on experimental measurements), this example code can be used almost as is.

For an introduction to the black-box optimization and FMQA, see "Black-Box Optimization Exploration of Model Superconducting Materials". For another advanced application case using FMQA, see "Black-Box Optimization of Airfoil Geometry by Fluid Flow Simulation".

This notebook contains the following sections.

• 1. Problem setting
  • 1.1. Reactor model and optimization target
  • 1.2. Description of the reactor simulator
  • 1.3. Implementation of the reactor simulator
  • 1.4. Simulation examples and definition of the objective function
• 2. FMQA program implementation
  • 2.1. Random seed initialization
  • 2.2. Configuration of Amplify client
  • 2.3. Implementing FM with PyTorch
  • 2.4. Construction of initial training data
  • 2.5. Execution class for FMQA cycle
• 3. Search for optimal operating conditions
  • 3.1. FMQA execution example
  • 3.2. Transition of objective function values during the optimization process
  • 3.3. Example output from this FMQA sample program

1. Problem setting¶

This section describes the problem definition in this example code and the simulation of the reactor used as the objective function. However, black-box optimization treats the objective function as a black box, so it is unnecessary to understand this simulation code.

1.1. Reactor model and optimization target¶

As shown in the figure above, we consider a chemical reactor controlled solely by the initial concentration distribution of the reactive substance A. In this reactor, the chemical reaction $A \rightarrow B$ occurs at a reaction rate based on the concentration distributions of A and B, producing the product B.

The goal of the present optimization is to maximize the total amount of B produced within a given production time $\Delta t_{prod}$ by appropriately determining the initial concentration distribution of A to improve productivity (think of it as a problem of minimizing the negative value of the total amount).

1.2. Description of the reactor simulator¶

This time, instead of considering reactors in actual three-dimensional (3D) space, we will focus on one-dimensional (1D) reactors for the sake of simulation cost. Such a reactor can be seen as, for example,

• A chemical reactor with a long and thin vessel, or
• A chemical reactor homogeneous in the y- and z-directions

Note that the essence of the present exercise, i.e., optimization for the reactor operating conditions, is independent of the dimension considered.

The chemical reaction and transport phenomena in the reactor are described by the following nonlinear partial differential equations, where $C_A$ denotes the concentration of substance A (reactant) and $C_B$ denotes the concentration of substance B (product).

$$\frac{\partial C_A}{\partial t} = \alpha \frac{\partial}{\partial x}\left(\frac{\partial C_A}{\partial x}\right) - \omega$$

$$\frac{\partial C_B}{\partial t} = \alpha \frac{\partial}{\partial x}\left(\frac{\partial C_B}{\partial x}\right) + \omega$$

$$\omega = R_r C_A (1-C_A) \exp{(-C_B)}$$

$$C_B = O \text{ at } t=0 \:\:\text{(initially, there is no B)}$$

Given an initial concentration distribution of A, $C_{A,0}$, the total production of B within the production time $\Delta t_{prod}$ corresponds to the spatial integral of $C_B$ in the reactor at $\Delta t_{prod}$ from the beginning of the simulation and can be obtained with the simulation class Reactor described next.

1.3. Implementation of the reactor simulator¶

The integrate function of the Reactor class simulates reactions in the chemical reactor. This simulator solves the above equations using the finite difference method, similar to commercially available physical simulation software using finite element or finite volume methods.

Here, the integrate function takes an initial concentration distribution of A, $C_{A,0}$, as an argument and returns the total amount of B produced (the spatial integral of the concentration distribution of B) obtained after the simulation is run for the physical time of $\Delta t_{prod}$.

1.4. Simulation examples and definition of the objective function¶

Now, let's execute the reaction simulator using the integrate function of the Reactor class. The function returns the total amount of B produced during the production time, given an initial concentration distribution $C_{A,0}$ of A (which is set by a random number for now). The first argument of the integrate function is a 1D binary array representing the distribution of $C_{A,0}$ in 1D space (i.e., $C_{A,0}$ takes either 0 or 1 in each coordinate). The second argument is an optional output flag for the result image (False by default).

Upon execution, the following result image is obtained according to $C_{A,0}$, which is determined based on a random number.

The resulting image above shows the initial concentration distribution of A $C_{A,0}$ (gray) and the concentration distribution of B at each time $C_B$ (red). The concentration distribution of B is shown at time $t=0$ (initial concentration distribution, bottom red line at $C_{B,0}=0$), $t=0.1 \Delta t_{prod}$ (second red line from the bottom), $t= 0.2 \Delta t_{prod}$ (third red line from the bottom), $t=0.4 \Delta t_{prod}$ (fourth red line from the bottom), $t = \Delta t_{prod}$ (production end time, top bold red line). In this optimization exercise, the goal is to maximize the integral value of the bold red line in the reactor.

At time 0, $C_B=C_{B,0}=0$ in the whole region, but with time, the chemical reaction proceeds, and B is produced ($C_B$ increases). The final concentration distribution of B at the end of the production time is shown by the bold red line. The time evolution of reactant A is not shown, but it is assumed that $C_{A}$ gradually decreases over time due to chemical reactions and molecular diffusion.

Since the random seed value is not fixed in the example code below, $C_{A,0}$ changes with each run, and the total production of B changes accordingly. Can you imagine what the initial concentration distribution $C_{A,0}$ of A that maximizes the total production of B would look like?

When the reactor domain is discretized with 100 finite volume, as in the above code, there are $2^{100} \sim 10^{30}$ possible values that the initial concentration distribution $C_{A,0}$ vector of A can take. Also, from the reaction rate equation $\omega$ in section 1.2, in order to maximize the production of B in time, it is not simply a matter of filling the entire reactor with A (i.e. $C_{A,0} =$[1, 1, 1,.... , 1]), but it is necessary to devise a way to fill it with as much A as possible, while appropriately arranging the local region where $C_A=0$. For such a system,

• The search space is large, and a full search is unrealistic due to the time cost of the simulation,
• The objective function is an unknown function (described by a nonlinear partial differential equation) and is a black box.

Therefore, the use of FMQA is considered as effective.

Although the material cost depends on the amount of A initially present in the reactor, we assume that the cost of reactant A is small compared to the price of product B and that the effect of the material cost on the overall cost is negligible.

As an objective function, we define a function my_obj_func that returns the negative value of the total production of B obtained from the simulation. This is because FMQA optimizes to minimize the value of the objective function.

2. FMQA program implementation¶

This section describes the program implementation of FMQA, which is identical to the implementation in "Black-Box Optimization with Quantum Annealing and Ising Machines", so please refer to that for details.

2.1．Random seed initialization¶

We define a function seed_everything() to initialize random seed values to ensure that the machine learning results do not change with each run.

2.2．Configuration of Amplify client¶

Here, we create an Amplify client and set the necessary parameters. In the following, we set the timeout for a single search by the Ising machine to 1 second.

2.3．Implementing FM with PyTorch¶

In this example code, FM is implemented with PyTorch. In the TorchFM class, we define the acquisition function $g(x)$ as a machine learning model. Each term in $g(x)$ corresponds directly to out_lin, out_1, out_2, and out_inter in the TorchFM class, as in the following equation.

$$\begin{aligned} g(\boldsymbol{x} | \boldsymbol{w}, \boldsymbol{v}) &= \underset{\color{red}{\mathtt{out\_lin}}}{\underline{ w_0 + \sum_{i=1}^n w_i x_i} } + \underset{\color{red}{\mathtt{out\_inter}}}{\underline{\frac{1}{2}\left(\underset{\color{red}{\mathtt{out\_1}}}{\underline{ \sum_{f=1}^k\left(\sum_{i=1}^n v_{i f} x_i\right)^2 }} - \underset{\color{red}{\mathtt{out\_2}}}{\underline{ \sum_{f=1}^k\sum_{i=1}^n v_{i f}^2 x_i^2 }} \right) }} \end{aligned}$$

Next, a function train() is defined to train the FM based on the training data sets. As in general machine learning methods, this function divides the data sets into training data and validation data, then optimizes the FM parameters using the training data, and validates the model during training using the validation data. The train() function returns the model with the highest prediction accuracy for the validation data.

2.4．Construction of initial training data¶

The gen_training_data function evaluates the objective function $f(\boldsymbol{x})$ against the input value $\boldsymbol{x}$ to produce $N_0$​​ input-output pairs (initial training data). The input value $\boldsymbol{x}$ can be determined in a variety of ways, such as by using a random number or a value suitable for machine learning based on prior knowledge. You can also build up the training data from the results of previous experiments or simulations.

2.5．Execution class for FMQA cycle¶

FMQA.cycle() executes an FMQA cycle that is performed for $N−N_0$​​ times using the pre-prepared initial training data. FMQA.step() is a function that executes only one FMQA cycle, and is called $N−N_0$​​ times by FMQA.cycle().

3．Search for optimal operating conditions¶

3.1. FMQA execution example¶

Now, using the reactor simulator introduced in section 1.4 as the objective function, we will perform optimization to maximize the total amount of B produced in a given time (minimize the negative value of the total amount) by the FMQA implemented in section 2.

The objective function is evaluated $N=30$ times, of which $N_0$=20 for initial data generation. Thus, in the example below, the FMQA cycle (machine learning, seeking the optimal solution in a QUBO manner, and the objective function evaluation) is performed $N-N_0=10$ times. With this setup, it takes approximately 1-5 minutes to complete the optimization.

3.2. Transition of objective function values during the optimization process¶

The following line displays the evolution of the objective function values during the optimization process (see the output example in "3.3. Example output from this FMQA sample program").

The initial $N_0$​ objective function values (blue line) are obtained from randomly generated input values during initial training data generation. The following red line shows the objective function values during the $N−N_0$​ FMQA optimization cycles.

The blue and red lines show how the smallest objective function value is successively updated from the currently optimal input value (red line) obtained by the FMQA optimization cycle.

In general, due to the principle of the heuristics algorithm employed in FixstarsClient, the solutions obtained are not perfectly reproducible, but when solved for the parameters in the above sample code, the resulting solution (initial concentration distribution of A) exceeds 0.8 for the total production of B. Compared to a random search (blue), this shows a substantial improvement in production.

3.3. Example output from this FMQA sample program¶

In general, due to the principle of the heuristics algorithm employed in FixstarsClient, the solutions obtained are not completely reproducible, but the following is a typical standard output and image output obtained when this sample code is executed. The values obtained may vary.

• Without changing the conditions in "FMQA execution example", the following standard output is sequentially output as the FMQA cycle progresses. The following figure is also output as a simulation result based on the optimized initial distribution of A.

  Generating 0-th training data set.
  Generating 10-th training data set.
  Starting FMQA cycles...
  FMQA Cycle #0 variable updated, pred_y=-0.0
  FMQA Cycle #1 variable updated, pred_y=-0.7341318540009673
  FMQA Cycle #2 variable updated, pred_y=-0.7836727189544249
  FMQA Cycle #3 variable updated, pred_y=-0.7862647410081264
  FMQA Cycle #4 
  FMQA Cycle #5 
  FMQA Cycle #6 
  FMQA Cycle #7 variable updated, pred_y=-0.8310535978823115
  FMQA Cycle #8 
  FMQA Cycle #9 
  pred x: [1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 0. 1. 1. 0. 1. 1. 1. 1. 1. 1.
  1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1.
  2. 1. 1. 1. 1. 1. 1. 0. 1. 1. 0. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.
  3. 1. 1. 1. 1. 1. 1. 1. 0. 1. 1. 1. 0. 1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1.
  4. 1. 1. 1.]
  minus_total=-8.31e-01, my_reactor.iter=997, my_reactor.cpu_time=3.7e-01s
  pred value: -0.8310535978823115
  

  

• The output image from fmqa_reactor.plot_history() described in "3.2. Transition of objective function values during optimization process" is as follows.