Black-Box Optimization of Model Superconducting Materials¶
To illustrate the effective use of black-box optimization, this sample code describes the optimization for superconducting materials composed of pseudo-materials as an example problem.
Although the present sample code performs material searches based on nonlinear algebraic models, it is possible to perform black-box optimization in the same steps using high-precision simulations or experimental measurement results instead of model algebraic expressions. Even in such cases, you can use this example almost as is.
For a basic introduction to black-box optimization and FMQA, see "Black-Box Optimization with Quantum Annealing and Ising Machines".
Also, more applied, complex optimization problems using black-box optimization are explained in:
- Black-Box Optimization of Operating Condition in a Chemical Reactor
- Black-Box Optimization of Airfoil Geometry by Fluid Flow Simulation
This notebook contains the following sections:
1. Problem setting¶
1.1. Search scenario for superconducting materials¶
Superconductivity technology is expected to be utilized in the fields of transportation, such as maglev trains, metrology, and energy. Various superconducting materials are currently being developed to realize superconductivity.
The temperature at which superconductivity is achieved (critical temperature) is generally around the absolute temperature of 0 K (Kelvin) for currently confirmed superconducting materials. Because of this, superconductivity requires costly cooling to be exploited, and its application in the real world is currently limited. Therefore, the search for high-temperature superconductors is a pressing issue.
Typically, the search for materials that realize superconductivity involves a trial-and-error process of selecting and synthesizing several materials, repeatedly evaluating the critical temperature of the synthesized materials by measurement, and identifying the material to be synthesized that achieves a higher critical temperature. This process of synthesis and critical temperature evaluation is considered time-consuming. For this search, a black box optimization method is used to find a combination of materials close to the optimal solution with a relatively small number of evaluations.
In this example, the search for superconducting materials consisting of pseudo materials is treated as an example to illustrate the material search by a black-box optimization method (FMQA), and a critical temperature model is used to evaluate the critical temperature. Note that the critical temperature models presented below and the combinations of materials obtained are not necessarily physically accurate, and thus the example serves for illustration purposes only.
1.2. Random number 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
1.3. Definition of critical temperature model¶
This example code selects a combination of several materials from $D$ types of materials and performs an optimization to maximize the critical temperature of the superconducting material produced by their synthesis.
In general, the critical temperature can be evaluated by experimental measurement, which requires a relatively large cost (time and money) each time the evaluation is performed.
In this example code, instead of measuring the critical temperature, the following critical
temperature
model supercon_temperature() is used for evaluation. However, this function is only a
substitute for experimental measurement, and its contents and parameters are treated as unknown, and
the
number of calls to supercon_temperature() is also treated as limited.
import numpy as np
import math
# A function to randomly generate various coefficient tables for critical temperature calculations.
def set_properties(size):
mu, sigma, ratio = 0.0, 1.0, 0.2
table1 = np.random.rand(size) * 1e5 * (0.1 * math.log(size) - 0.23)
table2 = np.random.lognormal(mu, sigma, size) * ratio
table3 = np.random.lognormal(mu, sigma, size) * ratio
return table1, table2, table3
# A model function which calculates the critical temperature of the superconducting material to be synthesized from a given combination of materials x (a one-dimensional array of numpy) and the physical properties of each material.
def supercon_temperature(x, debye_table, state_table, interaction_table):
debye_temperature = np.sum(x * debye_table) / np.sum(x)
state_density = np.sum(x * state_table) / np.sum(x)
interaction = np.sum(x * interaction_table) / np.sum(x)
crit_temp = debye_temperature * math.exp(-1.0 / state_density / interaction)
return crit_temp
In the following, the model supercon_temperature(x) for the critical temperature defined
above is used to evaluate the critical temperature of superconducting materials synthesized from a
random
selection of materials. Here, D is the number of materials to be selected, and the input
x is a vector of size D consisting of 0 or 1.
For example, in the case of selecting the first and last of five materials to be combined, the input
vector would be x = [1, 0, 0, 0, 0, 1]. In this case, there are $2^5-1=31$ possible
choices
(combinations).
For D = 100, the number of combinations is approximately $10^{30}$, and the full-search
method is considered difficult.
# Initialize random seed values
seed_everything()
# Size of input values (the number of pseudo-materials)
D = 100
# Property tables
debye_temperature_table, state_density_table, interaction_table = set_properties(D)
# Random searches: evaluate the supercon_temp() function for n_cycle times with random input x and output the obtained maximum and average critical temperatures.
n_cycle = 100
t_max = 0.0 # Variable to store the maximum value of critical temperature.
t_mean = 0.0 # Variable to calculate the average value of the critical temperature
for i in range(n_cycle):
x = np.random.randint(0, 2, D)
if np.sum(x) == 0:
continue
t_c = supercon_temperature(
x, debye_temperature_table, state_density_table, interaction_table
)
if t_max < t_c:
t_max = t_c
t_mean += t_c
t_mean /= n_cycle
print(f"Max. critical temperature: {t_max:.2f} K")
print(f"Mean critical temperature: {t_mean:.2f} K")
print(f"{n_cycle=}")
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. 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.client import FixstarsClient
client = FixstarsClient()
client.parameters.timeout = 1000 # Timeout 1s
# client.token = "xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx" # If you use Amplify in a local environment, enter the Amplify API token.
2.2. Implementing FM with PyTorch¶
Here, FM is implemented with PyTorch. In the TorchFM class, we define the acquisition
function $g(x)$ as a machine learning model.
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 right-hand 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
import copy
def train(
X,
y,
model_class=None,
model_params=None,
batch_size=1024,
epochs=3000,
criterion=None,
optimizer_class=None,
opt_params=None,
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)
if lr_sche_class is not None:
scheduler = lr_sche_class(optimizer, **lr_sche_params)
best_score = 1e18
for epoch 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 lr_sche_class is not None:
scheduler.step()
with torch.no_grad():
model.load_state_dict(best_model_wts)
model.eval()
return model
2.3. Construction of initial training data¶
The gen_training_data function evaluates the objective function $f(x)$ against the input
value $x$ to produce $N_0$ input-output 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.
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.4. 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().
from amplify import (
Solver,
BinarySymbolGenerator,
sum_poly,
BinaryMatrix,
BinaryQuadraticModel,
)
import matplotlib.pyplot as plt
import sys
class FMQA:
def __init__(self, D: int, N: int, N0: int, k: int, true_func, solver) -> None:
assert N0 < N
self.D = D
self.N = N
self.N0 = N0
self.k = k
self.true_func = true_func
self.solver = solver
self.y = None
# A member function that repeatedly performs (N-N0)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 input-output pair [x_hat, y_hat] to the training data set
X = np.vstack((X, x_hat))
y = np.append(y, y_hat)
# Copy the input-output 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 = BinarySymbolGenerator() # Declare a variable generator, BinaryPoly
q = gen.array(self.D) # Generate decision variables using BinaryPoly
cost = self.__FM_as_QUBO(
q, w0, w, v
) # Define FM as a QUBO equation from FM parameters
result = self.solver.solve(
cost
) # Pass the objective function to Amplify solver
if len(result.solutions) == 0:
raise RuntimeError("No solution was found.")
values = result.solutions[0].values
q_values = q.decode(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)
D = w.shape[0]
out_1 = sum_poly(self.k, lambda i: sum_poly(D, lambda j: x[j] * v[j, i]) ** 2)
# Note that x[j] = x[j]^2 since x[j] is a binary variable in the following equation.
out_2 = sum_poly(
self.k, lambda i: sum_poly(D, lambda j: x[j] * v[j, i] * v[j, i])
)
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 non-diagonal 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 i-th 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("i-th evaluation of f(x)", fontsize=18)
plt.ylabel("f(x)", fontsize=18)
plt.tick_params(labelsize=18)
return fig
3. FMQA execution example¶
3.1. Material search for the highest critical temperature¶
Now we will perform a material search using the implemented FMQA and model functions. Since we are maximizing the critical temperature estimated from the model, we implement the objective function to return the negative value of the critical temperature and perform the FMQA to minimize this value.
In the following, $N = 100$ and $N_0=60$. Thus, in the example below, $N-N_0=40$ cycles of FMQA (machine learning, search for the optimal solution with the quantum annealing or Ising machines, and evaluation of the objective function) are performed. Note that with this setup, it will take approximately 5-10 minutes to complete all FMQA cycles.
# Initialize random seed values
seed_everything()
# Size of input values (the number of pseudo-materials)
D = 100
# Property tables
debye_temperature_table, state_density_table, interaction_table = set_properties(D)
# Objective Function. We implement the objective function so that it returns a negative value of the critical temperature, and FMQA optimizes the choice of material to minimize this value
def true_func(x):
if np.sum(x) == 0:
return 0
return -supercon_temperature(
x, debye_temperature_table, state_density_table, interaction_table
)
N = 100 # Number of times the function can be evaluated
N0 = 60 # Number of samples of initial training data
k = 20 # Dimension of the vector in FM (hyperparameters)
# client: Amplify client created earlier
solver = Solver(client)
# Generate initial training data
X, y = gen_training_data(D, N0, true_func)
# Instantiate FMQA class
fmqa_solver = FMQA(D, N, N0, k, true_func, solver)
# Run FMQA cycle
pred_x = fmqa_solver.cycle(X, y)
# Output optimization results
print("pred x:", pred_x)
print("pred value:", true_func(pred_x))
3.2. Transition of objective function values during the FMQA optimization process¶
Plotted below are the $N_0$ objective function values obtained for randomly generated input values during initial training data generation and the evolution of objective function values during the FMQA optimization process for $N-N_0$ cycles.
The blue and red lines, respectively, show how the smallest objective function value is successively updated from the input values obtained by the FMQA optimization cycle (red line).
In general, due to the principle of the heuristics algorithm employed in FixstarsClient,
the
solutions obtained are not reproducible, but for the material choices obtained with the parameters
$N_0=60$ and $N-N_0=40$ in the sample code, the resulting critical temperature is approximately 50 K.
fig = fmqa_solver.plot_history()
3.3. Example output from the FMQA sample program¶
In general, due to the principle of the heuristics algorithm employed in FixstarsClient,
the
solutions obtained are not completely reproducible, but typical standard output and image output
obtained
when this sample code is used are shown below.
-
The following is a typical output obtained when this sample code is executed. When the FMQA code in Material search for the highest critical temperature" is executed under the given conditions, the following standard output is sequentially output as the FMQA cycle progresses.
Generating 0-th training data set. Generating 10-th training data set. Generating 20-th training data set. Generating 30-th training data set. Generating 40-th training data set. Generating 50-th training data set. Starting FMQA cycles... FMQA Cycle #0 variable updated, pred_y=-18.98476017536205 FMQA Cycle #1 FMQA Cycle #2 variable updated, pred_y=-25.897204545387414 FMQA Cycle #3 variable updated, pred_y=-30.641568733824826 FMQA Cycle #4 FMQA Cycle #5 variable updated, pred_y=-33.23380829087865 FMQA Cycle #6 FMQA Cycle #7 FMQA Cycle #8 variable updated, pred_y=-40.97929639761995 FMQA Cycle #9 FMQA Cycle #10 FMQA Cycle #11 FMQA Cycle #12 FMQA Cycle #13 FMQA Cycle #14 FMQA Cycle #15 FMQA Cycle #16 FMQA Cycle #17 FMQA Cycle #18 variable updated, pred_y=-42.00895340350797 FMQA Cycle #19 variable updated, pred_y=-47.787495086366945 FMQA Cycle #20 FMQA Cycle #21 variable updated, pred_y=-52.41427395241357 FMQA Cycle #22 FMQA Cycle #23 FMQA Cycle #24 FMQA Cycle #25 FMQA Cycle #26 FMQA Cycle #27 FMQA Cycle #28 FMQA Cycle #29 FMQA Cycle #30 FMQA Cycle #31 FMQA Cycle #32 FMQA Cycle #33 FMQA Cycle #34 FMQA Cycle #35 FMQA Cycle #36 FMQA Cycle #37 FMQA Cycle #38 variable updated, pred_y=-55.425491086604936 FMQA Cycle #39 pred x: [0. 0. 0. 1. 1. 0. 1. 0. 0. 1. 1. 1. 0. 0. 1. 0. 1. 1. 1. 1. 1. 0. 0. 1. 1. 1. 0. 0. 1. 1. 1. 1. 0. 1. 1. 1. 1. 0. 1. 0. 0. 0. 0. 1. 0. 0. 0. 0. 2. 0. 0. 0. 1. 1. 0. 1. 1. 1. 0. 1. 1. 1. 1. 1. 1. 0. 1. 1. 0. 0. 0. 1. 3. 1. 1. 1. 0. 0. 0. 1. 0. 1. 0. 1. 0. 0. 0. 0. 0. 1. 1. 0. 1. 1. 0. 0. 4. 0. 1. 1.] pred value: -55.425491086604936
-
The output image from
fmqa_reactor.plot_history()as described in "3.2. Transition of objective function values during FMQA optimization process" is as follows:
Exercises¶
Exercise 1¶
Using the random search approach, how many attempts will it take to find a combination of materials
that
achieves the same level of critical temperature as that obtained by the black box optimization? You
can
check this by varying the n_cycle in section 1.3.
-
Supplemental information
If $\tau_{ML}$ is the time required for one machine learning for FM, $\tau_{QA}$ is the time required to find the optimal solution, and $\tau_{eval}$ is the time required to evaluate the objective function, in general, the total time cost $c_t$ for the search can be described by:
$$ c_t = N_0 \cdot \tau_{eval} + (N - N_0) \cdot (\tau_{ML} + \tau_{QA} + \tau_{eval} ). $$
In this example code, $\tau_{eval}$ is relatively small because the model is used to evaluate the critical temperature. However, in general, for tasks that require black-box optimization, $\tau_{eval} \gg \tau_{ML}$ and $\tau_{eval} \gg _{QA}$. In that case, the total time cost is:
$$ c_t \sim N \cdot \tau_{eval}. $$
In the present case, for example, if one hour is needed for each material synthesis + critical temperature measurement, and if this is done independently 24/7, approximately 4 days are required for a search for $N=100$, and 1 year for $N=10000$. Therefore, to keep the optimization cost small, keeping the number of evaluations of the objective function $N$ small is a priority.
It will be shown that random searches generally do not approach or exceed the optimal solution by FMQA unless unrealistically-large $N$ is used (i.e., $c_t$ is enormous)
Exercise 2¶
Let us vary the hyperparameters related to the machine learning of FM. How will the optimization accuracy and computation time change?
- Hint: try to change the parameters in the model call in the
step()function of theFMQAclass excerpted below. (e.g., change the number of epochsepochto 1/10 of the original value)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}, )
Exercise 3¶
Think of a business or a problem in your immediate environment where black-box optimization can be utilized. What are the decision variables (input values) and the objective function? How would the objective function be evaluated?
Example: Optimization of fuel blends (hydrogen, natural gas, syngas, ammonia, steam, recirculated and flue gases, etc.) and thermochemical conditions in a next-generation gas turbine power plant. Reduce fuel procurement costs and pollutant generation while ensuring power output to meet daily demand. The objective function is the amount of pollutant generation, fuel cost, and $($electricity demand$-$power output$)^2$, and its evaluation is based on a (linear) calculation of fuel cost and full-scale simulation or measurement on an actual plant.
