Simulation device and simulation method

Abstract

To provide a simulation device that can efficiently arrive at a good solution by appropriately determining a range of a temperature parameter for searching a ground state of an interaction model using simulated annealing (SA).SOLUTION: A simulation device, which uses SA (simulated annealing) to find a solution to a combinatorial optimization problem, acquires a local optimal solution to an object function thereof by a local search method that changes a temperature parameter, calculates a statistic from the local optimum solution and a probability distribution function that a variable included in the object function follows, calculates, on the basis of the calculated statistic (standard deviation of the local optimum state), a maximum value (temperature parameter with a large standard deviation) and a minimum value (temperature parameter with a small standard deviation) that are a search range of the temperature parameter, and performs, in the search range of maximum and minimum values of the temperature parameter, a plurality of times of simulated annealing to find a good solution.SELECTED DRAWING: Figure 10

Description

[Technical field]

[0001] The present invention relates to a simulation device and a simulation method, and more particularly to a simulation device and a simulation method suitable for use in setting a temperature parameter when solving a combinatorial optimization problem by a simulated annealing technique. [Background technology]

[0002] Various physical and social phenomena can be expressed using interaction models. An interaction model is defined by multiple nodes that make up the model, nonlinear coefficients between the nodes (interaction coefficients), and, if necessary, linear coefficients that act on each node (bias coefficients and external magnetic field coefficients). In the fields of physics and social science, various models including the Ising model have been proposed, and all of them can be interpreted as a form of interaction model.

[0003] Finding the node state (called the ground state) that minimizes the energy of this interaction model can be applied to solving combinatorial optimization problems. Examples of combinatorial optimization problems include the problem of detecting cliques in social networks and portfolio optimization problems in the financial field. In the field of operations research, these are broadly classified as unconstrained binary quadratic optimization problems and mixed (variables contain a mixture of integers and real numbers) binary quadratic programming problems.

[0004] Simulated Annealing (hereinafter abbreviated as "SA") is an algorithm for solving the ground state search problem of an interaction model. This is a probabilistic method that gradually reduces the randomness of state updates, similar to the annealing of physical phenomena (for example, in metal processing, a method of stabilizing the shape of a metal by gradually cooling it after heating it).

[0005] In simulated annealing of physical phenomena, thermal fluctuations are gradually weakened, causing the object to converge to a more energetically stable state.

[0006] Similarly, simulated annealing optimization uses temperature as an external parameter, which is initially set to a high value and gradually lowered as the optimization progresses. At high temperatures, the algorithm is willing to accept a solution that worsens the objective function compared to the current solution with a certain probability, thus avoiding getting stuck in a local optimum and being unable to find a global optimum. Also, as the temperature decreases, the probability of accepting a solution that worsens the objective function decreases, narrowing the search region in which the optimal solution is located. This gradual cooling process allows simulated annealing to find a good solution effectively.

[0007] Regarding simulated annealing, a system for solving mixed binary quadratic programming problems including the ground state search problem of the Ising model is disclosed, for example, in Patent Document 1. Patent Document 1 describes a method for efficiently performing ground state search based on simulated annealing by converting the interaction relationship into a complete bipartite graph structure for mixed binary quadratic programming problems including the ground state search problem of the Ising model.

[0008] Furthermore, another technique for searching the ground state of an Ising model based on simulated annealing is disclosed, for example, in Patent Document 2. Patent Document 2 describes a method for determining the maximum value of the temperature parameter of simulated annealing by comparing the energy of a local optimum solution state and a state obtained by inverting the spin therefrom in relation to the ground state search of an Ising model based on simulated annealing. [Prior art documents] [Patent documents]

[0009] [Patent Document 1] International Publication No. 2021 / 220445 [Patent Document 2] JP 2020-173579 A Summary of the Invention [Problem to be solved by the invention]

[0010] As mentioned above, in simulated annealing, the hyperparameter that determines the strength of randomness is called the temperature parameter. In order to maximize the performance of simulated annealing, it is necessary to appropriately adjust the randomness, i.e., the temperature parameter.

[0011] If the range of the temperature parameter is too wide, the amount of calculation required to gradually lower the temperature parameter will increase, reducing the efficiency of the simulation.On the other hand, if the range of the temperature parameter is too narrow, there is a high possibility that a global optimum cannot be reached.

[0012] Patent Document 2 discloses a technology in which, in simulated annealing, a temperature correction unit calculates a correction temperature based on the maximum increase in the energy value corresponding to each of a plurality of state variables obtained by changing the state variables corresponding to each local solution by one bit (Fig. 10, paragraph

[0117] , Fig. 11).

[0013] However, although there is a description regarding a method for calculating the upper limit of the temperature parameter, there is no description regarding a method for calculating the lower limit.

[0014] In addition, there is no description of a method for determining the range of temperature parameters for a model in which the variables are multi-valued or continuous.

[0015] The object of the present invention is to provide a simulation device and a simulation method that enable efficient acquisition of good solutions by appropriately defining the range of temperature parameters for performing simulated annealing to search for the ground state of an interaction model, including cases where variables are multi-valued or continuous. [Means for solving the problem]

[0016] The configuration of the simulation device of the present invention is a simulation device that solves a combinatorial optimization problem by simulated annealing, in which a local optimal solution to the problem is obtained by a local search method, statistics are calculated from the local optimal solution and the probability distribution function to which the variables included in the objective function obey, and the maximum and minimum values ​​that form the search range of the temperature parameter are found based on the calculated statistics. Simulated annealing is then performed multiple times within the search range of the maximum and minimum values ​​of the temperature parameter, thereby efficiently finding a better solution. Effect of the Invention

[0017] According to the present invention, it is possible to provide a simulation device and a simulation method that enable a better solution to be obtained efficiently by appropriately defining the range of the temperature parameter for solving the ground state search problem of an interaction model by simulated annealing. [Brief description of the drawings]

[0018] [Figure 1] FIG. 1 is a schematic diagram of the energy landscape of the interaction model. [Diagram 2] FIG. 13 is a diagram showing the progress of a solution search when the temperature parameter T is gradually changed in SA. [Diagram 3] FIG. 1 is a diagram showing energy transitions between a state corresponding to a locally optimal solution of a ground state search problem of the Ising model and a state obtained by inverting the spin from the locally optimal solution. [Figure 4] FIG. 13 is a diagram showing the discrete probability distribution function that a binary variable follows depending on the value taken by a temperature parameter. [Diagram 5] FIG. 13 is a diagram showing the discrete probability distribution function that a continuous variable follows depending on the value of a temperature parameter. [Figure 6] FIG. 13 is a diagram showing state transitions during SA depending on the value of a temperature parameter. [Figure 7] FIG. 2 is a diagram illustrating a local search for obtaining a local optimum solution, thereby obtaining a plurality of local optimum solutions. [Figure 8] FIG. 2 is a functional configuration diagram of the simulation device. [Figure 9] FIG. 2 is a diagram showing the hardware and software configuration of the simulation device. [Figure 10] 13 is a flowchart showing a ground state search process. [Figure 11] 13 is a flowchart showing SA processing performed by the simulation device. [Figure 12] FIG. 2 is a data flow diagram in the simulation device. [Figure 13] FIG. 13 is a diagram showing an example of a simulation result screen. DETAILED DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0019] Hereinafter, an embodiment of the present invention will be described with reference to the drawings. The examples are illustrative for explaining the present invention, and are omitted and simplified as appropriate for clarity of explanation. The present invention can be implemented in various other forms. Unless otherwise specified, each component may be singular or plural.

[0020] In order to facilitate understanding of the invention, the position, size, shape, range, etc. of each component shown in the drawings may not represent the actual position, size, shape, range, etc. Therefore, the present invention is not necessarily limited to the position, size, shape, range, etc. disclosed in the drawings.

[0021] As examples of various information, the various information may be described using expressions such as "table," "list," and "queue," but the various information may be expressed using data structures other than these. For example, various information such as "XX table," "XX list," and "XX queue" may be expressed as "XX information." When describing identification information, expressions such as "identification information," "identifier," "name," "ID," and "number" are used, but these are mutually interchangeable.

[0022] When there are multiple components having the same or similar functions, they may be described by using the same reference numerals with different subscripts, or when there is no need to distinguish between these multiple components, the subscripts may be omitted.

[0023] In the embodiments, a process performed by executing a program may be described. Here, a computer executes the program using a processor (e.g., CPU, GPU), and performs the process defined by the program while using a storage resource (e.g., memory) and an interface device (e.g., communication port), etc. Therefore, the subject of the process performed by executing the program may be the processor. Similarly, the subject of the process performed by executing the program may be a controller, device, system, computer, or node having a processor. The subject of the process performed by executing the program may be a calculation unit, and may include a dedicated circuit that performs specific processing. Here, the dedicated circuit is, for example, an FPGA (Field Programmable Gate Array), an ASIC (Application Specific Integrated Circuit), or a CPLD (Complex Programmable Logic Device), etc.

[0024] The program may be installed in the computer from a program source. The program source may be, for example, a program distribution server or a computer-readable storage medium. When the program source is a program distribution server, the program distribution server may include a processor and a storage resource that stores the program to be distributed, and the processor of the program distribution server may distribute the program to be distributed to other computers. In addition, in the embodiment, two or more programs may be realized as one program, and one program may be realized as two or more programs.

[0025] First, the concept and processing that are the premise of the simulation of the present invention will be outlined.

[0026] The simulation device of this embodiment is a device that performs a simulation for solving a problem using an interaction model.

[0027] The interaction model is defined by multiple nodes that make up the model, nonlinear coefficients acting between the nodes, and linear coefficients acting on each node. Here, the variable x corresponding to each node i (i = 1 to N) i is a continuous variable x i ∈[-1,1] or a binary variable x i ∈{-1,1}. Based on the nonlinear and linear coefficients of the interaction model, the energy function H(x) (generally called Hamiltonian) is defined. In particular, when an interaction relationship is defined between two nodes, the energy function is a quadratic equation shown in the following (Equation 1).

[0028]

number

[0029] Here, J is a square matrix representing the parameters of the interactions between the nodes, and h is a vector representing the bias on the nodes (also called the external magnetic field).

[0030] In general, the interaction model is represented as an undirected graph, and the nonlinear terms are specified by the set of variables contained therein, so that the matrix J is a real symmetric matrix.

[0031] The interaction model is a concept that includes the Ising model. In the interaction model, all variables are treated as binary variables x i When we restrict it to ∈{-1,1}, we reduce to the Ising model, which is used, for example, in statistical mechanics as a lattice model to describe magnetic materials, with +1 / -1 corresponding to the up / down orientation of the spin.

[0032] In this embodiment, we consider the case where the above-mentioned continuous variables or binary variables correspond to each node, but by considering an expanded dimension of variables that take on multiple values ​​or continuous values ​​as variables corresponding to each node, models that appear in physics such as the Potts model, the XY model, and the Heisenberg model can also be regarded as interaction models.

[0033] The ground state search of the interaction model is an optimization problem to find the variable array x that minimizes the above energy function H(x). This is equivalent to a mixed binary quadratic programming problem with interval constraints on the variables.

[0034] Here, the energy transition of the interaction model, problems caused by the conventional technology, and the relationship between the statistics indicated by the SA model and the temperature parameters will be described with reference to FIG.

[0035] The horizontal axis of the graph is the variable array, and the vertical axis is the total energy of the system. In Figure 1, state A has the smallest energy, and updating -1 from state A to 1 results in state B with higher energy, and updating another 1 to -1 results in state C with lower energy.

[0036] In this embodiment, the ground state of the interaction model is searched for by Markov Chain Monte Carlo methods (hereinafter referred to as MCMC). In the stochastic transition in MCMC, a stochastic transition from the current state x to another state is repeated. There is a method of sequentially executing the stochastic state transition, called Gibbs sampling (also called the Heatbath method).

[0037] Let p(x) be the probability distribution function that state x follows, and let x be the state generated by the αth Gibbs sampling. (α) ={x1 (α) ,…,x N (α) In this case, the conditional probability of the α+1th state is generated by sequentially updating the variables with subscripts i=1, ..., N according to the following (Equation 2).

[0038]

number

[0039] When MCMC is run by introducing a parameter T that controls the above stochastic state transitions and gradually suppressing the transitions, it asymptotically converges to the ground state. Simulated annealing is a method that utilizes this to find the optimal solution to the minimization problem of the objective function.

[0040] Below, we will explain the transition of the energy correspondence in the state space when T is gradually changed in SA.

[0041] Corresponding to real-world simulated annealing, the parameter T is often called the temperature parameter. The larger the temperature parameter, the greater the randomness of the state transitions, and the smaller the temperature parameter, the more the state transitions behave like local search. In the following, we will refer to the maximum value of the temperature parameter as T H , the minimum value is T L As a method for scheduling the temperature parameters of SA, for example, T H From T LA geometric cooling method is commonly used, in which the temperature parameter is lowered by stepping between these values, which are interpolated using a geometric progression.

[0042] When applying MCMC or SA to an interaction model, the variable values ​​are probabilistically updated based on (Equation 2) or the like. Here, variables that do not have a nonlinear term between them are independent in the energy function formula and therefore the probability density function, and it is possible to simultaneously apply probabilistic state transitions based on (Equation 2). Therefore, by updating independent variables in parallel, it is possible to speed up the processing of MCMC or SA.

[0043] Such speedup can be achieved by using an interaction model in which the interaction relationship has a bipartite graph structure, i.e., each subgraph does not have an interaction relationship. For example, Patent Document 1 proposes a method for efficiently finding a ground state by solving the ground state search problem of an interaction model in which the interaction relationship has a complete bipartite graph structure after relaxing all variables to continuous values ​​in the ground state search problem of a fully connected interaction model.

[0044] As shown in the energy transition diagram in Figure 2, since the temperature parameter determines how the state space is explored, it is important to set the temperature parameter appropriately in order to maximize the performance of SA. H A method for determining T is shown. Figure 3 shows the change in the ground state and the energy transition at that time. T=T H In this case, it is desirable to explore the state space by transitioning to other states without remaining in the local optimum solution. Therefore, in the technology described in Patent Document 2, a local optimum solution is first obtained, and the energy difference with the state where the spin is inverted is calculated, and a temperature parameter T that exceeds the energy difference is set. H It is suggested to set it to .

[0045] However, in an interaction model that includes variables with multiple next-state candidates (for example, continuous variables as dealt with in Patent Document 1), multiple candidates for referring to the energy difference as shown in Figure 3 appear, making it impossible to estimate appropriate temperature parameters using the method of Patent Document 2.

[0046] In view of this background, the simulation method of this embodiment determines appropriate temperature parameters for simulated annealing for an interaction model having multi-valued or continuous-valued variables by referring to information on local optimum solutions and statistics of the probability distribution functions obeyed by the variables.

[0047] Here, the diagonal element J in the energy function of Eq. ii Here is an example for the case where x is >0. i If is a binary variable {-1,+1}, then J ii The term containing J ii x i 2 =J ii is a constant, and its contribution can be neglected.

[0048] In the simulation, the states are transitioned according to Gibbs sampling. By performing a sufficient number of state updates, the state {x} of the interaction model on the right-hand side of (Equation 1) appears with the probability distribution function p(x) based on the Boltzmann distribution shown in (Equation 3) below.

[0049]

number

[0050] Here, Z is the normalization factor of the Boltzmann distribution, which is called the partition function. By applying the Gibbs sampling shown in (Equation 2) to this Boltzmann distribution, it is possible to stochastically update x.

[0051] In the following, we explain the process of setting appropriate temperature parameters for simulated annealing based on the information of the local optimum and the statistics of the probability distribution function that the variables follow.

[0052] For this purpose, we define a vector quantity g shown in the following (Equation 4). This is related to the quantity obtained by differentiating the energy function H(x) shown in (Equation 1) with respect to x.

[0053]

number

[0054] Here, for each variable x i We assume that the probability distribution function of (Equation 3) is followed, and the standard deviation is used as an index to evaluate the degree of dispersion of the value. If the average value of a quantity · under the probability distribution function p(x) is E[·], then the standard deviation SD[x i ] is generally expressed as follows (Equation 5).

[0055]

number

[0056] First, the variable x i The case where is a binary variable will be described. From (Equation 3), the next state x i The probability distribution function for ∈{-1, +1} follows the discrete probability distribution function expressed by the following (Equation 6).

[0057]

number

[0058] If the standard deviation is calculated according to (Equation 5) for the probability distribution function of (Equation 6), it becomes as shown in the following (Equation 7).

[0059]

number

[0060] Figure 4 shows the discrete probability distribution function that the binary variable follows depending on the value of the temperature parameter.

[0061] When the temperature parameter is large (|g i / T|<<1), -1 and +1 appear with almost equal probability. In this case, the variable x i The average value of E[x i ] is almost 0 and the standard deviation SD[x i ] is almost 1.

[0062] On the other hand, when the temperature parameter is small (|g i / T|>>1) has g i Depending on the sign of x, it can be localized to -1 or +1. i The average value of E[x i ] approaches a localized value, and the standard deviation SD[x i ] is almost 0.

[0063] Next, the variable x i We will explain the case where is a continuous variable. From (Equation 3), the variable x i follows the continuous probability distribution function expressed by the following (Equation 8).

[0064]

number

[0065] This probability distribution function has a mean value μ=g defined on the closed interval [-1,+1]. i / J ii ,variance σ 2 =T / J ii This is a truncated normal distribution.

[0066] In this case, (a, b)=((-1-μ) / σ, (1-μ) / σ) is set, and the parameter indicating the range is replaced as shown in the following (Equation 9).

[0067]

number

[0068] This gives the standard deviation SD[x i ] is shown in the following (Equation 10).

[0069]

number

[0070] Figure 5 shows the discrete probability distribution function that the continuous variable follows depending on the value of the temperature parameter.

[0071] When the temperature parameter is large, the truncated normal distribution approaches the uniform distribution, and the x i The average value of E[x i ] is 0, and the standard deviation SD[x i ] approaches 1 / √3.

[0072] On the other hand, when the temperature parameter is small, |g i / J ii If |≧1, sign(g i / J ii ) to |g i / J ii If |<1, then g i / J ii to variable x i Since it is localized, the standard deviation SD[x i ] is almost always 0. Here, sign(·) is the sign function that returns the sign of the argument ·.

[0073] Based on the above, the concept of setting the range of the temperature parameter in SA in this embodiment will be described with reference to FIGS.

[0074] For interaction models that follow the probability distribution functions of (Eq. 6) and (Eq. 8), it is possible to randomly sample the next state from these distribution functions by using the inverse function method or the rejection sampling method.

[0075] Figure 6 shows how the state space moves by repeating random sampling of the next state according to (Equation 8) when variables x1 and x2 are continuous variables. When the temperature parameter is large (top of Figure 6), the state space moves significantly regardless of the energy value of the interaction model. On the other hand, when the temperature parameter is small (bottom of Figure 6), the energy value remains around a specific local optimum where the energy value is a minimum value.

[0076] In Figure 6, the darker the color, the lower the energy of the interaction model (the more stable it is). * The locus of is the search path for the value as a local optimum solution. The upper figure shows the situation when the temperature parameter is large, and the lower figure shows the situation when the temperature parameter is small.

[0077] T=T immediately after SA starts H In this case, it is desirable to have a state transition that does not remain at a local optimum. This corresponds to a situation where the probability distribution of a variable is spread, for example, the standard deviation is close to 1 to some extent, as shown in Figures 4 and 5.

[0078] On the other hand, when SA is about to end, T=T L In this case, it can be interpreted that it is desirable for there to be almost no state transitions from the local optimum and for the current solution to converge as the optimum solution. This corresponds to a situation in which the spread of the probability distribution of the variables, for example the standard deviation, is somewhat close to 0, as shown in Figures 4 and 5.

[0079] Therefore, in this embodiment, the maximum value T of the temperature parameter range is calculated in the following manner. H and the minimum value T L The aim is to be able to perform SA with appropriate temperature parameters by determining θ H and θL are the thresholds for the standard deviation when determining the maximum and minimum values ​​of the temperature parameters, respectively.

[0080] First, multiple local optimal solutions are generated using a local search method or SA with a reduced number of steps in MCMC, as shown in Figure 7. Next, the variables x i The spread of the distribution, standard deviation SD[x i ]. Then, calculate all SD[x i ] for SD[x i ]≧θ H The temperature at which this holds true is T=T H Then, in that temperature region, all variables will change state and not remain in the local optimum. Also, all SD[x i ] for SD[x i ]≦θ L The temperature at which the following holds is T=T L If this is the case, all the variables in that temperature region will barely move from the local optimum and the search will end.

[0081] All SD[x i ] for SD[x i ]=θ H If there are multiple temperatures for which this holds true, the largest value among them will be used as the maximum value of the range for the next SA.

[0082] On the other hand, all SD[x i ] for SD[x i ]=θ L If there are multiple temperatures for which the above holds, the smallest value among them will be used as the minimum value in the range for the next SA.

[0083] Here, all SD[x i ] was considered, but the SD[x i ] may be extremely large or small, so SD[x i ] outliers may be excluded.

[0084] In the above, the statistics are iBased on the standard deviation of H , minimum value T L In addition, the local optimum solution was updated with the state update. i Using statistics such as the expected absolute difference of H Or T L may be calculated.

[0085] Also, x before and after one state update i Instead of the absolute difference of H is determined so that the total travel distance in a given number of searches for the local optimum solution becomes small. L may be determined.

[0086] Next, the configuration of the simulation device of this embodiment and the data flow in the simulation will be described with reference to FIGS.

[0087] First, the functional configuration of the simulation device will be described with reference to FIG. As shown in FIG. 8, the simulation device 100 of this embodiment includes, as functional components, a model coefficient setting unit 110, an SA execution unit 120, a temperature parameter adjustment unit 130, and a storage unit 140.

[0088] The model coefficient setting unit 110 is a functional unit that sets coefficients such as a Hamiltonian that expresses a model, based on interaction model data 151 (described later).

[0089] The SA execution unit 120 is a functional unit that performs SA calculations. The SA execution unit 120 includes, as sub-functional units, a variable value initialization unit 121, a temperature parameter control unit 122, a state update control unit 123, and a variable value storage unit .

[0090] The variable value initialization unit 121 is a functional unit that initializes variables (such as the state values ​​of each node of the interaction model) used in the simulation.

[0091] The temperature parameter control unit 122 is a functional unit that controls the temperature parameter T in the SA.

[0092] The state update control unit 123 is a functional unit that executes calculations related to state update according to SA for the interaction model.

[0093] The variable value storage unit 124 is a functional unit that, when SA is completed, reads out the variable values ​​corresponding to the optimal solution and stores the read out values ​​in the auxiliary storage device 210.

[0094] The temperature parameter adjustment unit 130 adjusts the maximum value T H , minimum value T L This is a functional unit that calculates the value of

[0095] The storage unit 140 is a functional unit that stores data used by the simulation device 100 .

[0096] The storage unit 140 holds interaction model data 151, probability distribution function data 152, temperature parameter adjustment setting data 160, temperature parameter data 171, SA setting data other than temperature 172, and local optimum solution data 180. Details of each data will be described later.

[0097] Next, the hardware and software configuration of the simulation device will be described with reference to FIG. The simulation device 100 includes a processor 201 , a main memory device 202 , an auxiliary memory device 203 , an input device 204 , an output device 205 , a communication device 206 , and one or more simulation calculation devices 220 connected via a system bus 5 .

[0098] The simulation device 100 may be, for example, a general information processing device such as a PC, or may be realized in whole or in part by using virtual information processing resources such as a cloud server provided by a cloud system. The simulation device 100 may also be realized by, for example, a plurality of information processing devices connected to each other so as to be able to communicate with each other and to operate in cooperation with each other.

[0099] The processor 201 is, for example, a central processing unit (CPU) or a micro processing unit (MPU), and controls the entire simulation device 100 , refers to data in the main storage device 202 , and executes programs loaded into the main storage device 202 .

[0100] The main memory device 202 is a volatile memory device that stores programs and data, and is, for example, a semiconductor memory device such as a dynamic random access memory (DRAM).

[0101] The input device 203 is a device that accepts information input from a user, and is, for example, a keyboard, a mouse, a card reader, a touch panel, or the like.

[0102] The output device 204 is a device that outputs information to a user, and is, for example, a display device (LCD (Liquid Crystal Display), graphic card, etc.) that visualizes various information, an audio output device (speaker), a printer, etc.

[0103] The communication device 205 is a communication interface device that communicates with other devices, and is, for example, a network interface card (NIC), a wireless communication module, a universal serial interface (USB) module, a serial communication module, or the like.

[0104] The simulation calculation device 206 is a device that performs simulation calculation to execute a ground state search. The simulation calculation device 206 may take the form of an expansion card that is attached to the simulation device 100, such as a GPU (Graphics Processing Unit). The simulation calculation device 206 is configured by hardware such as a CMOS (Complementary Metal Oxide Semiconductor) circuit, an FPGA (Field Programmable Gate Array), and an ASIC (Application Specific Integrated Circuit). The simulation calculation device 206 includes a control device, a storage device, an interface for connecting to the system bus 5, and transmits and receives commands and information to and from the processor 201 via the system bus 5. The simulation calculation device 206 may be communicably connected to another simulation calculation device 206 via a communication line, for example, and may operate in cooperation with the other simulation calculation device 206. The functions realized by the simulation calculation device 206 may be realized by, for example, having the processor 201 execute a program.

[0105] The auxiliary storage device 210 is a non-volatile large-capacity storage device such as a hard disk drive, a solid state drive (SSD), etc. The programs and data stored in the auxiliary storage device 210 are loaded into the main storage device 202 as needed.

[0106] The auxiliary storage device 210 has installed therein a model coefficient setting program 211, a variable value initialization program 212, a temperature parameter control program 213, a state update control program 214, a variable value reading program 215, a temperature parameter adjustment program 216, and a simulation calculation device control program 220.

[0107] The model coefficient setting program 211, the variable value initialization program 212, the temperature parameter control program 213, the state update control program 214, the variable value reading program 215, and the temperature parameter adjustment program 216 are programs which respectively execute the functions of the model coefficient setting unit 110, the variable value initialization unit 121, the temperature parameter control unit 122, the state update control unit 123, the variable value storage unit 124, and the temperature parameter adjustment unit 130.

[0108] The simulation calculation device control program 220 is a program by which the processor 201 controls the simulation calculation device 200 .

[0109] In addition to the above programs, the auxiliary storage device 210 of the simulation device 100 also has installed therein programs such as an operating system, a file system, device drivers, and a DBMS (DataBase Management System), but these are omitted in FIG.

[0110] The auxiliary storage device 210 also stores interaction model data 151, probability distribution function data 152, temperature parameter adjustment setting data 160, temperature parameter data 171, SA setting data other than temperature 172, and local optimum solution data 180.

[0111] Furthermore, an optical storage device (such as a CD (Compact Disc) or a DVD (Digital Versatile Disc)) may be connected as the auxiliary storage device 210 so that external data can be read and saved.

[0112] Next, the processing performed by the simulation apparatus of this embodiment will be described with reference to FIGS.

[0113] The ground state search process shown in FIG. 10 is a process performed by the simulation device 100 to perform a ground state search of the interaction model by SA.

[0114] The ground state search process is started, for example, by receiving an instruction from a user via the input device 203 .

[0115] First, the model coefficient setting unit 110 sets information on the model coefficients in the simulation calculation device 206 (S101). The values ​​can also be set or edited by the user via the input device 203, for example.

[0116] Next, the variable value initialization unit 121 initializes the variable values ​​of the simulation calculation device 206, and obtains information on the local optimum solution according to an instruction from the temperature parameter adjustment unit 130 (S102). In addition, according to the obtained local optimum solution and information on the probability distribution function of the variable used in SA, T H and T L Determine (S103) Next, if necessary, the variable value initialization unit 121 reinitializes the variable values ​​of the simulation calculation device 206, and the temperature parameter control unit 122 and the state update control unit 123 give instructions to the simulation calculation device 206 to execute the SA processing (S104). Details of the SA processing will be described later with reference to FIG.

[0117] Next, the temperature parameter adjustment unit 130 determines the range of the temperature parameter (T H and T L ) is required (S105), and if it is determined that readjustment is required (S105: YES), the process returns to S103, and if it is determined that readjustment is not required (S105: NO), the process ends.

[0118] The temperature parameter adjustment unit 130 may use, for example, T H and T L and the newly obtained local optimum solution information is added to calculate T H and T L If there is a change, the SA is requested to perform the calculation again.

[0119] Next, the SA processing performed by the simulation device will be described with reference to FIG. This process corresponds to S104 in FIG.

[0120] First, the variable value initialization unit 121 initializes the variable values ​​of the simulation calculation device 206 (S201). For example, the initial value of the variables is selected by selecting one of the values ​​{-1, +1} with a probability of 50% for binary variables and by random sampling from a uniform distribution in the interval [-1, +1] for continuous variables.

[0121] Next, the temperature parameter control unit 122 updates the temperature parameters of the simulation calculation device 206 (S202). In the first update, the maximum value T H In the subsequent updates, the temperature parameter is set to the minimum value T L We will update this page accordingly.

[0122] The state update control unit 123 updates the state of the variables by random sampling based on the probability distribution function of (Equation 6) or (Equation 8) (S203). This process is performed for the variable x i The subscripts i in the formula are sequentially run from 1 to N, or the subscripts i in the formula are randomly run from 1 to N.

[0123] Next, the arithmetic unit program 902 judges whether or not the SA end condition is satisfied (for example, whether or not the state update is executed while changing the temperature parameter T a predetermined number of times) (S204). If it is judged that the SA end condition is satisfied (S204: YES), the process proceeds to S205. On the other hand, if it is judged that the stop condition is not satisfied (S204: NO), the process returns to S202.

[0124] Next, the variable value storage unit 124 reads out the value stored in the variable memory, and stores it in the auxiliary storage device 210 as the result of the ground state search (S205), and the process ends.

[0125] Next, the data flow in the simulation device will be described with reference to FIG. In FIG. 12, only the main functional components are shown in relation to the data.

[0126] The simulation device 100 receives interaction model related data 150, temperature parameter adjustment setting data 160, and SA setting data 170 from a user, and executes a series of simulation processes.

[0127] Based on the interaction model related data 150 from the user, interaction model data 151 and probability distribution function data 152 are held. The interaction model data 151 is data related to the values ​​of J and h in (Equation 1) and the type of variable (binary variable / continuous variable). The probability distribution function data 152 is data specifying a probability distribution function to which the state of the interaction model follows, such as the Boltzmann distribution in (Equation 3). In the above-mentioned embodiment, the probability distribution function of (Equation 6) or (Equation 8) is set by combining (Equation 1) corresponding to the interaction model data 151 and (Equation 3) corresponding to the probability distribution function data 152.

[0128] The temperature parameter adjustment setting data 160 is based on the local optimum solution described in the setting of the temperature parameter range. H and T L For example, when searching for a local optimum solution using SA, the threshold value (θ H , θ L 10. The temperature parameter adjustment setting data 160 is input to the temperature parameter adjustment unit 130.

[0129] The temperature parameter adjusting section 130 updates the temperature parameter data 171 based on the local optimum solution data 180 .

[0130] SA setting data 170 is held as temperature parameter data 171 and non-temperature SA setting data 172. Temperature parameter data 171 is data representing T in (Equation 6) and (Equation 8). Non-temperature SA setting data 172 is data relating to the number of times a variable is updated and the method of updating a temperature parameter, and is input to SA execution unit 120.

[0131] The SA execution unit 120 performs SA based on temperature parameter data 171 and SA setting data other than temperature 172, and updates the local optimum solution data 180.

[0132] The simulation device 100 sequentially outputs temperature parameter data 171 and local optimum solution data 180 during the process shown in FIG. 10, and can display to the user the process in which an appropriate temperature is determined. This enables the user to grasp the temperature parameters that are roughly appropriate for the problem and the provisional solution to the problem before the search is completed.

[0133] Next, the user interface provided by the simulation device will be described with reference to FIG.

[0134] A user interface when SA in this embodiment is applied to the traveling salesman problem is described below. The traveling salesman problem is a problem in which a set of cities and the distances between each city are given, and the shortest route (circuit) that returns to the starting city after visiting each city exactly once is found.

[0135] The simulation result screen 400 is made up of an optimum solution display area 410 and a temperature parameter display area 420 .

[0136] In the optimum solution display area 410, a graph 411 showing the order of the locations that the salesman will visit and a total travel distance 412 are displayed.

[0137] In the temperature parameter display area 420, the SA count 422 and the temperature parameter search range transition 421 are displayed.

[0138] As described above, according to the simulation device of this embodiment, it is possible to appropriately determine the range of the temperature parameter for searching the ground state of the interaction model by simulated annealing, based on the statistics of the already obtained local optimum solution, and to efficiently search for a better solution. [Explanation of symbols]

[0139] 100: Simulation device; 110: Model coefficient setting unit; 120: SA execution unit; 121: Variable value initialization unit; 122: Temperature parameter control unit; 123: State update control unit; 124: Variable value storage unit; 130: Temperature parameter adjustment unit; 140: Memory unit; 150...interaction model related data, 151...interaction model data, 152...probability distribution function data, 160...temperature parameter adjustment setting data, 170...SA setting data, 171...temperature parameter data, 172...SA setting data other than temperature, 180...local optimum solution data, 190...problem solution

Claims

1. A simulation device that uses simulated annealing to find a solution to a combinatorial optimization problem, A local optimum solution of the objective function of the problem is obtained by a local search method that varies the temperature parameter; Calculate statistics from the probability distribution function of the local optimum and the variables included in the objective function. A simulation device which determines maximum and minimum values ​​in a search range of a temperature parameter based on the calculated statistics, and performs simulated annealing multiple times in the search range of the maximum and minimum values ​​of the temperature parameter to determine a solution.

2. 2. The simulation device according to claim 1, wherein the statistics are a standard deviation of variables when a state is updated with respect to a state of a local optimum solution, and an expected value of a travel distance when the state is updated a predetermined number of times.

3. 2. The simulation device according to claim 1, wherein the statistical quantity is a standard deviation of a variable when a state is updated with respect to a state of a local optimum solution, and a temperature when the standard deviation is greater than a predetermined threshold is set as a maximum value within a search range of the temperature parameter, and a temperature when the standard deviation is smaller than the predetermined threshold is set as a minimum value within the search range of the temperature parameter.

4. 2. The simulation apparatus according to claim 1, wherein the objective function is a function that describes a state of an interaction model having a plurality of nodes.

5. The objective function is an energy function, 2. The simulation apparatus according to claim 1, wherein the local optimum solution is obtained as a state that minimizes the energy of the energy function.

6. 2. The simulation apparatus according to claim 1, wherein a state transition is made according to a probability distribution function including a temperature parameter as a variable, and a ground state search is performed by simulated annealing.

7. 6. The simulation apparatus according to claim 5, wherein the probability distribution function is a function based on the Boltzmann distribution.

8. 6. The simulation device according to claim 5, wherein the state transition is performed by a Markov chain Monte Carlo method.

9. 2. The simulation apparatus according to claim 1, wherein the objective function is a binary variable or a continuous variable.

10. A simulation method for finding a solution to a combinatorial optimization problem by simulated annealing using a simulation device, comprising the steps of: a step of the simulation device obtaining a locally optimal solution of the objective function by a local search method in which a temperature parameter is changed; a step of the simulation device calculating statistics from the local optimum solution and a probability distribution function to which a variable included in an objective function follows; determining maximum and minimum values ​​of a search range of a temperature parameter based on the calculated statistics by the simulation device; the simulation device performing simulated annealing multiple times in a search range of maximum and minimum values ​​of the temperature parameter to find a solution; The simulation method is characterized in that the simulation device sets the statistical quantity to a standard deviation for a local optimum state, and sets a temperature when the standard deviation for the local optimum state is greater than a predetermined threshold as a maximum value within the search range of the temperature parameter, and sets a temperature when the standard deviation for the local optimum state is smaller than the predetermined threshold as a minimum value within the search range of the temperature parameter.

Images