Processing device searching for an optimum value of a cost function by using dynamics

[0002] 1. Field of the Invention

[0003] The present invention relates to a processing device for searching for an optimum value of a cost function in an optimization problem such as an optimum structure problem, an optimum location problem, an optimum route problem, etc.

[0004] 2. Description of the Related Art

[0005] In recent years, the demand for solving an optimization problem has been increasing in various industrial fields. The optimization problem is a problem of searching for a state where a given cost function becomes a maximum, a minimum, a local maximum, or a local minimum. By changing the sign of the cost function, the problem of obtaining a maximum or a local maximum is changed to a problem of obtaining a minimum or a local minimum. Hereinafter, the optimization problem will be described mainly as a problem of obtaining a minimum or a local minimum.

[0006] For example, in a CAD (Computer-Aided Design), an evaluation value of intensity or stability is used as a cost function in order to improve the intensity of a building, a construction, etc., and the stability to an external force, and the structure corresponding to its optimum value is then obtained. Additionally, in a material design, the energy of an atomic/molecular level of a material is used as a cost function, and the structure corresponding to the lowest energy state is obtained. Also in a more general parameter fitting problem, a parameter combination which makes a cost optimum is obtained.

[0007] As the conventional methods for solving such optimization problems, there are two methods such as a descent method and a probabilistic method. A steepest descent method being a representative of the descent method is a method which calculates from a given state a direction where the value of a cost function drops, changes the state in the calculated direction, and obtains one candidate of the minimum value of the cost function. In this way, the local minimum value of the cost function can be obtained.

[0008] Additionally, algorithms representative of the probabilistic method include a random method, an SA (Simulated Annealing) method, and a GA (Genetic Algorithm) method.

[0009] The random method is a method which selects a state at random, and picks up a state where the value of a cost function is small. The SA method adopts a metropolis algorithm so as to determine a next state. According to this algorithm, if the cost of the next state drops, that state is adopted. If the cost rises, the state is adopted with a certain probability. This probability depends on a temperature parameter. The metropolis algorithm is used by setting a temperature somewhat high at first. Thereafter, the temperature is gradually lowered.

[0010] The GA method is an optimization method which simulates the mechanism of the evolution of living things. This method represents a state as a character string which is referred to as a chromosome, and optimizes respective genes by performing genetic operations such as selection, crossing-over, mutation, etc. for a chromosome aggregation.

[0011] However, the above described conventional methods have the following problems.

[0012] With the steepest descent method, if the state where a cost function becomes a local minimum is found during a search, the search falls into a trap and cannot get out of that state. Accordingly, the state where the cost function is the minimum cannot be necessarily found. The random method has a possibility that a strict solution can be found if the number of states is finite and small. However, if the number of states becomes large, this method does not work well.

[0013] With the SA method, a temperature must be lowered in order to improve a realization probability of a low cost value state according to the metropolis algorithm, which leads to a slow state change for a search. Accordingly, once the search is trapped into a local minimum state, it is impossible to get out of this state without waiting for a long time. Accordingly, the method for setting a temperature high at first, and for lowering the temperature gradually is considered to be effective. However, there is no deterministic and general implementation of this temperature scheduling method, and how to arrange a temperature schedule becomes a problem.

[0014] Furthermore, since the GA method is lack of a neighborhood search capability, it cannot find a more optimum state if it exists in a neighborhood, which frequently causes a near miss. Still further, a theoretical proof that an optimum state can be found is insufficient.

[0015] An object of the present invention is to provide a processing device for efficiently searching for an optimum value of a cost function regardless of its local form.

[0016] In a first aspect of the present invention, a processing device comprises an inputting device, a candidate searching device, and an outputting device.

[0017] The inputting device inputs a state where a problem is described, and information of a cost function giving the cost of the state. The candidate searching device obtains one or more candidates of an optimum state by using deterministic dynamics which generates a realization probability that takes a higher value at a state that attains a more optimum value of the cost function. The outputting device outputs the obtained candidates.

[0018] In a second aspect of the present invention, a processing device comprises an inputting device, a defining device, a setting device, a calculating device, and an outputting device.

[0019] The inputting device inputs a parameter. The defining device defines a cost function which gives a state cost. The setting device sets an adjustment a function. The calculating device numerically solves a differential equation, which is described by the parameter, the cost function, and the adjustment function, in terms of a value, and obtains one or more candidates of an optimum state. The outputting device makes an obtained calculation result visible, and presents it to a user.

[0020] FIG. 1 is a block diagram showing the principle of a processing device according to the present invention;

[0021] FIG. 2 shows the configuration of the processing device;

[0022] FIG. 3 shows input data;

[0023] FIG. 4 shows output data;

[0024] FIG. 5 is a flowchart showing the process performed by an inputting unit;

[0025] FIG. 6 is a flowchart showing the process performed by a cost defining unit;

[0026] FIG. 7 is a flowchart showing the process performed by a cost function value calculating unit;

[0027] FIG. 8 is a flowchart showing the process performed by a function generating unit;

[0028] FIG. 9 is a flowchart showing the process performed by an equation calculating unit;

[0029] FIG. 10 is a flowchart showing the processes performed by a numerical integration performing unit and a checking unit;

[0030] FIG. 11 is a flowchart showing the process performed by a candidate selecting unit;

[0031] FIG. 12 is a flowchart showing the process performed by a descending unit;

[0032] FIG. 13 shows a cost function;

[0033] FIG. 14 shows a coordinate value distribution;

[0034] FIG. 15 shows the realization frequency of a coordinate value;

[0035] FIG. 16 shows the realization frequency of a cost function value;

[0036] FIG. 17 is a block diagram showing the configuration of an information processing device;

[0037] FIG. 18 shows storage media;

[0038] FIG. 19 exemplifies a particle arrangement;

[0039] FIG. 20 shows a potential function;

[0040] FIG. 21 shows an initial arrangement;

[0041] FIG. 22 shows a first simulation result;

[0042] FIG. 23 shows a second simulation result;

[0043] FIG. 24 shows a third simulation result;

[0044] FIG. 25 shows the evaluation of the first simulation result;

[0045] FIG. 26 shows the evaluation of the third simulation result;

[0046] FIG. 27 shows first program code;

[0047] FIG. 28 shows second program code;

[0048] FIG. 29 shows third program code;

[0049] FIG. 30 shows fourth program code;

[0050] FIG. 31 shows fifth program code;

[0051] FIG. 32 shows sixth program code;

[0052] FIG. 33 shows seventh program code;

[0053] FIG. 34 shows eighth program code;

[0054] FIG. 35 shows ninth program code;

[0055] FIG. 36 shows tenth program code;

[0056] FIG. 37 shows eleventh program code;

[0057] FIG. 38 shows twelfth program code;

[0058] FIG. 39 shows thirteenth program code;

[0059] FIG. 40 shows fourteenth program code;

[0060] FIG. 41 shows fifteenth program code;

[0061] FIG. 42 shows sixteenth program code;

[0062] FIG. 43 shows seventeenth program code;

[0063] FIG. 44 shows eighteenth program code;

[0064] FIG. 45 shows nineteenth program code;

[0065] FIG. 46 shows twentieth program code;

[0066] FIG. 47 shows twenty-first program code;

[0067] FIG. 48 shows twenty-second program code;

[0068] FIG. 49 shows twenty-third program code;

[0069] FIG. 50 shows twenty-fourth program code;

[0070] FIG. 51 shows twenty-fifth program code;

[0071] FIG. 52 shows twenty-sixth program code; and

[0072] FIG. 53 shows twenty-seventh program code.

[0073] Preferred embodiments according to the present invention will be hereinafter described in detail, by referring to the drawings.

[0074] FIG. 1 is a block diagram showing the principle of a processing device according to the present invention. The processing device shown in this figure comprises an inputting device 1, a candidate searching device 2, and an outputting device 3.

[0075] The inputting device 1 inputs a state where a problem is described, and the information of a cost function giving the cost of the state. The candidate searching device 2 obtains one or more candidates of an optimum state by using deterministic dynamics which generates a realization probability that takes a higher value at a state that attains a more optimum value of the cost function. The outputting device 3 outputs the obtained candidates.

[0076] The candidate searching device 2 determines a set of variables representing a state and the calculation algorithm of dynamics by using the information input by the inputting device 1, and makes calculations according to the set of variables and the calculation algorithm. Then, the candidate searching device 2 obtains a set of variable values which will become an optimum state candidate corresponding to the optimum value of the cost function. The outputting device 3 outputs candidates obtained by the candidate searching device 2 as a search result.

[0077] Dynamics corresponds to time development of a coordinate (point) generated by a calculation based on an equation, etc. By performing the dynamics in a state space where a state is defined by coordinate variables, the whole of a given state space can be searched under an appropriate condition, so that the search is never trapped by a local minimum of a cost function, etc. Accordingly, calculations can be made regardless of a local form of a cost function.

[0078] Additionally, with the dynamics which generates a realization probability that takes a higher value at a state that attains a more optimum value of a cost function, the neighborhood of the optimum states can be searched well without decreasing a system temperature. Accordingly, there is no need to retard a state change in order to generates the realization probability that takes a higher value at a state that attains a more optimum value of the cost function unlike the conventional probabilistic method, thereby increasing processing efficiency.

[0079] As described above, the point of the present invention is to search for an optimum state by using the deterministic dynamics which generates a realization probability that takes a higher value at a state that attains a more optimum value of a cost function.

[0080] The processing device shown in FIG. 1 further comprises a storing device 4 and an optimum state searching device 5. The storing device 4 stores one or more candidates obtained by the candidate searching device 2. The optimum state searching device 5 obtains a state closer to the optimum state by changing the respective states of the stored candidates in the direction where the value of a cost function is improved.

[0081] The calculation made by the optimum state searching device 5 corresponds to, for example, a calculation of a descent method, and can obtain the state of a cost value closer to the, optimum value than the state of a candidate. By obtaining the state with the best cost value among the states thus obtained, a solution with high precision can be obtained.

[0082] For example, the inputting device 1 shown in FIG. 1 corresponds to an inputting unit 10 and an inputting device 53, which are respectively shown in FIGS. 2 and 17 and will be described later, and the candidate searching device 2 shown in FIG. 1 corresponds to a cost function value calculating unit 12, an equation calculating unit 14, a numerical integration performing unit 15, a frequency calculating unit 17, and a candidate selecting unit 18, which are shown in FIG. 2. Additionally, for example, the outputting device 3 shown in FIG. 1 corresponds to an outputting unit 19 and an outputting device 54, which are respectively shown in FIGS. 2 and 17, the storing device 4 corresponds to a memory 52 shown in FIG. 17, and the optimum state searching device 5 shown in FIG. 1 corresponds to a descending unit 20 shown in FIG. 2.

[0083] In this preferred embodiment, the minimum value of a cost function U(q) which is determined by a state q=(q₁, . . . , q_n) represented by “n” real variables, and its state q_min(or a value closer to the minimum value and its state) are searched. Note that U(q) is defined to be a function which is differentiable. By changing the sign of the cost function U(q), the problem of searching for the minimum value can be changed to the problem of searching for the maximum value. Here, the calculation having the following natures is aimed to be implemented.

[0084] (A) Having a neighborhood search capability.

[0085] (B) Not being trapped even if a local minimum state is found during a search.

[0086] (C) Enabling the realization probability of the state of a low cost value to be improved by setting parameters and adjustment functions. Furthermore, a search is prevented from slowing down like the SA method.

[0087] (D) Being given a theoretical proof that an optimum state can be found, covering how the realization probability of each state reaches after sufficiently long time and which conditions are required.

[0088] To implement such aims, several candidates of q_minare first obtained by using deterministic dynamics. Then, q_minor a state that takes a value closer to the minimum value is obtained by changing the states of the respective candidates with a descent method.

[0089] Here, the following ordinary differential equations are used as an example of the deterministic dynamics with which the candidates of q_minare obtained.

[0090] dq_i/dt=−(β/n)τ_K(p)p_i+( β/n)[B_i(x)∂Θ_Z(ζ,w)/∂w−∂B_i(x)/∂w](i=1, . . . , n)  (1)

[0091] dp_i/dt=(β/n)τ_U(q)F_i(q)−(β/n)[(1/Q)∂Θ_Z(ζ, w)/∂ζ+α∂Θ_Z(ζ,w)/∂w]p_i(i=1, . . . , n)  (2)

[0092] dζ/dt=[−(β/n)τ_K(p)K(p)−β]/Q+(β/n)[A(x)∂Θ_Z(ζ,w)/∂w−∂A(x)/∂w]  (3)

[0093] dw/dt=[−(β/n)τ_K(p)K(p)−β]α−(β/n)[A(x)∂Θ_Z(ζ,w)/∂ζ−∂A(x)/∂Θ]+(β/n)Σ[τ_U(q)F_i(q)B_i(x)+∂B_i(x)/∂q_i]  (4)

[0094] Notice that the following definitions are used here. $\Sigma \equiv \sum _{i=1}^n$  q≡(q₁, . . . , q_n)εD⊂Rⁿ  (5)

[0095] p≡(p₁, . . . , p_n)εRⁿ  (6)

[0096] x≡(q, p, ζ, w)εΓ≡D×Rⁿ×R²  (7)

[0097] U:Rⁿ⊃D→R  (8)

[0098] F_i(q)≡−∂U(q)/∂q_i  (9)

[0099] τ_U(q)≡[dΘ_U(u)/du]_u=U(q)  (10)

[0100] K(p)≡Σp_i²  (11)

[0101] τ_K(p)≡31 2[dΘ_K(k)/dk]_k=K(p)  (12)

[0102] Θ_U:R⊃D_U≡U(D)→R  (13)

[0103] Θ_K:[0,∞)→R  (14)

[0104] Θ_Z:R²→R  (15)

[0105] A:Γ→R  (16)

[0106] B=(B₁, . . . ,B_n):Γ→Rⁿ  (17)

[0107] The ordinary differential equations (1) through (4) correspond to equations extended from Newton equations. It is supposed that “q_i” is a coordinate and, “p_i” corresponds to a momentum. Additionally, “ζ” and “w” are introduced as extended variables. “Θ_U”, “Θ_K”, “Θ_Z”, “A”, and “B” are smooth functions which are introduced as adjustment functions, while “β”, “Q”, and “α” are introduced as adjustment parameters. “β” serves as a set temperature. These adjustment functions and parameters can be arbitrarily set. At this time, the equations (1) through (4) have the following characteristics.

[0108] The left side of the equation (2) represents a time change in momentum, and “F_i(q)” in the first term on the right side represents a force obtained by attaching a negative sign to the differentiation of a cost function. In other words, the direction of a state change (direction of an acceleration) includes the component of the differentiation of the cost function with the reverse sign. Accordingly, if the cost function increases, the state changes in the direction reverse thereto. If the cost function decreases, the state changes in that direction. As described above, the equation (2) represents the nature in which if a local minimum exists, an approach is made to that direction, and has the nature described in (A).

[0109] Furthermore, the second term on the right side of the equation (2) represents a frictional force which is proportional to “p_i”, and the coefficient of “p_i” described by Θ_Z(ζ,w) corresponds to a frictional coefficient. The left side of the equation (3) represents a time change of “ζ”, the first term—(β/n)τ_K(p)K(p) in the first term on the right side represents system temperature, and K(p) represents the kinetic energy of the system.

[0110] Here, assume that ∂Θ_Z(ζ,w)/∂ζ is an increasing function of “ζ”, the system temperature exceeds the set temperature “β” and the time change of “ζ” becomes positive. In this case, the frictional force on the right side of the equation (2) increases and the momentum decreases, so that the system temperature drops. Inversely assume that the system temperature becomes lower than the set temperature “β”, and the time change of “ζ” becomes negative. In this case, the frictional force on the right side of the equation (2) decreases and the momentum increases, so that the system temperature rises. The frictional force in the equation (2) may be of a positive or a negative value.

[0111] Accordingly, the equations (2) and (3) represent the nature which attempts to make the system temperature approach the set temperature “β”. Actually, it can be proved that the time average of the system temperature becomes “β” under a suitable condition. Therefore, if thermal vibrations are given to the system by making an adjustment such as a setting of “β” to a large value, it becomes possible to get out of a local minimum state. As described above, the equations (2) and (3) have the above described nature (B).

[0112] Furthermore, if some conditions which are considered to be satisfied in normal cases are satisfied, the probability that U(q) takes a value between “u1” and “u2” after a long time (theoretically, an infinite time) elapses can be represented by the adjustment function Θ_Uwhich generates τ_U(q) in the equation (2). To be more specific, if the system has ergodicity, it can be proved that this probability is given by the following equation under an additional mathematical condition. $\begin{array}{cc}\underset{\tau \to \infty }{\lim }\frac{1}{\tau }{\int }_0^{\tau }{\chi }_S\left(T_t\left(x\right)\right)dt={\int }_{\mathrm{u1}}^{\mathrm{u2}}ku\left(u\right)\Omega \left(u\right)du& \left(18\right)\end{array}$  S≡{qεD¦u1≦U(q)≦u2}×Rⁿ⁺²  (19)

[0113] k_U(u)=exp[−Θ_u(u)]  (20)

[0114] Here, the left side of the equation (18) represents a limiting value of the definite integral with respect to time “t”, while its right side represents the definite integral with respect to the value “u” of the cost function U(q). T_t(x) on the left side represents the time development (flow) of a coordinate “x”, and χ_s(x) is a function which becomes “1”if “x” is an element of a set “S”, and becomes “0” if “x” is not an element of the set “S”.

[0115] Additionally, Ω on the right side is a function which represents the density of the number of states of the cost function. k_U(u)Ω(u) represents the probability density with which the cost value U(q) of the system becomes “u”. The probability represented by the equation (18) represents the ratio of the time period during which the system remains in a state where U(q) takes a value in a range between “u1” and “u2”. This probability may be referred to also as a visit frequency of an orbit to that range.

[0116] The equation (18) is proved from the fact that the system has an invariant measure given by the density function ρ(x) in the following equation under a suitable condition.

[0117] ρ(x)≡exp[−{Θ_U(U(q))+Θ_K(K(p))+Θ_Z(ζ, w)}]  (21)

[0118] where ρ(x) represents the probability density with which the state of x=(q, p, ζ, w) is realized.

[0119] To implement the above described nature (C), “k_U” may be set so that “k_UΩ” in the equation (18) becomes a maximum in the neighborhood of the minimum value of “u”. The density function Ω rapidly increases with the increase of “u” in normal cases. Therefore, if “k_U” is suitably set, for example, as a function which rapidly decreases with the increase of “u”, it is possible to make the peak of “k_UΩ” approach the minimum value of “u”. To make the function “k_U” rapidly decrease, Θ_Umay be made to rapidly increase. As a result, the visit frequency of the state of a low cost value can be improved after a long time elapses.

[0120] Because the function “k_U” can be set independently of the set temperature “β” which determines the proceeding speed (that is, the search speed) of the system, it is theoretically guaranteed to prevent the search speed from degrading while increasing the realization probability of a state of a low cost value. Accordingly, the search speed does not slow down as a result of dropping the temperature in order to increase the realization probability of the state of a low cost value, unlike the conventional SA method. As described above, the equation (2) has the above described natures (C) and (D).

[0121] In other words, examples of the ordinary differential equations having the invariant measure which realizes the probability represented by the equation (18) are the equation (1) through (4). The equations (1) through (4) have an advantage to flexibly cope with an individually given problem by suitably setting the adjustment parameters and functions. However, the dynamics which realizes the probability represented by the equation (18) is not limited to the equations (1) through (4), and may be implemented as a different equation. For example, the number of extended variables to be introduced is not required to be two (“ζ” and “w”) in all cases.

[0122] Explained next is the processing device which searches for an optimum value of a cost function by using the equations (1) through (4). This processing device is configured, for example, by using a computer. This device solves the equations (1) through (4) with a suitable numerical integration method, and obtains candidates of q_min. Normally, however, convergence into an optimum solution is not guaranteed. Therefore, the processing device terminates a search by setting a suitable termination condition. Then, the device descends the cost value with a suitable descent method by using the respectively obtained candidates of q_minas initial values, and obtains a more optimum state.

[0123] FIG. 2 shows the configuration of such a processing device. The processing device shown in this figure comprises an inputting unit 10, a cost defining unit 11, a cost function value calculating unit 12, a function generating unit 13, an equation calculating unit 14, a numerical integration performing unit 15, a checking unit 16, a frequency calculating unit 17, a candidate selecting unit 18, an outputting unit 19, and a descending unit 20.

[0124] The inputting unit 10 inputs input data 21. The cost defining unit 11 sets a cost function and its partial derivative. The cost function value calculating unit 12 calculates the value of the function set by the cost defining unit 11 at a time “t”. The function generating unit 13 generates a new adjustment function depending on need. The equation calculating unit 14 performs the calculations of the equations (1) through (4) at the time “t” by using the information from the cost function value calculating unit 12 and the function generating unit 13. The numerical integration performing unit 15 performs numerical integration by using the results of the calculations made by the equation calculating unit 14. The checking unit 16 checks the result of the integration.

[0125] The frequency calculating unit 17 calculates the realization frequencies of a cost function “U”, kinetic energy “K” (temperature), a coordinate value “q_obs” to be observed, etc. The candidate selecting unit 18 selects candidates of an optimum state “q_min” based on these frequencies. The outputting unit 19 outputs as output data 22 the result of the integration and that of the calculations made by the frequency calculating unit 17 to a file.

[0126] The numerical integration performing unit 15 checks whether or not a termination condition is satisfied. If the termination condition is not satisfied, the numerical integration performing unit 15 advances the time “t” by Δt, and performs the next integration. Thereafter, if the termination condition is satisfied, the numerical integration performing unit 15 terminates the integration.

[0127] Then, the output data 22 is made visible on a display screen, and the plurality of obtained optimum state candidates 23 are passed to the descending unit 20. The descending unit 20 searches for “q_min” with a descent method based on input data 21 and the optimum state candidates 23, and outputs the obtained state as an optimum state 24. Also the optimum state 24 can be made visible on the display screen.

[0128] In FIG. 2, the numerical integration performing unit 15, the checking unit 16, the frequency calculating unit 17, the candidate selecting unit 18, the outputting unit 19, and the descending unit 20 have universal capabilities which do not depend on a given problem.

[0129] FIG. 3 shows the input data 21. In this input data 21, parameters 31 are parameters required for defining a cost function, and degree of freedom 32 is the number “n” of state variables which describe a given problem.

[0130] Additionally, simulation conditions 33 include the number of steps, a time interval, termination conditions, an output specification parameter, etc. The number of steps represents the number of times that the numerical integration and the descent method are repeated. The time interval represents an interval Δt of the numerical integration. The termination conditions represent the termination conditions of the numerical integration and the descent method. The output specification parameter is a parameter which specifies the output interval of the output data 22, etc. As the termination conditions, for example, the following conditions are used.

[0131] (a) Terminating the calculation when the calculation time or the number of processing steps reaches a predetermined value.

[0132] (b) Terminating the calculation when states of a cost value which is lower than a predetermined target value of U(q), the number of which is predetermined number or more, are obtained.

[0133] Furthermore, frequency calculation parameters 34 include a discretization interval, a variable scope parameter, etc. The discretization interval represents an interval between values when the realization frequency of the value of a given variable or function is calculated. The variable scope parameter represents the calculation range of the value of the variable or function.

[0134] Still further, adjustment parameters 35 are the values of the above described parameters “β”, “Q” and “β”. Selection conditions of the adjustment functions 36 are the conditions for setting the above described functions “Θ_U”, “Θ_K”, Θ_Z, “A”, and “B”. Various adjustment functions are stored as built-in functions in the processing device beforehand. If the identification number of a desired function is input as a selection condition, the specified built-in function is automatically set. Or, if the definition of a new function is input as a selection condition, the new adjustment function is set.

[0135] A boundary condition 37 represents a boundary condition for a domain “D” of a cost function. For example, if a torus is specified as a boundary condition, the processing device recognizes the area “D” to be successive in a torus state, and performs the numerical integration.

[0136] FIG. 4 shows the output data 22. In this output data 22, a time change of a state variable 41 represents a change of a variable “q”, and a time change of other variable values 41 represent a change of variables other than “q”. A time change of an optimum cost value 43 represents a change of an optimum value of a cost value obtained by a search.

[0137] Additionally, a frequency of a cost function value 44 represents a realization frequency of the cost value, which is totaled for each discretization interval. A frequency of a system temperature 45 represents a realization frequency of the temperature, which is totaled for each discretization interval. A frequency of a coordinate to be observed 46 represents a realization frequency of a coordinate value “q_obs”, which is totaled for each discretization interval.

[0138] Next, the processes performed by the processing device shown in FIG. 2 are explained in detail by referring to FIGS. 5 through 12.

[0139] FIG. 5 is a flowchart showing the process performed by the inputting unit 10. The inputting unit 10 first inputs input data 21 (step S11), performs the process for defining the initial values of state variables and other variables (step S12), and prepares for a frequency calculation (step S13).

[0140] In the initial value definition process, the inputting unit 10 inquires of a user as to whether or not to perform automatic generation (step S14). If the automatic generation is instructed, the inputting unit 10 generates the initial values of respective variables with a predetermined method (step S15), and terminates the process. If the automatic generation is not instructed, the inputting unit 10 reads the initial values from a predetermined external file, assigns the read values to the respective variables (step S16), and terminates the process.

[0141] FIG. 6 is a flowchart showing the process performed by the cost defining unit 11. The cost defining unit 11 first defines a cost function and its partial derivative based on the parameters 31 (step S21), and inquires of a user as to whether or not to modify the cost function (step S22). If modification is instructed, the cost defining unit 11 modifies the cost function (step S23), and terminates the process. If the modification is not instructed, the cost defining unit 11 terminates the process without modifying the cost function. As an example of cost function modification, the operation for excluding the part of a relatively large value of a defined cost function from a search target is considered.

[0142] FIG. 7 is a flowchart showing the process performed by the cost function value calculating unit 12. The cost function value calculating unit 12 receives a cost function from the cost defining unit 11, and calculates the value of the cost function based on an updated state “q” at a time “t” (step S31). Then, the cost function value calculating unit 12 calculates the value of the partial derivative of the cost function (step S32), and terminates the process.

[0143] FIG. 8 is a flowchart showing the process performed by the function generating unit 13. The function generating unit 13 first determines whether or not to newly generate an adjustment function based on the selection conditions of adjustment functions 36 (step S41). If the selection conditions 36 instructs new generation, the function generating unit 13 generates an adjustment function based on input information (step S42), and terminates the process. If the new generation is not instructed, the function generating unit 13 terminates the process without generating an adjustment function.

[0144] FIG. 9 is a flowchart showing the process performed by the equation calculating unit 14. The equation calculating unit 14 calculates the variables such as a temperature, etc., which are required for calculating the equations (1) through (4) (step S51), calculates the right sides of the equations (1) through (4) by using the calculation results (step S52), and terminates the process.

[0145] FIG. 10 is a flowchart showing the processes performed by the numerical integration performing unit 15 and the checking unit 16. The numerical integration performing unit 15 performs numerical integration with a Runge-Kutta method, a Gear method, or a different method (step S60), and performs the process regarding the boundary condition 37 (step S61). Then, the checking unit 16 checks whether or not a numerical error occurs (step S62).

[0146] If no numerical error occurs, the checking unit 16 requests the frequency calculating unit 17 to perform the succeeding process (step S63), and terminates the process. If a numerical value error occurs, the checking unit 16 terminates the numerical integration regardless of whether or not a termination condition is satisfied (step S64), and terminates the process.

[0147] FIG. 11 is a flowchart showing the process performed by the candidate selecting unit 18. The candidate selecting unit 18 selects a plurality of optimum state candidates from among the states obtained so far based on the results of the calculations made by the frequency calculating unit 17 and by the cost function value calculating unit 12 (step S71). Then, the candidate selecting unit 18 stores the selected states as candidates of “q_min” and the corresponding cost values as candidates of “U_min” (step S72), and terminates the process.

[0148] FIG. 12 is a flowchart showing the process performed by the descending unit 20. The descending unit 20 first advances the calculation step of a predetermined descent method by recognizing a given candidate of “q_min” as an initial state (step S81), and calculates the value of a cost function (step S82). Then, the descending unit 20 outputs the calculation result to a file (step S83), and checks whether or not a termination condition is satisfied (step S84).

[0149] If the termination condition is not satisfied, the descending unit 20 repeats the operations in and after step S81. If the termination condition is satisfied, the descending unit 20 terminates the descent method (step S85). Then, the descending unit 20 makes the output data visible on a display screen (step S86), outputs the obtained optimum state (step S87), and terminates the process.

[0150] Described next is the result of simulation of the numerical integration with a 4th Runge-Kutta method by setting n=2 and representing a cost function U(q) with a sum of seven two-dimensional Gaussian functions. A search area “D” is assumed to be a two-dimensional torus, and adjustment functions are set as follows.

[0151] A=0

[0152] B=0

[0153] Θ_U(u)=(1/2T′)u²

[0154] Θ_K(k)=(1/2T)k

[0155] Θ_Z(ζ, w)=(1/2T)(Qζ²+α′w²)  (22)

[0156] Here, “T′”, “T”, and “α′” are parameters which determine the adjustment functions, and are set as follows.

[0157] T′=10.0

[0158] T=10.0

[0159] α′=0.0

[0160] Additionally, the adjustment parameters are set as follows.

[0161] β=nT=20.0

[0162] Q=0.001

[0163] α=0.0

[0164] Furthermore, the initial conditions are set to be q₁=q₂=p₂=ζ=w=0.0, and p₁=T^0.5, the number of steps of the numerical integration is set to be 10000000, and the time interval of the numerical integration is set to be 0.0001. Since the amount of data may sometimes become significant if data are output in all of the steps, the data output interval is set to be 1000 steps.

[0165] At this time, U(q) is given as a function shown in FIG. 13, and includes a plurality of local minimum values. The smallest value among these values becomes an optimum value of a cost function. After the numerical integration is performed by a specified number of steps, the coordinate value distribution shown in FIG. 14 is obtained and the realization frequency of a coordinate value becomes the one shown in FIG. 15. In this figure, the frequency becomes high at the positions corresponding to the minima values shown in FIG. 13. Additionally, the realization frequency of the cost function becomes the one shown in FIG. 16. The frequency peaks shown in FIG. 16 correspond to the local minimum values shown in FIG. 13.

[0166] In the above described preferred embodiments, the cost function is defined to be a function of “n” real variables, which is differentiable. Also in a problem described by discrete variables, a similar search can be performed if a suitable cost function is set. As discrete optimization problems, for example, an optimum location problem, an optimum route problem, an optimum network problem, an optimum flow problem, an optimum efficiency problem, etc.

[0167] The optimum location problem is a problem of optimizing the location of a facility, etc. in an urban design. The optimum route problem is a problem of optimizing a route in car navigation, an electric circuit, etc.

[0168] Additionally, the optimum network problem is a problem of optimizing gas or water pipe laying, electrical wiring, communications network, etc. The optimum flow problem is a problem of optimizing a traffic flow on a road, a data flow on a network, etc. The optimum efficiency problem is a problem of optimizing an efficiency in fields such as science, engineering, economy, business, etc.

[0169] The above described processing device can be configured by using an information processing device (computer) shown in FIG. 17. The information processing device shown in this figure comprises a CPU (Central Processing Unit) 51, a memory 52, an input device 53, an output device 54, an external storage device 55, a medium driving device 56, and a network connecting device 57, which are interconnected by a bus 58.

[0170] The memory 52 includes, for example, a ROM (Read Only Memory), a RAM (Random Access Memory), etc., and stores programs and data used for processing. The CPU 51 performs necessary processing by executing the programs with the memory 52.

[0171] The inputting unit 10, the cost defining unit 11, the cost function value calculating unit 12, the function generating unit 13, the equation calculating unit 14, the numeric integration performing unit 15, the checking unit 16, the frequency calculating unit 17, the candidate selecting unit 18, the outputting unit 19, and the descending unit 20, which are shown in FIG. 2, correspond to software components stored in specific program code segments in the memory 52, and are implemented by a program composed of one or more instructions.

[0172] The input device 53 is, for example, a keyboard, a pointing device, a touch panel, etc., and is used to input the instruction or the information from a user. The output device 54 is, for example, a display, a printer, a speaker, etc., and is used to output an inquiry to a user or a process result.

[0173] The external storage device 55 is, for example, a magnetic disk device, an optical disk device, a magneto-optical disk device, etc. The information processing device may store the above described programs and data in the external storage device 55, and can use them by loading into the memory 52 depending on need. Additionally, the external storage device 55 may be used also as a database storing cost functions, adjustment functions, etc.

[0174] The medium driving device 56 drives a portable storage medium 59, and accesses its stored contents. As the portable storage medium 59, an arbitrary computer-readable storage medium such as a memory card, a floppy disk, a CD-ROM (Compact Disk-Read Only Memory), an optical disk, a magneto-optical disk, etc. is used. The above described programs and data may be stored in the portable storage medium 59, and can be used by being loading into the memory 52 as occasion demands.

[0175] The network connecting device 57 communicates with an external device via an arbitrary network (line) such as a LAN (Local Area Network), and performs data conversion accompanying a communication. The information processing device may receive the above described programs and data from an external device via the network connecting device 57, and can use them by loading into the memory 52 depending on need.

[0176] FIG. 18 shows computer-readable storage media which can provide programs and data to the information processing device shown in FIG. 17. The programs and data stored in the portable storage medium 59 or an external database 60 are loaded into the memory 52. The CPU 51 performs necessary processing by executing the programs with the data.

[0177] Next, the result of the simulation where the search processing according to the present invention is applied to the problem of optimizing the arrangement of a multi-particle system such as molecules, atoms, etc. In this case, the particle arrangement corresponding to the minimum energy state of a particle system is searched.

[0178] Assuming that the space where N particles are arranged is an m-dimensional space, the number of degrees of freedom n becomes Nm. At this time, the positional vectors r₁, . . . , r_Nare given by the following equations.

[0179] r₁≡(q₁, . . . , q_m),

[0180] r₂≡(q_m+1, . . . , q_2m),

[0181] r_N≡(q_(N−1)m+1, . . . , q_Nm)  (23)

[0182] Additionally, assuming that the search space of each positional vector is the whole of the m-dimensional space, the Nm-dimensional arrangement space is represented as follows. $\begin{array}{cc}{R}^{\mathrm{Nm}}-\underset{i<j}{U}\left\{\left(r_1,\dots ,r_N\right)|r_i=r_j\right\}& \left(24\right)\end{array}$

[0183] For example, if N=100 and m=2, the particle arrangement shown in FIG. 19 can possibly exist. In FIG. 19, 100 particles are arranged on the two-dimensional plane. Furthermore, as the cost function U(q), the potential energy function of the particle system represented by the following equation is used. $\begin{array}{cc}U\left(q\right)=\sum _{1\le i,j\le N}\phi \left(\Vert r_i-r_j\Vert \right)& \left(25\right)\end{array}$

[0184] where ∥r_i−r_j∥ represents the distance between “i”th and “j”th particles, and a function φ(∥r_i−r_j∥) represents a pair potential. If the pair potential is set on a Lennard-Jones type, a function φ(r) is given by the following equation.

[0185] φ(r)=4ε((σ/r)¹²−(σ/r)⁶)  (26)

[0186] If “Å” is used as the unit of the distance between particles and ε=1.7258 and σ=3.446 are set, the function form of φ(r) becomes the one shown in FIG. 20. Additionally, the adjustment functions and the adjustment parameters are set as follows.

[0187] A=0

[0188] B=0

[0189] Θ_U(u)={(U_*−U_M)/aT′} {(u−U_M)/(U_*−U_M)}^a

[0190] Θ_K(k)=(1/2β)k

[0191] Θ_Z(ζ,w)=(1/2β)(Qζ²+α′w²)

[0192] T′=20.0

[0193] α′=α=0

[0194] β=nT=4000

[0195] Q=0.001  (27)

[0196] Note that, however, the following three types of combinations are used as U_*, U_M, and a. $\mathrm{C1}\text{:}\begin{array}{c}{U}_{*}=8079.0+U_{1b}\\ U_M=7000.0+U_{1b}\\ a=4.0\end{array}$ $\mathrm{C2}\text{:}\begin{array}{c}{U}_{*}=8079.0+U_{1b}\\ U_M=7000.0+U_{1b}\\ a=2.0\end{array}$ $\mathrm{C3}\text{:}\begin{array}{c}{U}_{*}=1.0+U_{1b}\\ U_M=U_{1b}\\ a=2.0\end{array}$  U_1b=−εN(N−1)/2=−8542.7  (28)

[0197] Furthermore, the numerical integration is performed with the 4th Runge-Kutta method under the condition that the particle arrangement shown in FIG. 21 is used as the initial arrangement, the number of steps of the numerical integration is 50,000 to 200,000, and the time interval of the numerical integration is 0.0001.

[0198] As a result, the particle arrangements shown in FIG. 22, 23 and 24 are obtained as the candidates of the minimum energy states in the above described cases of C1, C2, and C3. The values of the cost function in the arrangements shown in FIGS. 22, 23, and 24 are respectively −478.52, −479.67, and −463.45.

[0199] In FIGS. 22 and 24, if adjacent particles are linked by straight lines and evaluated, the results become the ones shown in FIGS. 25 and 26. Something like a lattice defect occurs as indicated by a symbol “∘” in the arrangement shown in FIG. 26. In the meantime, such a defect is not found in the arrangement in FIG. 25.

[0200] FIGS. 27 through 53 exemplify the program codes of the search processing using the cost function shown in FIG. 13. FIGS. 27 through 36 among these figures correspond to the main routine. The program code shown in FIG. 27 defines the variables used for the processing. The program code shown in FIG. 28 defines input data shown in FIG. 3. The program code shown in FIG. 29 defines constants.

[0201] Here, NFREE, NSTEP, DT, NOUT, NUMFUNC, and IPBC respectively correspond to the degree of freedom “n”, the number of steps, the time interval, the output specification parameter, the number of a built-in function, and the boundary condition. Here, n is equal to 2.

[0202] Additionally, TEMP, QMASS, ALPHA1, ALPHA2, TEMP2, and BETA respectively correspond to the above described parameters T, Q, α, α′, T′, and β. Also UEX is given as one of the adjustment parameters. Here, the functions represented by the equations (22) are used as the adjustment functions.

[0203] Furthermore, Q(1), Q(2), P(1), P(2), Z(1), and Z(2) respectively correspond to the above described variables q₁, q₂, p₁, p₂, ζ, and w. UMIN, DOSUU(*,1), DOSUU(*,2), and DOSUU2(*,*,1) respectively correspond to the optimum cost value, the frequency of the cost function value, the frequency of the system temperature, and the frequency of the coordinate to be observed shown in FIG. 4.

[0204] The program code shown in FIG. 30 corresponds to the process in step S12 shown in FIG. 5, and defines the initial values of variables. The program code shown in FIG. 31 corresponds to the process in step S13 shown in FIG. 5, and prepares the histogram for a frequency calculation. The program code shown in FIG. 32 records input data. Here, DR1 and DR2 corresponds to the discretization interval shown in FIG. 3. UINI, UFIN, P21NI, and P2FIN corresponds to the variable scope parameters shown in FIG. 3. The program code shown in FIG. 33 corresponds to the loop process of the numerical integration shown in FIG. 2, and calls a subroutine AMDYN.

[0205] The program code shown in FIG. 34 makes output data visible, and outputs frequency data. The program code shown in FIG. 35 further outputs a cost function graph and its partial differentiation graph. The program code shown in FIG. 36 calls a subroutine DESCENT of a descent method, and terminates the process.

[0206] FIGS. 37 through 43 show the subroutine AMDYN of the numerical integration, which is called in FIG. 33. The program code shown in FIG. 37 defines the variables used for the process. The program code shown in FIG. 38 corresponds to the process in step S60 shown in FIG. 10, and performs the numerical integration by using a subroutine EQ1 with the 4th Runge-Kutta method. The program code shown in FIG. 39 corresponds to the process in step S61 shown in FIG. 10, and processes a boundary condition. Here, a torus is used as the boundary condition. However, a different boundary condition may be specified. The program code shown in FIG. 40 calculates the values of variables.

[0207] The program code shown in FIG. 41 corresponds to the process in step S63 shown in FIG. 10, and calculates frequencies. The program code shown in FIG. 42 corresponds to the processes in steps S71 and S72, and selects an optimum state candidate. Here, the state where the cost function becomes a minimum is selected as an optimum state. Alternatively, also a semi-optimum state, etc. may be selected. The program code shown in FIG. 43 outputs a calculation result as output data.

[0208] FIGS. 44 through 47 show the subroutine EQ1 of an equation calculation, which is called in FIG. 38. The program code shown in FIG. 44 defines the variables used for the process. The program code shown in FIG. 45 calls a subroutine COSTFUNC for calculating the value of a cost function. The program code shown in FIG. 46 corresponds to the process in step S51 shown in FIG. 9, and calculates a system temperature. The program code shown in FIG. 47 corresponds to the process in step S52 shown in FIG. 9, and calculates the right side of the equation.

[0209] FIGS. 48 through 52 show the subroutine COSTFUNC for calculating the value of a cost function, which is called in FIGS. 37, 40, and 45. This subroutine corresponds to the process of FIG. 7. First of all, the program code shown in FIG. 48 defines the variables used for the process. The program code shown in FIG. 49 specifies a built-in function.

[0210] In this way, the value of the cost function corresponding to the program code shown in FIG. 50 or 51 is calculated. Here, since NUMFUNC=1, the cost function shown in FIG. 50 is used. The program code shown in FIG. 52 records input data.

[0211] FIG. 53 shows a subroutine DESCENT of the descent process, which is called in FIG. 36. This program code corresponds to the process in FIG. 12, and obtains an optimum state with the descent method.

[0212] According to the present invention, optimum state candidates can be sequentially searched by using the capability for local searching and the capability for avoiding a trap into a local minimum state in the processing for searching the state corresponding to the optimum value of a cost function. These two capabilities may be implemented by suitably setting adjustment functions, not in accordance with the nature of a force function or a temperature control.

[0213] Additionally, the realization probability of the state of a low cost value can be improved not directly dependent on the value of the temperature parameter for controlling a search speed. Therefore, the temperature parameter can be increased to a value such that no fault is caused by a numerical value calculation error, whereby the processing speed promises to be improved.