BlackBox Optimization of the Airfoil Geometry with Fluid Flow Simulation¶
Various types of airfoils are used in aircraft, transportation equipment, and various energy devices, and the crosssectional shape (airfoil shape) has a significant impact on the performance and efficiency of these industrial products.
Generally, airfoil optimization can be performed by various methods (trial and error) to maximize lift (or downforce) and minimize drag. However, evaluating the airfoil profile (computing/measuring lift and drag forces) requires experiments or fluid simulations, and the evaluation cost (in terms of time and money) per evaluation is considered relatively high. Therefore, it is difficult to optimize an airfoil profile based on multiple parameters in a full search manner.
In this tutorial, blackbox optimization based on a factorization machine and simulated annealing will be used to optimize an airfoil profile with a limited number of evaluations. In the optimization process, the ratio of lift $F_L$ and drag $F_D$ acting on an airfoil situated in a flow field (liftdrag ratio) is maximized. The liftdrag ratio $r_{LD}$ is defined as:
$$ r_{LD} = F_L/F_D $$
There are several ways to construct the airfoil shape. Here, the airfoil shape is constructed by applying a mapping called the Jukowski transformation to a circle in the complex plane. The negative value of the liftdrag ratio ($r_{LD}$) is used as the objective function, and fluid flow simulation is used to obtain (evaluate) the value of the objective function.
In addition, optimization is performed for nonbinary decision variables, which in this case take integers as input variables, which are then used to express discrete real numbers. To describe the objective function and the acquisition function in QUBO formulas, it is necessary to convert the nonbinary variables into binary decision variables. For this purpose, this example program uses onehot encoding and also introduces the necessary encoder/decoder implementation.
For basic knowledge of blackbox optimization and FMQA, see "BlackBox Optimization with Quantum Annealing/Ising Machines". For the other practical application case of blackbox optimization, see "BlackBox Optimization of Operating Condition in a Chemical Reactor".
This notebook consists of the following sections.
 1. Airfoil definition and fluid flow simulation
 2. FMQA program implementation (nonbinary decision variables)
 3. Search for the optimal airfoil profile
$^*$ The execution time is limited to about 20 minutes in this online demo and tutorial environment. If necessary, please run the program with a smaller number of optimization cycles $N$ (see below) or copy this example program to your own environment. In this case, you should also download the following libraries related to airfoil generation and fluid flow simulation as well, save them in the following directory structure, and execute the example code.
├ fmqa_3_aerofoil.ipynb
(this sample program)
└ utils/
├ __init__.py
(blank)
├ joukowski.py
└ lbm.py
1. Airfoil definition and fluid flow simulation¶
The following airfoil shape generator and fluid analysis classes are provided in this example program.
 Airfoil generator class:
joukowski.WingGenerator
 Fluid flow simulator class:
lbm.Solver
In this section, we will use the prepared classes to simulate fluid flow around an airfoil and explain how to obtain the liftdrag ratio $r_{L/D}$, which is directly related to the objective function in this blackbox optimization example.
1.1. Airfoil generator class¶
This example program uses the Joukowski transform as the method for constructing the airfoil. The Joukowski transform is,
$$ z(x,y) = \zeta(\xi, \eta) + \frac{c^2}{\zeta(\xi, \eta)} $$
and represents a map between two complex planes in $\zeta\rightarrow z$. The airfoil profile in this example program is constructed by the Jukowski transform, based on a circle in $\zeta$ whose center is $(\xi_0, \eta_0)$ and passes through the point $(c,0)$, which is then rotated by $\alpha$ in the $z$ plane. The $\xi_0$, $\eta_0$, and $\alpha$ are parameters, corresponding to the wing thickness, warping, and angle of attack, respectively.
Let us generate an airfoil using joukowski.WingGenerator
, which is an airfoil generator
class. The parameters $(\xi_0, \eta_0, \alpha)$ for the mapping are real numbers, but for the sake of
discrete optimization, we will choose each parameter from a finite number of candidates. Such an
operation
is called integer indexing.
from utils import joukowski
# Instantiate the airfoil generator class, joukowski.WingGenerator
wing_generator = joukowski.WingGenerator()
# Display integer index table of airfoil generator parameters
wing_generator.show_index_table()
With wing_generator.show_index_table()
, you can check the correspondence between the
integer
indices used for airfoil generation and the parameter values $(\xi_0, \eta_0, \alpha)$. For example,
$\xi_0=5.50$, $\eta_0=3.75$, $\alpha=20.00$ deg can be specified as $\xi_0[10]=5.50$,
$\eta_0[15]=3.75$,
$\alpha[20]=20.00$ deg using integer indices. Using these parameters, the airfoil profile will be
generated and displayed as follows.
Let us construct an airfoil using different integer index values and see the mapping parameters and the airfoil profile.
# Perform airfoil generation (mapping transformation). xi0, eta0, and alpha are specified by integer indices. See the output of wing_generator.show_index_table() for the index table.
wing = wing_generator.generate_wing(10, 15, 20)
# Show the constructed airfoil profile
wing.draw()
# Show the mapping parameters for transformation
wing_generator.show_parameters("My wing")
1.2. Fluid flow simulator class¶
In this tutorial, we will use a fluid flow simulator to evaluate airfoil profiles. Various fluid
simulation methods have been developed, but in this tutorial, we will use a twodimensional fluid flow
simulation method based on the Lattice Boltzmann method due to its computational cost and ease of
setting
boundary conditions. A Lattice Boltzmann solver is available from the fluid flow simulator class
lbm.Solver
.
As an example, the flow around the airfoil constructed above is simulated as follows. Note that under the given conditions, one simulation will take about 20 seconds of computation time. After the simulation, the flow velocity vector (black arrow) and the fluid force acting on the airfoil (combined lift and drag force; red arrow) are shown in the visualization. The color indicates the vorticity.
This time, we aim to maximize the ratio of the lift force $F_L$ and drag force $F_D$ (liftdrag ratio) $r_{LD} = F_L/F_D$, so we consider the negative value of the liftdrag ratio $r_{LD}$ as the objective function and perform the optimization to minimize it. The evaluated value of $r_{LD}$ will be shown below at the end of the simulation.
from utils import lbm
# Simulation grid size (nx: horizontal, ny: vertical)
nx, ny = 400, 160
# Instantiate the lbm.Solver class.
solver = lbm.Solver(nx, ny)
# Add the airfoil, wing, constructed above as a boundary condition (profile model)
solver.add_model(wing)
# Run 3000timestep simulation
solver.integrate(steps=3000)
# Calculate the negative value of the liftdrag ratio, the objective function, from the simulation results (for this FMQA, minimize this value)
cost = solver.fy_average / solver.fx_average
# Print the evaluated value of the objective function
print(f"{cost=:.2f}")
2. FMQA program implementation (nonbinary decision variables)¶
In this section, we will implement the blackbox optimization (FMQA). Much of the implementation of the FMQA part can be found in "BlackBox Optimization with Quantum Annealing/Ising Machines", so please refer to that document for details.
Please note that the main difference is that the present blackbox optimization considers integer vectors as input values to the objective function instead of binary variables. For the handling of integer input values in FMQA, see "Onehot encoder/decoder".
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¶
This example code implements FM with PyTorch. In the TorchFM
class, we define the
acquisition function $g(x)$ as a machine learning model. Each term in $g(x)$ directly corresponds 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 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.
To handle nonbinary integer input values, we use the input_generator
function, which
generates an integer input vector with random numbers and converts it to a binary vector based on
onehot
encoding. The function input_generator
is defined in the onehot encoder/decoder class
EncDecOneHot
, which will be introduced later.
def gen_training_data_flow_onehot(D: int, N0: int, true_func, input_generator):
X = np.zeros((N0, D))
y = np.zeros(N0)
print("")
print("###################################################")
print(" Making initial training data")
print("###################################################")
for i in range(N0):
print(f"Generating training data set #{i}.")
# Randomly determine input vector consisting of (0 or 1)
x = input_generator()
# If the same input vector has already exist, regenerate the input values
is_identical = True
while is_identical:
is_identical = False
for j in range(i + 1):
if np.all(x == X[j, :]):
x = input_generator()
is_identical = True
break
# Evaluate input value x with objective function f(x) and copy inputoutput pair (x,f(x)) to teacher data
X[i, :] = x
y[i] = true_func(x)
print(f"")
return X, y
2.5. Execution class for FMQA cycle¶
FMQA.cycle()
executes an FMQA cycle that is run 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()
.
For FMQA with nonbinary variables, each variable must be encoded as a binary variable. The
FMQA_NB
class assumes onehot encoding. The encoder and decoder of nonbinary variables
are
specified at the instance of this class as encoder_decoder.encoder
and
encoder_decoder.decoder
, respectively. Also, nbins_array
is a list of the
total
(variable) number of bins used to encode each nonbinary variable.
from amplify import VariableGenerator, solve, one_hot, sum
import matplotlib.pyplot as plt
import sys
# FMQA for nonbinary inputs
class FMQA_NB:
def __init__(self, D, N, N0, k, true_func, client, encoder_decoder) > 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.nbins_array = encoder_decoder.nbins_array
self.encoder = encoder_decoder.encoder
self.decoder = encoder_decoder.decoder
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:
# Weight for the onehot constraints
constraint_weight = max([1, int(np.abs(y).max() + 0.5)])
print("")
print("###################################################")
print(f" Starting FMQA cycles... {constraint_weight=}")
print("###################################################")
pred_x = X[0]
pred_y = 1e18
for i in range(self.N  self.N0):
print(f"FMQA Cycle #{i} ")
try:
x_hat = self.step(X, y, constraint_weight)
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 the 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, :]):
# Decode binary vector to integer vector and copy to inputss
inputs = self.decoder(x_hat)
# Get the surrounding values (inputs is an integer vector)
inputs += np.random.randint(1, 2, len(inputs))
for i_inp in range(len(inputs)):
if inputs[i_inp] < 0:
inputs[i_inp] = 0
elif inputs[i_inp] > self.nbins_array[i_inp]  1:
inputs[i_inp] = self.nbins_array[i_inp]  1
# Encode from integer vector to binary vector and copy to x_hat
x_hat = self.encoder(inputs)
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 the 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("")
print(f"")
self.y = y
return pred_x
# Member function to perform one FMQA cycle
def step(self, X, y, constraint_weight) > 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,
# Set scheduler to reduce learning rate with number of epochs
opt_params={"lr": 0.5},
lr_sche_class=torch.optim.lr_scheduler.StepLR,
lr_sche_params={"step_size": 50, "gamma": 0.9},
)
# 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
cost = self.__FM_as_QUBO(q, w0, w, v) # Define FM as a QUBO equation
# Construct onehot constraints for nonbinary variables
constraint_list = []
ista = 0
iend = 0
for i in range(len(self.nbins_array)):
iend += self.nbins_array[i]
constraint_list.append(one_hot(q[ista:iend]))
ista = iend
constraints = sum(constraint_list)
# Add up the objective function and constraints, and pass it to the Fixstars Amplify solver
model = cost + constraint_weight * constraints
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
"""The sum_poly used in __FM_as_QUBO above is inefficient in terms of computation speed and
memory. In the case of FM, where the interaction terms of the decision variables are generally
nonzero, the following implementation using BinaryMatrix is more efficient. Here, the quadratic
terms in BinaryMatrix correspond to the nondiagonal terms represented by the upper triangular
matrix, so x(1/2) for the quadratic terms in the FM formula is unnecessary. Also, although x is
taken as an argument just to match the function signature with __FM_as_QUBO above (implementation
using sum_poly), it is not needed in this implementation using BinaryMatrix.
def __FM_as_QUBO(self, x, w0, w, v):
out_1_matrix = v @ v.T
out_2_matrix = np.diag((v * v).sum(axis=1))
matrix = BinaryMatrix(out_1_matrix  out_2_matrix + np.diag(w))
# Do not forget to put the constant term w0 in the second argument of BinaryQuadraticModel.
model = BinaryQuadraticModel(matrix, w0)
return model
"""
# 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
2.6. Onehot encoder/decoder¶
When using the generate_wing
method of the airfoil generator class, an integer index
(e.g.
generate_wing(10, 15, 20)
) is required as input. On the other hand, the input to QUBO
must be
a binary variable, so encoding from integer indices to binary variables is required. In this tutorial,
we
will use a onehot encoding. For example, in onehot encoding with a 4bit representation of an
integer
that can take values from 1 to 4,
 1 is encoded to
[1, 0, 0, 0]
 2 is encoded to
[0, 1, 0, 0]
 3 is encoded to
[0, 0, 1, 0]
 4 is encoded to
[0, 0, 0, 1]
The following EncDecOneHot
class defines functions for converting integer index to
binary
variable and vice versa (encoder()
, decoder()
) and functions for generating
an
integer index input vector with random numbers and converting it to a binary variable matrix
gen_random_input()
.
class EncDecOneHot:
def __init__(self, D: int, nbins_array):
self.D = D
self.nbins_array = nbins_array
def gen_random_input(self):
ista = 0
iend = 0
x = np.zeros(self.D, int)
for i in range(len(self.nbins_array)):
iend += self.nbins_array[i]
idx = np.random.randint(ista, iend, 1)[0]
x[idx] = 1
ista = iend
return x
def encoder(self, inputs):
if len(inputs) != len(self.nbins_array):
raise RuntimeError("inputs should be the same length as nbins_array!")
x = np.zeros(self.D, int)
i_offset = 0
for i in range(len(inputs)):
x[inputs[i] + i_offset] = 1
i_offset += self.nbins_array[i]
return x
def decoder(self, x):
bit_locs = np.where(x == 1)[0]
if len(bit_locs) != len(self.nbins_array):
raise RuntimeError("bit_locs should be the same length as nbins_array!")
i_offset = 0
for i in range(1, len(bit_locs)):
i_offset += self.nbins_array[i  1]
bit_locs[i] = i_offset
return bit_locs
3. Search for the optimal airfoil profile¶
3.1. FMQA execution example¶
Now, with the fluid flow simulator introduced in Section 1 as the objective function, the optimization to minimize the negative value of the liftdrag ratio is performed by FMQA implemented in Section 2.
In the following code, the number of times $N$ the objective function can be evaluated is set to 5, of which $N_0=3$ is the number of evaluations to generate the initial training data. This setting is for a test purpose and you may need to use more appropriate $N$ and $N_0$ for meaningful optimization.
In the example output shown in "3.3. Example output from this FMQA sample program", $N=40$ and $N_0=20$ are used so that meaningful optimization can be achieved. Note that with this setup, it takes approximately 20 minutes to complete all FMQA cycles.
# Initialize random seed value
seed_everything()
# Instantiate airfoil generator class, joukowski.WingGenerator
wing_generator = joukowski.WingGenerator()
# Display integer index table of airfoil generator parameters
wing_generator.show_index_table()
# Onehot encoding bin size for each variable
nbins_array = np.array(
[wing_generator.nbins_xi0, wing_generator.nbins_eta0, wing_generator.nbins_alpha]
)
# Overall size of binary variables (= size of binary decision variables handled by QUBO)
D = wing_generator.nbins_xi0 + wing_generator.nbins_eta0 + wing_generator.nbins_alpha
# Instantiate onehot encoderdecoder class (EncDecOneHot)
enc_dec_one_hot = EncDecOneHot(D, nbins_array)
# Counter for the number of objective function evaluations
n_eval = 0
# Objective function (construct airfoil → run simulation → compute objective function value)
def obj_func(x):
global n_eval
# Decode binary variables to integers (mapping parameter indices)
i_xi0, i_eta0, i_alpha = enc_dec_one_hot.decoder(x)
# Perform airfoil generation (mapping transformation). xi0, eta0, alpha are integer indices.
# See the output of wing_generator.show_index_table() for table of indices
wing = wing_generator.generate_wing(i_xi0, i_eta0, i_alpha)
# Show the mapping parameters for the transformation
wing_generator.show_parameters("My wing")
# Simulation grid size (nx: horizontal, ny: vertical)
nx, ny = 400, 160
# Instantiate the lbm.Solver class.
solver = lbm.Solver(nx, ny)
# Add the airfoil, wing, constructed above as a boundary condition (airfoil model)
solver.add_model(wing)
# Run 5000timestep simulation
solver.integrate(steps=5000)
# Negative value of the liftdrag ratio
cost = solver.fy_average / solver.fx_average
# If the fluid simulation diverges
if np.isnan(cost):
cost = 0
print(f"#{n_eval} evaluation finished. cost:{cost:.2f}")
n_eval += 1
return cost
N = 5 # 40 is used for Sec. 3.3. # Number of times the objective function can be evaluated
N0 = 3 # 20 is used for Sec. 3.3. # Number of samples of initial training data
k = 20 # Dimension of vector in FM (hyperparameters)
# Generate initial training data
X, y = gen_training_data_flow_onehot(D, N0, obj_func, enc_dec_one_hot.gen_random_input)
# Instantiate FMQA
fmqa_solver = FMQA_NB(D, N, N0, k, obj_func, client, enc_dec_one_hot)
# Perform FMQA cycle
pred_x = fmqa_solver.cycle(X, y)
# Output optimization results
print("###################################################")
print(" Optimal input values and corresponding output")
print("###################################################")
print("pred x:", pred_x)
print("pred value:", obj_func(pred_x))
3.2. Transition of objective function values during the optimization process¶
The following line plots the evolution of the objective function values during the optimization process. The initial $N_0$ objective function values (blue line) are obtained from randomly generated input values during the 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 when large enough $N$ is considered.
Here, input values $\hat{x}$ can sometimes be obtained that produce relatively large objective function values, even during the FMQA cycle. There are several reasons for this, but the main ones are
 If the current training data contains the same input value as the newly searched input value, the
neighboring input value $\hat{x}$ is ramdomly searched and the objective function is evaluated (see
FMQA_NB.cycle
),  The accuracy of the acquisition function $g(\boldsymbol{x})$ at a certain FMQA cycle is not high enough,
In any case, we were able to efficiently perform airfoil optimization by blackbox optimization!
fig = fmqa_solver.plot_history()
3.3. Example output from this FMQA sample program¶
In general, due to the principle of the heuristic algorithm used in FixstarsClient
, the
solutions obtained are not completely reproducible, but the following is a typical standard output and
image output obtained when running this example code with the following $N$ and $N_0$. The values
obtained
may vary.
N = 40 # Number of times the objective function can be evaluated N0 = 20 # Number of samples of initial training data

Executing "FMQA execution example" will output followings as the FMQA cycle progresses.

First, the airfoil profiles constructed for the input values based on random numbers generated during the initial training data generation and their fluid flow simulation results are sequentially output ($N_0$ times).
[Index table] 0 1 2 3 4 5 6 7 8 9 10 11 xi0[i] 1.0 1.45 1.9 2.35 2.8 3.25 3.7 4.15 4.6 5.05 5.5 5.95 \ eta0[i] 0.0 0.25 0.5 0.75 1.0 1.25 1.5 1.75 2.0 2.25 2.5 2.75 alpha[i] 0.0 1.0 2.0 3.0 4.0 5.0 6.0 7.0 8.0 9.0 10.0 11.0 12 13 14 15 16 17 18 19 20 21 22 xi0[i] 6.4 6.85 7.3 7.75 8.2 8.65 9.1 9.55    \ eta0[i] 3.0 3.25 3.5 3.75 4.0 4.25 4.5 4.75 5.0 5.25 5.5 alpha[i] 12.0 13.0 14.0 15.0 16.0 17.0 18.0 19.0 20.0 21.0 22.0 23 24 25 26 27 28 29 30 31 32 33 xi0[i]            \ eta0[i] 5.75 6.0 6.25 6.5 6.75 7.0 7.25 7.5 7.75 8.0 8.25 alpha[i] 23.0 24.0 25.0 26.0 27.0 28.0 29.0 30.0 31.0 32.0 33.0 34 35 36 37 38 39 xi0[i]       eta0[i] 8.5 8.75 9.0 9.25 9.5 9.75 alpha[i] 34.0 35.0 36.0 37.0 38.0 39.0 ################################################### Making initial training data ################################################### Generating training data set #0. My wing: xi0[12]=6.40, eta0[0]=0.00, alpha[3]=3.00deg
#0 evaluation finished. cost:0.53  Generating training data set #1. My wing: xi0[3]=2.35, eta0[39]=9.75, alpha[9]=9.00deg
#1 evaluation finished. cost:1.86 
・
・
・ 
After the initial training data construction is completed, an $NN_0$ FMQA cycle begins. During the cycle, the following output is sequentially displayed for each FMQA attempt.
################################################### Starting FMQA cycles... constraint_weight=2 ################################################### FMQA Cycle #0 My wing: xi0[6]=3.70, eta0[23]=5.75, alpha[24]=24.00deg
#20 evaluation finished. cost:2.09 variable updated, pred_y=2.0937129113571267  FMQA Cycle #1 My wing: xi0[3]=2.35, eta0[23]=5.75, alpha[24]=24.00deg
#21 evaluation finished. cost:1.95 
・
・
・ 
After $NN_0$ FMQA cycles, the simulation results of a flow around the airfoil and the objective function values as well as the obtained optimal binary input vectors are shown as follows.
################################################### Optimal input values and corresponding output ################################################### pred x: [0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0] My wing: xi0[3]=2.35, eta0[39]=9.75, alpha[15]=15.00deg
#40 evaluation finished. cost:2.74 pred value: 2.735594392729547


The output image of
fmqa_reactor.plot_history()
described in "3.2. Transition of objective function values during the optimization process shows how the minimum negative values of the liftdrag ratio are updated one after the other. The liftdrag ratio of the optimal airfoil obtained by the blackbox optimization is 2.74, which is 17% higher than the maximum liftdrag ratio (2.35) of the airfoil obtained by the random search (initial training data construction) process.
3.4. Commentary¶
If this example program is executed as is, the resulting optimal liftdrag ratio $r_{LD}$ is approximately 2.7. The liftdrag ratio of a typical aircraft is about 15 in cruise, which is much higher than that of the optimum airfoil obtained by this blackbox optimization. Two reasons for this are as follows.

Jukowski transformation
For the sake of simplicity, the Joukowski transformation, which can define the airfoil shape with three parameters, was used in this study, but the airfoil shape obtained by this transformation is not necessarily the one that achieves a high liftdrag ratio.

Aspect ratio of the airfoil
To increase the liftdrag ratio, the aspect ratio should be increased (thinner and longer airfoils), but simulating thinner and longer objects requires higher resolution, which in turn requires higher computational cost. In this case, the default range of possible mapping parameters has been adjusted accordingly to allow for relatively thick airfoils (small aspect ratio) in relation to the length of the airfoil.
In any case, within the framework of the aspect ratio and Joukowski transformation considered in this example program, we could perform blackbox optimization to find the optimal airfoil profile.