Black-Box Optimization of Design Parameters for Mixing Device¶

Mixing devices are essential in various fields, including chemical processes and food production. While their design and operating conditions significantly impact mixing efficiency (product quality), their determination often relies on experience and trial-and-error.

This tutorial introduces an approach to finding optimal conditions for five parameters involved in the device design using Black-Box Optimization (BBO) with Ising machine utilization. While the specific meanings of these parameters are not stated, they should be considered as "design choices" that influence the performance of the mixing device. As in actual design and control, we start from a state where "it is unclear what impact each design parameter has on the device's performance".

In this tutorial, we will evaluate the mixing device's performance (in this case, the degree of concentration variation after mixing) for different parameter combinations using the simulator. Understanding the details of the physical model or the simulator's contents is unnecessary. The primary focus of black-box optimization is on efficiently finding the best design parameters without knowing the shape and structure of the objective function.

Let's try black-box optimization by utilizing Ising machines in an engineering problem setting.

For a basic knowledge of black-box optimization based on machine learning and Ising machines used in this sample code, see "Black Box Optimization with Quantum Annealing Ising Machines". For other black-box optimization examples, see here.

This sample program (optimal design parameters by black-box optimization) consists of the following:

• 1. Objective function overview
  • 1.1. Mixing simulator overview
  • 1.2. Implementing the black-box objective function
• 2. Black-box optimization implementation (integer variables)
  • 2.1. Defining the decision variable class
    • 2.1.1. What is domain wall encoding?
    • 2.1.2. Integer decision variable class IntegerVariable
  • 2.2. FM model implementation
  • 2.3. Machine learning function implementation
  • 2.4. Solver client setup
  • 2.5. Optimization using the Ising machine
  • 2.6. Generating initial training data
• 3. Optimization of design parameters
• 4. Evaluating optimization results and history
  • 4.1. Plotting results
  • 4.2. FMQA example run

---

*In this online demo and tutorial environment, the continuous execution time is limited to about 20 minutes. If you expect the execution time to exceed 20 minutes, for example, when trying optimization by changing conditions, please copy this sample program to your local environment before execution. In that case, download the following libraries related to the MAS traffic simulator as appropriate and save them in the following directory structure.

├ fmqa_5_mixing.ipynb（this program）
└ utils/
     ├ __init__.py （blank file）
     └ mixing.py

1. Objective function overview¶

1.1. Mixing simulator overview¶

In this tutorial, we use MixingSimulator as the stirring simulator. Below is a conceptual diagram of the stirrer.

The stirrer considered in the simulator has a total of five design parameters: x0, x1, x2, x3, x4. Inputting these parameters into the simulator determines the specifications of the stirrer. The initial state is set with a substance (black) added to the liquid, and stirring is simulated over a specific period based on these stirrer specifications. The result of the simulation is the degree of mixing (standard deviation of concentration).

Here, due to the constraints of the stirrer, the value ranges for these parameters are as follows:

Below are examples of how to use the MixingSimulator stirring simulator.

The simulate() method of MixingSimulator executes the stirring simulation and returns the standard deviation of the concentration c_std in the final concentration field. If the stirring is perfectly uniform, c_std will be zero. Additionally, the plot_evolution() method displays the time-series change of the concentration field.

The following simulation results show that the stirring process gradually progresses over time, and the standard deviation of the concentration distribution also decreases. Finally, a standard deviation of concentration of 0.031 is obtained, which indicates a relatively non-uniform concentration distribution.

1.2. Implementing the black-box objective function¶

Based on the simulator MixingSimulator explained above, we will implement a black-box function. The black-box function blackbox below takes five design parameters as arguments, executes a simulation for a specific duration, and returns the result (the standard deviation of the final substance concentration).

The goal of this tutorial is to find a set of design parameters that minimizes the standard deviation of the concentration after mixing.

2. Black-box optimization implementation (integer variables)¶

This section describes the program implementation of FMQA, a black-box optimization method that utilizes the Ising machine.

The processing flow of FMQA is described in this tutorial, but it is performed according to the following cycle.

In black-box optimization, including FMQA, the objective function is treated as a black box, allowing it to be generally applied without modifying the program implementation itself. Therefore, it is not always necessary to understand the program implementation of black-box optimization explained here.

2.1. Defining the decision variable class¶

Since the Ising machine used during FMQA can directly handle only binary decision variables, appropriate encoding is necessary when considering non-binary decision variables such as integers and real numbers. In this tutorial, we consider integer decision classes with domain wall encoding and address non-binary decision variables.

2.1.1. What is domain wall encoding?¶

Domain-Wall Encoding is a method for converting discrete, non-binary variables (e.g., integer variables that can take values from 0 to k) into binary variables.

Assume an integer variable $x$ takes $(k+1)$ possible values in $\{0, 1, \dots, k\}$.

In Domain-Wall Encoding, values are represented using $k$ binary variables $q_1, q_2, \dots, q_k$ as follows:

• When $x = i$, the first $i$ elements of the variable sequence $\boldsymbol{q}$ are 1, and the rest are 0.
• In other words, the position where the switch from 1 to 0 occurs (the domain wall) represents the value of $x$.

Example: When $k = 4$¶

Value of $x$	Encoded bit string $\boldsymbol{q}$
0	[0, 0, 0, 0]
1	[1, 0, 0, 0]
2	[1, 1, 0, 0]
3	[1, 1, 1, 0]
4	[1, 1, 1, 1]

Thus, the representation of integer values is uniquely determined by the "wall position (domain wall)". The above explanation applies to integer variables with a minimum value of 0. Still, it can also be extended to integer variables with a non-zero minimum value or real variables by applying a discretization technique.

2.1.2. Integer decision variable class IntegerVariable¶

The IntegerVariable class, to be implemented below, efficiently represents integer decision variables by internally using the Amplify SDK's binary variables and applying domain wall encoding.

The core functionality of the IntegerVariable class is the conversion between integer values and their binary representations.

• encode(self, x: int) -> np.ndarray: This method converts the integer value x passed as an argument into a binary vector (NumPy array) based on the corresponding domain wall encoding. For example, encoding x=3 for a variable with bounds=(0, 5) will return an array like [1., 1., 1., 0., 0.].

• decode(self, x: np.ndarray) -> int: This method takes the binary variable results (NumPy array) computed by the Ising machine and decodes them back into their original integer values. Following the domain wall encoding, it reconstructs the original integer value by counting the number of ones in the array and adding the lower bound.

Additionally, the IntegerVariable class provides the following two properties:

• constraint: This property returns the domain wall constraint considered for this integer variable. In optimization, this constraint needs to be taken into account. This ensures that the Ising machine handles the associated binary variables (below binary_variables) according to the domain wall encoding.

• binary_variables: This property returns the vector of Amplify SDK binary variables that constitute this integer variable. It can be used when you want to reference the integer variable directly as a combination of binary variables.

The implementation below also considers Variables, which collectively manages multiple IntegerVariable objects.

2.2. FM model implementation¶

We will implement a TorchFM class using PyTorch to define the FM model, similar to how a standard machine learning model is typically constructed. The FM model is a machine learning model represented by the following polynomial. Here, $\boldsymbol{x}$ represents the variables, $d$ is a constant representing the length of the input to the black-box function, $\boldsymbol{v}$, $\boldsymbol{w}$, and $w_0$ are the model coefficients (weights and biases in machine learning terms), and $k$ is a hyperparameter representing the size of the parameters.

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

2.3. Machine learning function implementation¶

Next, we will implement the function train for training the TorchFM model defined above. This process is also handled in the same way as standard machine learning. Still, as an essential performance metric for FM models in black-box optimization, we separately display the correlation coefficient between the predicted values of the trained model and the actual values.

2.4. Solver client setup¶

We will configure the Ising machine (solver client) for use in black-box optimization. In this tutorial, we use Fixstars Amplify Annealing Engine (Amplify AE). If running in a local environment, please enter the API token for Fixstars Amplify, which can be obtained for free.

2.5. Optimization using the Ising machine¶

In this section, we define a function to optimize the trained FM model above. The procedure is as follows:

1. Based on the model coefficients of the trained FM model, construct a QUBO model mathematically equivalent to the FM model using Amplify SDK's decision variables.
2. Optimize the constructed QUBO model using an Ising machine.
3. Return the optimization results.

2.6. Generating initial training data¶

We evaluate the black-box function using randomly generated input vectors. The n_0 input-output pairs obtained in this way are adopted as the initial training data.

3. Optimization of design parameters¶

We will execute black-box optimization consisting of n cycles using the functions and classes implemented so far. In the code below, we set the number of times the objective function is evaluated to n = 5. This setting is for minimal operation confirmation in this demo/tutorial environment with a limited execution time. For execution conditions and examples of actual black-box optimization, please refer to "FMQA Execution Examples".

For execution, we first encode the input vectors of the initial training data and convert them into binary values. While the FM model's training data considers the binarized x_encoded, it evaluates the black-box function using the decoded integer decision variables x_hat_decoded.

4. Evaluating optimization results and history¶

4.1. Plotting results¶

The following functions plot the transition of objective function evaluation values during the initial learning data generation process and the optimization cycle.

4.2. FMQA example run¶

In general, due to the nature of the heuristic algorithm employed by AmplifyAEClient, there is no absolute reproducibility in the obtained solutions. However, we present a typical execution result obtained when running the sample code below.

The figure below shows an example of execution when n = 50, illustrating the optimization history (the transition of black-box function values and the history of best solution updates).

While there are increases and decreases in the objective function values in each optimization cycle, on average, the design parameters that stirr more uniformly are explored as the optimization cycles progress. In this execution example, a concentration standard deviation of 0.014 was obtained as the final best solution, and the stirring process at those design parameters is shown below.