BlackBox Optimization of Operating Condition in a Chemical Reactor¶
As a practical example of blackbox 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 blackbox optimization and FMQA, see "BlackBox Optimization Exploration of Model Superconducting Materials". For another advanced application case using FMQA, see "BlackBox Optimization of Airfoil Geometry by Fluid Flow Simulation".
This notebook contains the following sections.
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, blackbox 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 threedimensional (3D) space, we will focus on onedimensional (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 zdirections
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 (1C_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}$.
import matplotlib.pyplot as plt
import numpy as np
import time
# A class to solve the transport equations (partial differential equations) for the concentration of a given reactive substance (A>B) by finite difference methods.
class Reactor:
def __init__(self, nfolds=5, alpha=1.0, dt_prod=1e3, rr=1e4):
self.nfolds = nfolds # Parameters that determine the spatial resolution of the simulation domain
self.alpha = alpha # Molecular diffusion coefficient
self.dt_prod = dt_prod # Predetermined production time
self.rr = rr # Coefficients related to the reaction rate
# Function to compute the secondorder derivative of the distribution f by secondorder central differencing (assuming periodic boundary conditions)
def __dfdx2(self, f, dfdx2):
dfdx2[0] = (f[1]  2 * f[0] + f[self.nx  1]) / self.dx / self.dx
dfdx2[self.nx  1] = (
(f[0]  2 * f[self.nx  1] + f[self.nx  2]) / self.dx / self.dx
)
dfdx2[1 : self.nx  1] = (
np.array([f[i + 1]  2 * f[i] + f[i  1] for i in range(1, self.nx  1)])
/ self.dx
/ self.dx
)
return dfdx2
# Determine the initial conditions
def __init_field(self, x):
self.nx = self.nfolds * len(
x
) # Number of mesh points used in the spatial discretization
self.dx = 1.0 / (self.nx  1) # Mesh point spacing
self.concn_A = np.zeros(self.nx) # Concentration distribution of A
self.concn_B = np.zeros(self.nx) # Concentration distribution of B
self.x_cord = np.array(
[i / self.nx  0.5 for i in range(self.nx)]
) # Coordinate of discrete points
self.concn_A = np.array(
[float(x[i]) for i in range(len(x)) for j in range(self.nfolds)]
) # Generate the initial field for A
# Function to evolve the transport equation in time by dt_prod physical time according to the initial distribution init_A of A and return the total amount of B produced
def integrate(self, init_A, fig=False):
self.__init_field(init_A)
start = time.perf_counter()
omega = np.zeros(self.nx)
dfdx2 = np.zeros(self.nx)
dt = 0.25 * self.dx * self.dx / self.alpha # Time step width in Eulerian method
lts = int(self.dt_prod / dt)
if fig: # Plot of reaction progress
fig = plt.figure(figsize=(6, 4))
plt.tick_params(labelsize=16)
plt.xlabel("x", fontsize=16)
plt.ylabel("Concentration", fontsize=18)
plt.plot(
self.x_cord,
self.concn_A,
linestyle="",
linewidth=1,
color=[0.6, 0.6, 0.6],
label="$C_{A,0}$",
)
self.iter = 0
while self.iter * dt < self.dt_prod:
if fig and any(
[
self.iter == i
for i in [0, int(0.1 * lts), int(0.2 * lts), int(0.4 * lts)]
]
): # Plot of reaction progress
plt.plot(
self.x_cord, self.concn_B, linestyle="", linewidth=2, color="r"
)
omega = (
self.rr * np.exp(self.concn_B) * self.concn_A * (1.0  self.concn_A)
) # Reaction rate
self.concn_A = (
self.concn_A
+ (self.alpha * self.__dfdx2(self.concn_A, dfdx2)  omega) * dt
) # Time advancement for the concentration of A
self.concn_B = (
self.concn_B
+ (self.alpha * self.__dfdx2(self.concn_B, dfdx2) + omega) * dt
) # Time advancement for the concentration of B
self.iter += 1
if fig: # Plot of reaction progress
plt.plot(
self.x_cord,
self.concn_B,
linestyle="",
linewidth=4,
color="r",
label="$C_B$",
)
plt.legend(fontsize=16)
self.cpu_time = time.perf_counter()  start # Measure the computation time
return (
np.sum(self.concn_B) * self.dx
) # Simplified spatial integration of the concentration of B
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?
# Generate 1D binary array representing C_{A,0} by random numbers
# Discretize 1D space into 100 finite regions to represent the initial distribution of C_A
# Random seed is not fixed, so the contents of c_a0 change with each run
c_a0 = np.random.randint(0, 2, 100)
# Reaction simulation is performed by the `integrate`` function of the Reactor class (option to output result image)
amount_B = Reactor().integrate(c_a0, fig=True)
# Output the total amount of B produced
print(f"Total amount of B: {amount_B:.2f}")
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.
# Objective function (returns the negative value of the total amount of product B produced)
def my_obj_func(init_A, fig=False):
my_reactor = Reactor()
minus_total = my_reactor.integrate(init_A, fig)
if fig:
# (Optional) Displays objective function value, number of integrations in integrate(), and CPU time required for the simulation
print(f"{minus_total=:.2e}, {my_reactor.iter=}, {my_reactor.cpu_time=:.1e}s")
return minus_total
# Example 1: objective function value when a certain binary vector c_a0 (defined in the cell above) is fed
amount_B = my_obj_func(c_a0)
print(f"{amount_B=:.2f}")
# Example 2: Objective function value when a certain binary vector c_a0 (defined in the cell above) is fed (log and image as well)
amount_B = my_obj_func(c_a0, fig=True)
print(f"{amount_B=:.2f}")
2. FMQA program implementation¶
This section describes the program implementation of FMQA, which is identical to the implementation in "BlackBox 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.
import os
import torch
import numpy as np
def seed_everything(seed=0):
os.environ["PYTHONHASHSEED"] = str(seed)
np.random.seed(seed)
torch.manual_seed(seed)
torch.cuda.manual_seed(seed)
torch.backends.cudnn.deterministic = True
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.
from amplify import FixstarsClient
from datetime import timedelta
client = FixstarsClient()
client.parameters.timeout = timedelta(milliseconds=1000) # timeout is 1000 ms
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
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} $$
import torch.nn as nn
class TorchFM(nn.Module):
def __init__(self, d: int, k: int):
super().__init__()
self.V = nn.Parameter(torch.randn(d, k), requires_grad=True)
self.lin = nn.Linear(
d, 1
) # The first and second terms on the righthand side are fully connected network
def forward(self, x):
out_1 = torch.matmul(x, self.V).pow(2).sum(1, keepdim=True)
out_2 = torch.matmul(x.pow(2), self.V.pow(2)).sum(1, keepdim=True)
out_inter = 0.5 * (out_1  out_2)
out_lin = self.lin(x)
out = out_inter + out_lin
return out
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.
from torch.utils.data import TensorDataset, DataLoader
from sklearn.model_selection import train_test_split
from typing import Type
import copy
def train(
X: np.ndarray,
y: np.ndarray,
model_class: Type[nn.Module],
model_params: dict[str, int  float],
batch_size=1024,
epochs=3000,
criterion=nn.MSELoss(),
optimizer_class=torch.optim.AdamW,
opt_params={"lr": 1},
lr_sche_class=None,
lr_sche_params=None,
):
X_tensor, y_tensor = (
torch.from_numpy(X).float(),
torch.from_numpy(y).float(),
)
indices = np.array(range(X.shape[0]))
indices_train, indices_valid = train_test_split(
indices, test_size=0.2, random_state=42
)
train_set = TensorDataset(X_tensor[indices_train], y_tensor[indices_train])
valid_set = TensorDataset(X_tensor[indices_valid], y_tensor[indices_valid])
loaders = {
"train": DataLoader(train_set, batch_size=batch_size, shuffle=True),
"valid": DataLoader(valid_set, batch_size=batch_size, shuffle=False),
}
model = model_class(**model_params)
best_model_wts = copy.deepcopy(model.state_dict())
optimizer = optimizer_class(model.parameters(), **opt_params)
scheduler = None
if lr_sche_class is not None:
scheduler = lr_sche_class(optimizer, **lr_sche_params)
best_score = 1e18
for _ in range(epochs):
losses = {"train": 0.0, "valid": 0.0}
for phase in ["train", "valid"]:
if phase == "train":
model.train()
else:
model.eval()
for batch_x, batch_y in loaders[phase]:
optimizer.zero_grad()
out = model(batch_x).T[0]
loss = criterion(out, batch_y)
losses[phase] += loss.item() * batch_x.size(0)
with torch.set_grad_enabled(phase == "train"):
if phase == "train":
loss.backward()
optimizer.step()
losses[phase] /= len(loaders[phase].dataset)
with torch.no_grad():
model.eval()
if best_score > losses["valid"]:
best_model_wts = copy.deepcopy(model.state_dict())
best_score = losses["valid"]
if scheduler is not None:
scheduler.step()
with torch.no_grad():
model.load_state_dict(best_model_wts)
model.eval()
return model
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$ inputoutput 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.
def gen_training_data(D: int, N0: int, true_func):
assert N0 < 2**D
# N0 input values are obtained using random numbers
X = np.random.randint(0, 2, size=(N0, D))
# Remove duplicate input values and add new input values using random numbers
X = np.unique(X, axis=0)
while X.shape[0] != N0:
X = np.vstack((X, np.random.randint(0, 2, size=(N0  X.shape[0], D))))
X = np.unique(X, axis=0)
y = np.zeros(N0)
# Obtain output values corresponding to N0 input values by evaluating the objective function, true_func
for i in range(N0):
if i % 10 == 0:
print(f"Generating {i}th training data set.")
y[i] = true_func(X[i])
return X, y
2.5．Execution class for FMQA cycle¶
FMQA.cycle()
executes an FMQA cycle that is performed for $N−N_0$ times using the
preprepared initial training data. FMQA.step()
is a function that executes only one FMQA
cycle, and is called $N−N_0$ times by FMQA.cycle()
.
from amplify import VariableGenerator, solve
import matplotlib.pyplot as plt
import sys
class FMQA:
def __init__(self, D: int, N: int, N0: int, k: int, true_func, client) > None:
assert N0 < N
self.D = D
self.N = N
self.N0 = N0
self.k = k
self.true_func = true_func
self.client = client
self.y = None
# A member function that repeatedly performs (NN0)x FMQA based on the training data with adding new training data
def cycle(self, X, y, log=False) > np.ndarray:
print(f"Starting FMQA cycles...")
pred_x = X[0]
pred_y = 1e18
for i in range(self.N  self.N0):
print(f"FMQA Cycle #{i} ", end="")
try:
x_hat = self.step(X, y)
except RuntimeError:
sys.exit(f"Unknown error, i = {i}")
# If an input value identical to the found x_hat already exists in the current training data set, a neighboring value is used as a new x_hat.
is_identical = True
while is_identical:
is_identical = False
for j in range(i + self.N0):
if np.all(x_hat == X[j, :]):
change_id = np.random.randint(0, self.D, 1)
x_hat[change_id.item()] = 1  x_hat[change_id.item()]
if log:
print(f"{i=}, Identical x is found, {x_hat=}")
is_identical = True
break
# Evaluate objective function f() with x_hat
y_hat = self.true_func(x_hat)
# Add an inputoutput pair [x_hat, y_hat] to the training data set
X = np.vstack((X, x_hat))
y = np.append(y, y_hat)
# Copy the inputoutput pair to [pred_x, pred_y] when the evaluated value of the objective function updates the minimum value
if pred_y > y_hat:
pred_y = y_hat
pred_x = x_hat
print(f"variable updated, {pred_y=}")
else:
print("")
# Exit the "for" statement if all inputs have been fully explored
if len(y) >= 2**self.D:
print(f"Fully searched at {i=}. Terminating FMQA cycles.")
break
self.y = y
return pred_x
# Member function to perform one FMQA cycle
def step(self, X, y) > np.ndarray:
# Train FM
model = train(
X,
y,
model_class=TorchFM,
model_params={"d": self.D, "k": self.k},
batch_size=8,
epochs=2000,
criterion=nn.MSELoss(),
optimizer_class=torch.optim.AdamW,
opt_params={"lr": 1},
)
# Extract FM parameters from the trained FM model
v, w, w0 = list(model.parameters())
v = v.detach().numpy()
w = w.detach().numpy()[0]
w0 = w0.detach().numpy()[0]
# Solve a QUBO problem using a quantum annealing or Ising machine
gen = VariableGenerator() # Declare a variable generator
q = gen.array("Binary", self.D) # Generate binary decision variables
model = self.__FM_as_QUBO(q, w0, w, v) # Define FM as a QUBO equation
result = solve(model, self.client) # Pass the objective function to Amplify
if len(result.solutions) == 0:
raise RuntimeError("No solution was found.")
q_values = q.evaluate(result.best.values)
return q_values
# A function that defines FM as a QUBO equation from FM parameters. As with the previously defined TorchFM class, the formula is written as per the acquisition function form of g(x).
def __FM_as_QUBO(self, x, w0, w, v):
lin = w0 + (x.T @ w)
out_1 = np.array([(x * v[:, i]).sum() ** 2 for i in range(self.k)]).sum()
# Note that x[j] = x[j]^2 since x[j] is a binary variable in the following equation.
out_2 = np.array([(x * v[:, i] * v[:, i]).sum() for i in range(self.k)]).sum()
return lin + (out_1  out_2) / 2
# A function to plot the history of ith objective function evaluations performed within the initial training data construction (blue) and during FMQA cycles (red).
def plot_history(self):
assert self.y is not None
fig = plt.figure(figsize=(6, 4))
plt.plot(
[i for i in range(self.N0)],
self.y[: self.N0],
marker="o",
linestyle="",
color="b",
) # Objective function evaluation values at the time of initial training data generation (random process)
plt.plot(
[i for i in range(self.N0, self.N)],
self.y[self.N0 :],
marker="o",
linestyle="",
color="r",
) # Objective function evaluation values during the FMQA cycles (FMQA cycle process)
plt.xlabel("ith evaluation of f(x)", fontsize=18)
plt.ylabel("f(x)", fontsize=18)
plt.tick_params(labelsize=18)
return fig
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 $NN_0=10$ times. With this setup, it takes approximately 15 minutes to complete the optimization.
# Initialize random seed values
seed_everything()
D = 100 # Size of input values (problem size)
N = 30 # Number of times the function can be evaluated
N0 = 20 # Number of samples of initial training data
k = 20 # Dimension of the vector in FM (hyperparameters)
# Generate initial training data
X, y = gen_training_data(D, N0, my_obj_func)
# Instantiate FMQA class
fmqa_reactor = FMQA(D, N, N0, k, my_obj_func, client)
# Run FMQA cycle
pred_x = fmqa_reactor.cycle(X, y, log=True)
# Output optimization results
print("pred x:", pred_x)
print("pred value:", my_obj_func(pred_x, fig=True))
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.
fig = fmqa_reactor.plot_history()
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 0th training data set. Generating 10th 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.31e01, my_reactor.iter=997, my_reactor.cpu_time=3.7e01s 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.