WO2023238606A1 - 評価装置、評価方法、およびプログラム - Google Patents

評価装置、評価方法、およびプログラム Download PDF

Info

Publication number
WO2023238606A1
WO2023238606A1 PCT/JP2023/018114 JP2023018114W WO2023238606A1 WO 2023238606 A1 WO2023238606 A1 WO 2023238606A1 JP 2023018114 W JP2023018114 W JP 2023018114W WO 2023238606 A1 WO2023238606 A1 WO 2023238606A1
Authority
WO
WIPO (PCT)
Prior art keywords
characteristic
data
control
evaluation value
point
Prior art date
Application number
PCT/JP2023/018114
Other languages
English (en)
French (fr)
Inventor
幹生 潮田
伸夫 原
Original Assignee
パナソニックＩｐマネジメント株式会社
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by パナソニックＩｐマネジメント株式会社 filed Critical パナソニックＩｐマネジメント株式会社
Publication of WO2023238606A1 publication Critical patent/WO2023238606A1/ja

Links

Images

Classifications

    • G—PHYSICS
    • G05—CONTROLLING; REGULATING
    • G05B—CONTROL OR REGULATING SYSTEMS IN GENERAL; FUNCTIONAL ELEMENTS OF SUCH SYSTEMS; MONITORING OR TESTING ARRANGEMENTS FOR SUCH SYSTEMS OR ELEMENTS
    • G05B13/00—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion
    • G05B13/02—Adaptive control systems, i.e. systems automatically adjusting themselves to have a performance which is optimum according to some preassigned criterion electric
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F30/00—Computer-aided design [CAD]
    • G06F30/20—Design optimisation, verification or simulation
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N20/00—Machine learning
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06N—COMPUTING ARRANGEMENTS BASED ON SPECIFIC COMPUTATIONAL MODELS
    • G06N7/00—Computing arrangements based on specific mathematical models
    • G06N7/01—Probabilistic graphical models, e.g. probabilistic networks
    • G—PHYSICS
    • G06—COMPUTING; CALCULATING OR COUNTING
    • G06Q—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES; SYSTEMS OR METHODS SPECIALLY ADAPTED FOR ADMINISTRATIVE, COMMERCIAL, FINANCIAL, MANAGERIAL OR SUPERVISORY PURPOSES, NOT OTHERWISE PROVIDED FOR
    • G06Q50/00—Information and communication technology [ICT] specially adapted for implementation of business processes of specific business sectors, e.g. utilities or tourism
    • G06Q50/04—Manufacturing
    • Y—GENERAL TAGGING OF NEW TECHNOLOGICAL DEVELOPMENTS; GENERAL TAGGING OF CROSS-SECTIONAL TECHNOLOGIES SPANNING OVER SEVERAL SECTIONS OF THE IPC; TECHNICAL SUBJECTS COVERED BY FORMER USPC CROSS-REFERENCE ART COLLECTIONS [XRACs] AND DIGESTS
    • Y02—TECHNOLOGIES OR APPLICATIONS FOR MITIGATION OR ADAPTATION AGAINST CLIMATE CHANGE
    • Y02P—CLIMATE CHANGE MITIGATION TECHNOLOGIES IN THE PRODUCTION OR PROCESSING OF GOODS
    • Y02P90/00—Enabling technologies with a potential contribution to greenhouse gas [GHG] emissions mitigation
    • Y02P90/30—Computing systems specially adapted for manufacturing

Definitions

  • the present disclosure relates to a technology for efficiently controlling product characteristic values within specifications in mass production of general industrial products.
  • the optimal solution for process conditions can often be searched for using mathematical optimization techniques when the relationship between process conditions and product characteristics can be expressed by a physical formula. However, if the relationship is unknown, a set of combinations of process condition values (i.e., control points) is selected to perform actual control (i.e., production). Then, as a control result, a combination of product characteristic values (that is, a characteristic point) corresponding to the control point is obtained. By repeating such control, it is possible to search for the optimal solution for the process conditions.
  • Patent Document 1 discloses a method of systematically adjusting a gain value using PID control, which is one of modern control methods.
  • Bayesian optimization is an optimization method that assumes a Gaussian process as a mathematical model that expresses the correspondence between input and output.
  • Bayesian optimization in a system where one set of input-output relationships is obtained each time control is performed, each time a control result is obtained, it is calculated based on the joint distribution of correspondence relationships that can be calculated from known control results. , calculate a predicted distribution marginalized by known control results.
  • an evaluation standard called an acquisition function is used to select the next optimal control condition. This enables quantitative evaluation that does not depend on the ability of the analyst, and can also contribute to the automation of optimal solution search work.
  • FIG. 13 is a diagram for explaining processing by the evaluation value calculation unit according to the embodiment.
  • FIG. 14 is a diagram illustrating an example of predicted distribution data according to the embodiment.
  • FIG. 15A is a diagram illustrating an example of an improvement area according to the embodiment.
  • FIG. 15B is a diagram showing another example of the improvement area according to the embodiment.
  • FIG. 16A is a diagram for explaining a method of calculating the volume of an improved region according to the embodiment.
  • FIG. 16B is a diagram illustrating an example of dividing the entire characteristic space into a plurality of small regions according to the embodiment.
  • FIG. 16C is a diagram illustrating an example of a lower end point and an upper end point of a small region according to the embodiment.
  • FIG. 17 is a diagram illustrating an example of a characteristic space divided into regions when the region reduction rule according to the embodiment is applied.
  • FIG. 18A is a diagram illustrating another example of a characteristic space divided into regions when the region reduction rule according to the embodiment is applied.
  • FIG. 18B is a diagram illustrating another example of the characteristic space divided into regions when the region reduction rule according to the embodiment is applied.
  • FIG. 18C is a diagram illustrating another example of the characteristic space divided into regions when the region reduction rule according to the embodiment is applied.
  • FIG. 19 is a diagram for explaining a Pareto boundary calculation method to which the area reduction rule according to the embodiment is applied.
  • FIG. 20 is a diagram for explaining a Pareto boundary calculation method to which the area reduction rule according to the embodiment is applied.
  • FIG. 21 is a diagram for explaining a Pareto boundary calculation method to which the area reduction rule according to the embodiment is applied.
  • FIG. 22 is a diagram for explaining a Pareto boundary calculation method to which the area reduction rule according to the embodiment is applied.
  • FIG. 23 is a diagram illustrating an example of a Pareto boundary when there are no constraint conditions according to the embodiment.
  • FIG. 24 is a diagram showing an example of an improvement area under the Pareto boundary shown in FIG. 23.
  • FIG. 25 is a diagram illustrating an example of a Pareto boundary when there are constraint conditions according to the embodiment.
  • FIG. 26 is a diagram showing an example of an improvement area defined when there are constraint conditions under the Pareto boundary shown in FIG. 25.
  • FIG. 32B is a diagram for explaining a control image of control method 1 in which the optimization target value is set at the center of the standard according to the embodiment.
  • FIG. 33 is a diagram conceptually showing an overshoot phenomenon and a hunting phenomenon.
  • FIG. 34 is a diagram conceptually showing the noise canceling effect, which is a characteristic of the Kalman filter.
  • FIG. 35A is a diagram for explaining a control image of control method 2 in which the optimization target value is set to the upper and lower limits of the standard according to the embodiment.
  • FIG. 35B is a diagram for explaining a control image of control method 2 in which the optimization target value is set to the upper and lower limits of the standard according to the embodiment.
  • FIG. 35A is a diagram for explaining a control image of control method 2 in which the optimization target value is set to the upper and lower limits of the standard according to the embodiment.
  • FIG. 35B is a diagram for explaining a control image of control method 2 in which the optimization target value is set to the upper and lower limits of the standard according
  • FIG. 36A is a diagram for explaining a control image of control method 3 in which the optimization target value according to the embodiment is gradually moved to the center of the standard.
  • FIG. 36B is a diagram for explaining a control image of control method 3 in which the optimization target value according to the embodiment is gradually moved to the center of the standard.
  • FIG. 37 is a diagram for explaining that there are multiple optimization target value candidates when the number of product characteristics is two or more dimensions according to the embodiment.
  • FIG. 38A is a diagram illustrating a combined active area when the center does not coincide with the standard center according to an embodiment.
  • FIG. 38B is a diagram illustrating a combined active area when the center coincides with the standard center according to the embodiment.
  • standard ranges may be set as constraints on the values of product characteristics. For example, ⁇ I want the battery capacity to fall within the standard range of 1850 to 1950 [mAh]'' or ⁇ The minimum value of the lifespan is 3 years, the longer the better (in other words, the minimum value of the standard range of lifespan is 3 years, and the maximum value is within the standard range such as ⁇ + ⁇ )''.
  • ⁇ I want the battery capacity to fall within the standard range of 1850 to 1950 [mAh]'' or ⁇ The minimum value of the lifespan is 3 years, the longer the better (in other words, the minimum value of the standard range of lifespan is 3 years, and the maximum value is within the standard range such as ⁇ + ⁇ )''.
  • Non-Patent Document 2 discloses EHVIC (Expected Hypervolume Improvement with Constraints), which is an extension of EHVI when there are constraint conditions.
  • the acquisition function design method described in Non-Patent Document 2 comprehensively indexes the probability of falling within the standard range and the amount of improvement, and evaluates all candidate control points. Therefore, there is a high possibility that search efficiency will improve (that is, a true optimal solution will be found).
  • Patent Document 1 discloses a method of systematically adjusting the gain value using PID control, which is one of the modern control methods.
  • the method disclosed in Patent Document 1 is a rule-based gain adjustment method, quantitative evaluation is difficult, and good control results that can suppress overshoot or hunting phenomena in controlling time-series data at mass production sites cannot be achieved. There is no guarantee that it will be possible. For this reason, in controlling time-series data at mass production sites, parameter evaluation becomes dependent on humans.
  • the Kalman filter is widely known as a calculation method for estimating states that change over time.
  • the Kalman filter is a calculation method that estimates invisible states inside a system using a mathematical model called a state space model. Therefore, the Kalman filter can be applied to control time-series data at mass production sites based on estimated state information.
  • an objective of the evaluation device of the present disclosure is to be able to quantitatively evaluate candidate control points that can suppress the occurrence of an overshoot phenomenon or a hunting phenomenon for a time-series data control problem.
  • An evaluation device includes a plurality of characteristic points corresponding to a plurality of candidate control points at a second time following the first time, based on known characteristic points corresponding to a controlled control point at a first time.
  • an evaluation device for evaluating unknown characteristic points of by Bayesian optimization the first receiving means acquiring control result data indicating controlled control points at the first time and known characteristic points at the first time; and the unknown characteristic point indicates a value of one or more product characteristics, at least one product characteristic has an optimization purpose, and a second receiving means for acquiring purpose data indicating the optimization purpose; a third receiving means for acquiring constraint data indicating a constraint given to the at least one product characteristic; and a method for dividing a characteristic space expressed by the at least one product characteristic; fourth receiving means for acquiring area reduction rule data indicating a dimension to be reduced for each area of the characteristic space divided by the dividing method; the control result data, the objective data, the constraint data, and the area reduction rule; Calculation means for calculating evaluation values for each of
  • the optimal candidate control point to be set next can be quantitatively calculated, and highly accurate control can be expected to be achieved regardless of the ability of the analyst.
  • the constraint condition may be at least one constraint range.
  • the optimization objective may include a first objective of keeping the product characteristics within one of the at least one constraint range, and a second objective of minimizing or maximizing the product characteristics.
  • the calculation means may calculate the evaluation value by performing different weighting processes for each of at least one product characteristic in the following cases (i) to (iii).
  • the plurality of constraint ranges are a standard range and a management range included in the standard range.
  • the first case is that the section is within the standard range and outside the control range
  • the optimization objective is the first objective
  • the first case is that the section is within the control range
  • the cases are divided into a second case where the optimization objective is the first objective, and a second case where the optimization objective is the first objective. Further, for example, weighting processing using a larger weight is performed in the second case than in the first case.
  • there are a plurality of constraint ranges and by further dividing the case (ii) into a plurality of cases, it is possible to weight each of the plurality of constraint ranges in stages. Therefore, even if the value of the product characteristic falls within the standard range and is desired to fall within the control range as much as possible, the evaluation value can be appropriately calculated. As a result, the scope of application to optimization problems can be further expanded.
  • the method may further include candidate control point creation means for creating the plurality of candidate control points by combining values that satisfy predetermined conditions of each of the plurality of process conditions.
  • the predetermined condition is that the sum of the values of each ratio variable of a plurality of process conditions is 1.
  • the ratio variable is a blending ratio of materials, such as compounds, corresponding to process conditions. Therefore, for each combination of compounding ratios of multiple types of compounds, an evaluation value for that combination can be calculated. As a result, it is possible to appropriately search for an optimal solution for at least one product property of the synthetic material obtained by blending these compounds.
  • the calculation means may calculate a predicted distribution at the plurality of candidate control points using a Kalman filter, and use the calculated predicted distribution to calculate the evaluation value.
  • a candidate control point to be set next as a control point can be selected based on the evaluation value.
  • the evaluation value calculated for each candidate control point is output, the user of the evaluation device can select the candidate control point as the next control point based on those evaluation values, and use that control point.
  • the characteristic points obtained by the controlled control can be used to calculate the evaluation value of each candidate control point.
  • the characteristic space is divided into an area within the standard range and an area outside the standard range, depending on the set standard range.
  • a constraint condition that is a standard range may be given to a product characteristic that has an optimization purpose.
  • the constraint condition is a condition given to a product characteristic, and includes, for example, a constraint range that specifies a range of values of a product characteristic as a condition. Examples of the constraint range include a standard range defined by product characteristic standards and a management range that can be set as appropriate by the user.
  • the evaluation device 100 describes the correspondence between candidate control points and characteristic points using a Kalman filter.
  • the Kalman filter used in the evaluation device 100 will be described below.
  • the mathematical model of the Kalman filter is to express the temporally changing internal state of the system using (Equation 1), where X (t) is the internal state of the system at time t, and Y (t) is the observed quantity at time t. can be expressed.
  • A, B, and C represent matrices that define the transformation.
  • v (t) and w (t) represent Gaussian noise at time t.
  • the average and variance of Gaussian noise can be set as appropriate by the analyst, such as being set to 0 and 1, for example.
  • Equation 1 is also called a state equation, and describes the temporal evolution of the internal state of the system that is not observed.
  • the lower equation of (Equation 1) is also called the observation equation, and describes the conversion from the internal state of the system to the observable quantity that we observe. If all internal state quantities of the system are observed, C may be used as a unit matrix. Since the values of each element of A, B, and C are usually unknown, it is also possible to proceed while estimating them sequentially using a time series analysis method such as an autoregressive model.
  • the Kalman filter formulated based on the above (Equation 1) is called a linear Gaussian Kalman filter.
  • the internal state amount at the next point in time predicted from the observed state amount is derived using a prediction distribution as shown in (Equation 2).
  • the multidimensional normal distribution has the property of preserving normality even if conditioning operations are performed on some dimensions.
  • the predicted distribution of the target characteristic Y as a product characteristic given the control factor X as a process condition can be derived as a normal distribution.
  • the mean and variance of the normal distribution are given by (Equation 5).
  • FIG. 3A shows the elements of matrix A defining the transformation and its submatrix AY .
  • FIG. 3B shows the elements of the matrix P- and its submatrices P - XX , P - XY , P - YY , P - YX .
  • FIG. 4 is a diagram showing the configuration of evaluation device 100 according to this embodiment.
  • the evaluation device 100 includes an input section 101a, a communication section 101b, an arithmetic circuit 102, a memory 103, a display section 104, and a storage section 105.
  • the input unit 101a is an HMI (Human Machine Interface) that accepts input operations by the user.
  • the input unit 101a is, for example, a keyboard, a mouse, a touch sensor, a touch pad, or the like.
  • the input unit 101a accepts setting information 210 as input from the user.
  • the setting information 210 includes process condition data 211, objective data 212, constraint data 213, and area reduction rule data 214.
  • the process condition data 211 is, for example, data indicating possible values of the process conditions, as shown in FIG. 2(a).
  • the values of the process conditions may be continuous values or discrete values.
  • the objective data 212 is, for example, data indicating an optimization objective of product characteristics such as minimization/maximization.
  • the constraint condition data 213 is, for example, data indicating a constraint condition such as a constraint range.
  • the region reduction rule data 214 is data indicating a rule for calculating Pareto boundaries, and changes the method for calculating the amount of improvement.
  • the region reduction rule data 214 indicates a method of dividing a characteristic space expressed by at least two product characteristics, and sets the dimension in which the active region is to be reduced to the region of the characteristic space divided by the dividing method. Shown below. Details will be described later.
  • the communication unit 101b connects to other devices by wire or wirelessly, and transmits and receives data to and from the other devices.
  • the communication unit 101b receives characteristic point data 201 from another device (for example, a control device).
  • the display unit 104 displays images, characters, etc.
  • the display unit 104 is, for example, a liquid crystal display, a plasma display, an organic EL (Electro-Luminescence) display, or the like.
  • the display section 104 may be a touch panel integrated with the input section 101a.
  • the storage unit 105 stores a program (ie, a computer program) 200 in which instructions to the arithmetic circuit 102 are written and various data.
  • the storage unit 105 is a nonvolatile recording medium, such as a magnetic storage device such as a hard disk, a semiconductor memory such as an SSD (Solid State Drive), or an optical disk.
  • the program 200 and various data may be provided to the evaluation device 100 from the other devices mentioned above via the communication unit 101b and stored in the storage unit 105, for example.
  • the storage unit 105 stores candidate control point data 221, control result data 222, predicted distribution data 223, and evaluation value data 224 as various data.
  • the candidate control point data 221 is data indicating each candidate control point.
  • each candidate control point is expressed by a combination of values of the first process condition and the second process condition.
  • the candidate control point data 221 may be data in a table format in which combinations of values of the first process condition and the second process condition are listed. A specific example of such candidate control point data 221 will be described in detail using FIG. 11A and FIG. 11B.
  • the control result data 222 is data indicating one or more control points used for control and characteristic points corresponding to each of the one or more control points.
  • the control result data 222 is a combination of a control point on the control space shown in FIG. 2(a) and a characteristic point on the characteristic space shown in FIG. 2(b) obtained by control using the control point. shows.
  • the control point is expressed by a combination of values of the first process condition and the second process condition.
  • a characteristic point is expressed by a combination of values of a first product characteristic and a second product characteristic.
  • the control result data 222 may be data in a table format in which combinations of control points and characteristic points are listed. A specific example of this control result data 222 will be explained in detail using FIG. 12.
  • the predicted distribution data 223 is data indicating the predicted distribution of all candidate control points indicated by the candidate control point data 221.
  • the predicted distribution is a distribution determined by (Equation 3) based on the Kalman filter as described above, and is expressed by, for example, the mean and variance.
  • the predicted distribution data 223 may be data in a table format that shows the predicted distribution of the first product characteristic and the predicted distribution of the second product characteristic in association with each other for each candidate control point. A specific example of this predicted distribution data 223 will be described in detail using FIG. 14.
  • the evaluation value data 224 is data indicating evaluation values for each of a plurality of candidate control points, as shown in FIG. 1, for example.
  • the evaluation value data 224 may be data in a table format that shows evaluation values associated with each of a plurality of candidate control points. Other specific examples of this evaluation value data 224 will be described in detail using FIG. 28 and the like.
  • FIG. 5 is a block diagram showing the functional configuration of the arithmetic circuit 102.
  • the reception control unit 10 receives characteristic point data 201, process condition data 211, objective data 212, constraint data 213, and area reduction rule data 214 via the input unit 101a or the communication unit 101b. For example, when the characteristic point data 201 is input by the user's input operation to the input unit 101a, the reception control unit 10 controls the storage unit 105 by associating the characteristic points indicated in the characteristic point data 201 with control points. Write to result data 222. As a result, the control result data 222 is updated. When this control result data 222 is updated, the reception control section 10 causes the evaluation value calculation section 12 to execute a process using the updated control result data 222. That is, the reception control unit 10 causes the evaluation value calculation unit 12 to calculate the evaluation value.
  • the reception control unit 10 may cause the evaluation value calculation unit 12 to start calculating the evaluation value in response to another trigger. For example, if the control result data 222 is already stored in the storage unit 105, the reception control unit 10, triggered by the user's input of the control point level, causes the evaluation value calculation unit 12 to start calculating the evaluation value. You may let them.
  • the level of the control point is, for example, the minimum value, maximum value, and discrete width of the values that the process conditions can take. That is, when the level of control points is input by the user and the candidate control point data 221 is generated based on the level, the reception control unit 10 generates the candidate control point data 221, the control result data 222, and the area reduction rule data.
  • the evaluation value calculation unit 12 is caused to start calculating an evaluation value based on 214.
  • the reception control unit 10 triggered by the user's input of the control result data 222, instructs the evaluation value calculation unit 12 to calculate the evaluation value. You may start it.
  • the reception control unit 10 instructs the evaluation value calculation unit 12 to calculate an evaluation value based on the control result data 222, candidate control point data 221, and area reduction rule data 214. and start it.
  • the reception control section 10 if it has the candidate control point data 221 and the control result data 222, it causes the evaluation value calculation section 12 to start calculating an evaluation value based on them.
  • the reception control unit 10 triggered by the user's input of the candidate control point data 221, instructs the evaluation value calculation unit 12 to calculate the evaluation value. You may start it.
  • the reception control unit 10 triggered by the input of the start instruction by the user, sends the evaluation value calculation unit 12 to the evaluation value calculation unit 12. You may also start the calculation of .
  • the evaluation value calculation unit 12 reads candidate control point data 221 and control result data 222 from the storage unit 105, generates predicted distribution data 223 based on these data, and stores the predicted distribution data 223 in the storage unit 105. . Further, the evaluation value calculation unit 12 generates evaluation value data 224 based on the predicted distribution data 223, the objective data 212, the constraint data 213, and the area reduction rule data 214 acquired by the reception control unit 10. , and stores the evaluation value data 224 in the storage unit 105.
  • FIG. 6 is a diagram showing an example of a reception image displayed on the display unit 104 to receive input of the setting information 210.
  • process condition data 211 corresponding to the input results is input to the evaluation apparatus 100.
  • FIG. 7 is a diagram showing an example of a reception image displayed on the display unit 104 to receive input of the area reduction rule data 214.
  • the evaluation device 100 when a radio button indicating "apply” is selected by the user's input operation on the input unit 101a, the evaluation device 100 applies the area reduction rule in the evaluation value calculation process. Furthermore, if the radio button indicating "apply” is selected by the user's input operation on the input unit 101a, the evaluation device 100 performs region division including the empty set in the evaluation value calculation process. Compute the volume of the improved region by applying the region reduction rule.
  • FIG. 8 is a diagram showing an example of the area reduction rule data 214.
  • the area reduction rule data 214 input by the area reduction rule area 330 of the reception image 300 in FIG. 7 includes changes in the definition of the Pareto boundary and in the definition of the active area, as shown in FIG. 8, for example. , and that region segmentation is applied. This changes the method of calculating the amount of improvement. Note that a specific example of the changed definition and applied area division will be described later, so a description thereof will be omitted here.
  • the discrete variable does not have a size relationship or numerical magnitude, such as "apple, tangerine, banana" or "with catalyst, without catalyst.”
  • the ratio variable indicates, for example, the blending ratio of the materials under the first process condition or the second process condition in a synthetic material produced by combining the materials under the first process condition and the material under the second process condition. .
  • the objective data 212 indicates, for example, an optimization objective for the first product characteristic and an optimization objective for the second product characteristic.
  • the constraint data 213 input through the product characteristic area 320 of the reception image 300 in FIG. 6 may indicate the standard range of the first product characteristic and the standard range of the second product characteristic.
  • the objective data 212 may indicate "within standard range” as the optimization objective for the first product characteristic, and may indicate "minimize” as the optimization objective for the second product characteristic.
  • FIG. 9A is a diagram illustrating an example of a standard range.
  • the standard range indicated by the constraint data 213 is expressed as a rectangular range on the characteristic space, as shown in FIG. 9A, for example.
  • the shape of the standard range is a rectangle, but it may have another shape.
  • the shape of the standard range may be any shape as long as it is possible to implement the calculation of the evaluation value described below.
  • FIG. 9B is a diagram showing another example of the standard range.
  • the standard range may be circular, for example, as shown in FIG. 9B.
  • the standard range in the characteristic space of the first product characteristic and the second product characteristic is expressed by the center (20, 20) and radius 10 of a circle.
  • the shape of the standard range may be a shape other than a circle, such as an ellipse or a star shape.
  • the evaluation value calculation unit 12 may calculate the evaluation value of each candidate control point based on the standard range of a shape different from a rectangle.
  • evaluation values are calculated based on standard ranges such as circular, oval, and star-shaped in the characteristic space, so the scope of application is further expanded without being limited to cases where the shape of the standard range is rectangular. be able to.
  • FIG. 10 is a flowchart showing the processing operation of the evaluation device 100 according to the present embodiment.
  • the candidate control point generation unit 11 generates candidate control point data 221 using the process condition data 211 (step S21).
  • the reception control unit 10 acquires the target data 212 (step S22). That is, the reception control unit 10 executes a second reception step of acquiring purpose data 212 indicating the optimization purpose.
  • the unknown characteristic point indicates the value of one or more product characteristics, and at least one product characteristic has an optimization objective.
  • the reception control unit 10 acquires constraint data 213 (step S23). That is, the reception control unit 10 executes a third reception step of acquiring constraint data 213 indicating a constraint given to at least one product characteristic. Further, the reception control unit 10 obtains area reduction rule data 214 indicating rules for calculating Pareto boundaries (step S24).
  • the reception control unit 10 specifies a region that indicates a division method of a characteristic space expressed by at least one product characteristic, and indicates a dimension in which the active region is to be reduced for each region of the characteristic space divided by the division method.
  • a fourth receiving step of acquiring reduction rule data 214 is executed. Further, the reception control unit 10 reads the control result data 222 from the storage unit 105 (step S25). That is, the reception control unit 10 executes a first reception step of acquiring control result data 222 indicating the controlled control point at the first time and the known characteristic point at the first time. Note that if the control result data 224 does not indicate any characteristic points, the processes of steps S25 to S27 including this step S25 are skipped.
  • the evaluation value calculation unit 12 calculates the evaluation value of each candidate control point based on the objective data 212, constraint data 213, area reduction rule data 214, candidate control point data 221, and control result data 222 (Ste S26). That is, the evaluation value calculation unit 12 executes a calculation step of calculating evaluation values for each of the plurality of unknown characteristic points based on the data. Specifically, the evaluation value calculation unit 12 calculates the evaluation value of each candidate control point shown in the candidate control point data 221 using an improvement amount calculation method changed based on the area reduction rule. Furthermore, in this calculation step, the evaluation value calculation unit 12 may assign weighting to the evaluation value for at least one product characteristic according to the degree of compliance with the constraint conditions. Then, the evaluation value calculation unit 12 generates evaluation value data 224 indicating the calculated evaluation value of each candidate control point.
  • the evaluation value output unit 13 outputs the evaluation value calculated in step S5, that is, the evaluation value data 224, to the display unit 104 (step S27). That is, the evaluation value output unit 13 executes an output step of outputting the evaluation value. Thereby, the evaluation value data 224 is displayed on the display unit 104, for example.
  • the reception control unit 10 acquires an operation signal from the input unit 101a in response to the user's input operation to the input unit 101a.
  • the operation signal indicates the end of the search for the optimal solution or the continuation of the search for the optimal solution.
  • the search for the optimal solution is a process of calculating and outputting the evaluation value of each candidate control point based on the new control result.
  • the reception control unit 10 determines whether the operation signal indicates the end of the search for the optimal solution or a continuation (step S28).
  • the reception control unit 10 determines that the operation signal indicates the end of the search for the optimal solution ("end" in step S28), it ends all processing.
  • the reception control unit 10 determines that the operation signal indicates the continuation of the search for the optimal solution (“Continue” in step S28)
  • the reception control unit 10 stores the candidate control point selected as the next control point as the control result in the storage unit 105. Write to data 222.
  • the reception control unit 10 selects a candidate control point from the evaluation value data 224 as the next control point.
  • the reception control unit 10 writes the candidate control points selected in this way into the control result data 222.
  • FIG. 11A is a diagram showing an example of candidate control point data 221.
  • the candidate control point generation unit 11 generates candidate control point data 221 shown in FIG. 11A based on the process condition data 211. For example, if each value of all the process conditions indicated by the process condition data 211 is the value of a continuous variable and there is no constraint regarding the value, the candidate control point creation unit 11 determines the value of each process condition. Each of all combinations is created as a candidate control point.
  • the process condition data 211 has the continuous variable values "10, 20, 30, 40, 50" of the first process condition and the continuous variable values "100, 200, 300, 400, 500" of the second process condition. Let us show that.
  • the candidate control point creation unit 11 combines the value "10" of the first process condition and the value "100” of the second process condition, the value "10” of the first process condition and the value of the second process condition. All combinations, such as the combination with "200", are created as candidate control points.
  • the candidate control point creation unit 11 associates a control point number with the created candidate control point, and generates candidate control point data 221 indicating the candidate control point with which the control point number is associated.
  • the candidate control point data 221 includes candidate control points (10, 100) associated with control point number "1" and candidates associated with control point number "2", as shown in FIG. 11A.
  • a control point (10,200), a candidate control point (10,300) associated with control point number "3", etc. are shown. Note that the first component of these candidate control points indicates the value of the first process condition, and the second component indicates the value of the second process condition.
  • FIG. 11B is a diagram showing another example of candidate control point data 221.
  • the candidate control point creation unit 11 generates candidate control point data 221 shown in FIG. 11B based on the process condition data 211.
  • the process condition data 211 indicates "0.0, 0.2, 0.4, 0.6, 0.8, 1.0" as the value of the ratio variable for the second process condition
  • the third process condition Assume that the values of the ratio variables are "0.0, 0.2, 0.4, 0.6, 0.8, 1.0".
  • the combination of values of these ratio variables corresponds to the above-mentioned blending ratio of the first compound and the second compound.
  • the candidate control point creation unit 11 adjusts the value of the first process condition and the second process condition so that the sum of the value of the ratio variable of the second process condition and the value of the ratio variable of the third process condition satisfies 1.
  • a combination of the process condition value and the third process condition value may be created as a candidate control point.
  • the candidate control point creation unit 11 creates a combination of the value "10" of the first process condition, the value "0.2" of the second process condition, and the value "0.8" of the third process condition, etc. Combinations of values such that the sum of ratio variable values satisfies 1 may be created as candidate control points.
  • the candidate control point creation unit 11 associates a control point number with the created candidate control point, and generates candidate control point data 221 indicating the candidate control point with which the control point number is associated.
  • candidate control point data 221 includes candidate control points (10, 0.0, 1.0) associated with control point number “1" and control point number "2". ”, the candidate control point (10, 0.2, 0.8) associated with the control point number “3”, the candidate control point (10, 0.4, 0.6) associated with the control point number “3”, etc.
  • the first component of these candidate control points indicates the value of the first process condition
  • the second component indicates the value of the second process condition
  • the third component indicates the value of the third process condition.
  • the candidate control point creation unit 11 creates each of the plurality of candidate control points using a predetermined condition for each of the plurality of process conditions.
  • the candidate control point is created by combining values that satisfy the following.
  • the predetermined condition is that the sum of the values of the ratio variables of the plurality of process conditions is 1, as shown in FIG. 11B.
  • the ratio variable is a blending ratio of materials, such as compounds, corresponding to process conditions. Therefore, for each combination of compounding ratios of multiple types of compounds, an evaluation value for that combination can be calculated. As a result, it is possible to appropriately search for optimal solutions for one or more product properties of the synthetic material obtained by blending these compounds.
  • FIG. 12 is a diagram showing an example of the control result data 222.
  • the evaluation value calculation unit 12 reads the control result data 222 stored in the storage unit 105 in order to calculate the evaluation value.
  • this control result data 222 indicates, for each control number, the control point used in the control identified by that control number and the characteristic point that is the control result obtained by that control. .
  • a control point is expressed by a combination of values of each process condition.
  • the control point is expressed by a combination of values that is a combination of the value "10" of the first process condition and the value "100" of the second process condition.
  • a characteristic point is expressed by a combination of values of each product characteristic obtained through control.
  • the value of the product characteristic is hereinafter also referred to as a product characteristic value.
  • the characteristic point is expressed by a combination of the value "8" of the first product characteristic and the value "0.0" of the second product characteristic.
  • control result data 222 includes control points (10, 100) and characteristic points (8, 0.0) associated with control number "1" and control number "2". ” control point (10,500) and characteristic point (40, 1.6), control point (50,100) and characteristic point (40, 1.6) associated with control number “3”, etc. shows.
  • FIG. 13 is a diagram for explaining processing by the evaluation value calculation unit 12.
  • the evaluation value calculation unit 12 generates predicted distribution data 223 based on the candidate control point data 221 generated by the candidate control point generation unit 11 and the control result data 222 stored in the storage unit 105.
  • the evaluation value calculation unit 12 then generates objective data 212 indicating the optimization purpose of each product characteristic, constraint condition data 213 indicating the standard range of each product characteristic, and area reduction rule data 214 indicating the rule for calculating the Pareto boundary.
  • Evaluation value data 224 is generated based on the predicted distribution data 223 and the estimated distribution data 223 .
  • control result data 222 corresponds to one or more control points that are one or more candidate control points that have already been used for control among the plurality of candidate control points, and each of the one or more control points.
  • the evaluation value calculation unit 12 describes the correspondence between candidate control points and characteristic points using the Kalman filter described above.
  • the evaluation value calculation unit 12 calculates an evaluation value based on an evaluation standard called an acquisition function in Bayesian optimization.
  • the above-mentioned predicted distribution is used to calculate this evaluation value.
  • the acquisition function in this embodiment is an acquisition function in Bayesian optimization with constraint conditions.
  • EHVI is defined as the expected value of the amount of improvement in the prediction distribution for each candidate control point, as shown in (Equation 6) below.
  • a candidate control point with a larger value obtained by EHVI has a larger expected value of improvement amount, and represents a control point to be executed next.
  • R minimize represents a region in which all product characteristics y 1 to y Dminimize , whose optimization objective is minimization, are within the standard range.
  • R range represents a region in which product characteristics y Dminimize+1 to y D whose optimization objective is within the standard range are all within the standard range.
  • each region of R minimize and R range is expressed by a function indicating the shape of the standard range corresponding to the region. As shown in FIG. 9B, if the shape of the standard range is a circle, each region of R minimize and R range is expressed by a function representing the circle. Further, if the shape of the standard range is star-shaped, each region of R minimize and R range is expressed by a function indicating the star shape.
  • y new,minimize represents a vector obtained by extracting each dimension of a product characteristic whose optimization objective is minimization from all dimensions of the characteristic point y new .
  • y new,range represents a vector obtained by extracting each dimension of a product characteristic whose optimization objective is within the standard range from all dimensions of the characteristic point y new .
  • IC(y new ) is the amount of improvement when there is a constraint condition, and represents the volume of the area surrounded by the existing Pareto boundary and the newly determined Pareto boundary.
  • the existing Pareto boundary is a boundary determined from at least one Pareto point existing within the specification range and the respective coordinates of the specification range.
  • the evaluation value calculation unit 12 calculates the amount of improvement (i.e., IC(y new )), which is the volume of the improvement area, as shown in FIG. 16A, for the dimension of the product characteristic whose optimization objective is minimization. That is, the evaluation value calculation unit 12 divides the improvement area into a plurality of small areas at the respective coordinates of the Pareto point and the new characteristic point, calculates the expected volume of each small area, and then sums the expected values. The amount of improvement (ie, IC(y new )) is calculated by taking . Furthermore, the evaluation value calculation unit 12 calculates the probability that each product characteristic value falls within the standard range for a dimension of the product characteristic whose optimization objective is within the standard range.
  • IC(y new ) the amount of improvement
  • FIG. 16B is a diagram showing an example of dividing the entire characteristic space into a plurality of small regions.
  • the priority may be the reciprocal of the weighting coefficient cd . Unless otherwise specified, that is, if each dimension d has the same priority, all c d of each dimension d is set to 1, for example.
  • the volume within the constraint range in the characteristic space is calculated as the optimization improvement amount, and the improvement is calculated.
  • the evaluation value can be appropriately calculated from the amount.
  • the improvement region is divided into multiple small regions, the expected volume of each small region is calculated, and It is necessary to calculate the sum of the expected values of . More generally, when the number of product characteristics, that is, the number of dimensions, is D, the improvement area can be expressed as a sum area of a plurality of D-dimensional hypercuboids. Therefore, when calculating the acquisition function for Bayesian optimization when there are constraints, the improvement region is divided into multiple D-dimensional hypercuboids, the expected value of the volume of each hypercuboid is calculated, and then those expected values are calculated.
  • the number of hypercuboids is determined by using D, which is the number of product characteristics (number of dimensions), and the number of Pareto points, N pareto , which is the number of Pareto points among the observed characteristic points.
  • the region reduction rules indicate a method for dividing the characteristic space into a predetermined number of regions and a method for calculating Pareto boundaries. In the following description, in order to simplify the explanation, it will be assumed that no standard range is set.
  • the entire characteristic space is divided into D+1 regions (region division) for D product characteristics, that is, D-dimensional product characteristics. Any method may be used for region division. It is assumed that the regions of the divided characteristic space are named region 1, region 2, . . . , region D, region D+1 in order.
  • FIG. 17 is a diagram showing an example of a characteristic space that is divided into regions when the region reduction rule is applied.
  • the characteristic space is divided into three regions, namely region 1, region 2, and region 3.
  • a third region consisting of the range of the first product characteristic from “- ⁇ " to "10” and the range of the second product characteristic from “- ⁇ ” to “10” is shown in the characteristic space.
  • the first region is an area below the straight line defined by a 45-degree inclination in the region excluding the third region, and the straight line defined by the 45-degree inclination in the region excluding the third region.
  • a second region which is also an upper region, is shown.
  • an empty set may be set in D+1 regions of the divided characteristic space, but an empty set for all D+1 regions is nonsense, so at least one Let be a non-empty set.
  • FIGS. 18A to 18C are diagrams showing other examples of the characteristic space divided into regions when the region reduction rule according to the present embodiment is applied.
  • FIG. 18A shows a case where the characteristic space is divided into two regions by setting region 3 to an empty set for two-dimensional product characteristics. More specifically, as shown in FIG. 18A, since the third region is an empty set, the first region, which is a region below a straight line defined by a 45-degree inclination in the characteristic space, and the characteristic The area is divided into a second area which is an area above a straight line defined by an inclination of 45 degrees in space.
  • the example shown in FIG. 18C shows a case where the characteristic space is divided into two regions by setting region 3 to an empty set for two-dimensional product characteristics. That is, in the example shown in FIG. 18C, since the third region is an empty set, the entire characteristic space is divided into the first region and the second region. More specifically, as shown in FIG. 18C, in the characteristic space, there is a first region that is a region represented by the center (10, 10) of a circle and a radius 5, and a region other than the first region in the characteristic space. The area is divided into a second area and a second area.
  • an empty set may be set in the D+1 regions of the region-divided characteristic space.
  • any point on the characteristic space is assigned to any one area other than the empty set.
  • the Pareto point is a characteristic point that is provisionally considered a Pareto solution at this point, and is also called a non-inferior solution.
  • the optimization objective of each of the first product characteristic and the second product characteristic is minimization.
  • a Pareto point is a characteristic point for which there is no other characteristic point for which the value of either the first product characteristic or the second product characteristic is smaller than that point compared to all other observed characteristic points. Become.
  • a Pareto boundary is a boundary determined from the coordinates of at least one Pareto point.
  • the above-mentioned Pareto boundary was a boundary line determined by extending and connecting the coordinates of Pareto points in the direction in which the values of the first product characteristic and the second product characteristic are large.
  • the Pareto boundary calculation method is changed. That is, when the area reduction rule is applied, the definition of the Pareto boundary is changed by determining the dimension by which the active area is reduced for each area.
  • FIGS. 19 to 22 are diagrams for explaining a method of calculating Pareto boundaries to which the area reduction rule is applied.
  • FIGS. 19 to 22(a) show examples of the Pareto boundary calculation method before the definition change, that is, the Pareto boundary calculation method to which the area reduction rule is not applied.
  • FIGS. 19 to 22(b) show examples of the Pareto boundary calculation method after the definition change, that is, the Pareto boundary calculation method to which the area reduction rule is applied.
  • the optimization objective of each of the first product characteristic and the second product characteristic is minimization in a characteristic space composed of two-dimensional product characteristics, for example.
  • FIGS. 19 to 22(b) it is assumed that the characteristic space is divided into region 1 and region 2 by a straight line passing through the origin and defined by an inclination of 45 degrees.
  • the new characteristic point y new(1) becomes a Pareto point.
  • the coordinates of the new characteristic point y new(1) are set in the direction of the first product characteristic and the second product characteristic.
  • the boundary line determined by connecting along is the Pareto boundary.
  • the new characteristic point y new(1) is located in area 2, so the new characteristic point y new(1 ) passing through the coordinates of the first product characteristic and parallel to the axis of the second product characteristic is the Pareto boundary.
  • the coordinate y new1 (1) of the first product characteristic of the new characteristic point y new (1) is larger than the coordinate y new2 (1) of the second product characteristic of the new characteristic point y new (1) . Therefore, the area to the right of the new characteristic point y new1(1) is set as an inactive area. This can also be expressed as reducing (reducing) the active area at the coordinates y new1(1) of the new characteristic point y new(1) .
  • a boundary line formed by a line parallel to the axis of the second product characteristic is a Pareto boundary.
  • the coordinate y new2(2) of the second product characteristic of the new characteristic point y new (2) is larger than the coordinate y new1(2) of the first product characteristic of the new characteristic point y new(2) . Therefore, the area to the right of the new characteristic point y new1 (1) or the area above the new characteristic point y new2 (2) is set as an inactive area. Furthermore, this can also be expressed as further reducing (reducing) the active area at the coordinates y new2(2) of the new characteristic point y new(2) .
  • the new characteristic point y new(3) when the third new characteristic point y new(3) is obtained, the new characteristic point y new(3) does not become a Pareto point. In this case, the Pareto boundary will not be changed, as shown in FIGS. 21(a) and (b). In other words, if the new characteristic point y new(3) is not included in the active area but is included in the inactive area, it does not become a Pareto point, so the Pareto boundary is not changed.
  • the new characteristic point y new(4) becomes a Pareto point.
  • the new characteristic point y new(4) is included in the active region, so the Pareto boundary is changed.
  • the new characteristic point y new(4) is included in the inactive area, so the Pareto boundary will not be changed.
  • the definition of the Pareto boundary is changed by determining the dimension by which the active area is reduced for each area.
  • FIG. 23 is a diagram showing an example of a Pareto boundary when there are no constraint conditions.
  • FIG. 23 shows an example of a Pareto boundary calculated when the area is divided as shown in FIG. 17.
  • the optimization objective of each of the first product characteristic and the second product characteristic is minimization in the characteristic space composed of, for example, two-dimensional product characteristics.
  • region 1 and region 2 are used for calculating the Pareto boundary, while region 3 is not used for calculating the Pareto boundary.
  • the initial value of y'd is the upper limit of each standard.
  • the area expressed by (Equation 9) becomes the Pareto boundary when the area reduction rule is applied.
  • D represents the number of product characteristics (number of dimensions)
  • FIG. 24 is a diagram showing an example of an improvement area under the Pareto boundary shown in FIG. 23.
  • the acquisition function (i.e. EHVIC with constraints) when the region reduction rule is applied uses a prediction distribution calculated based on a Kalman filter for each candidate control point, similar to EHVI. defined. More specifically, the acquisition function when the area reduction rule is applied can be defined by the expected value of the improvement amount, as shown in (Equation 7). Then, the quality of the next control point to be executed is evaluated based on the expected value of the amount of improvement.
  • the evaluation value calculation unit 12 calculates the characteristic point that has the smallest coordinate y'd among the characteristic points included in the area d and is within the standard range.
  • the boundary defined by the coordinate y' d is calculated as the active boundary.
  • the evaluation value calculation unit 12 can calculate the evaluation of the acquisition function using (Formula 7) and (Formula 8).
  • the evaluation value calculation unit 12 further sets an intermediate value between 0 and 1, such as 0.5, for l(yd, yd') in (Equation 8). By doing so, the evaluation value can be calculated. Note that 0.5 is just an example, and other numerical values may be used.
  • the user can determine whether to continue or end the search for the optimal solution. Furthermore, when continuing the search for the optimal solution, the user selects the next one from all displayed control point numbers, that is, all candidate control points, based on each displayed evaluation value and each ranking Candidate control points to be control points can be selected. For example, the user selects the candidate control point corresponding to the largest evaluation value (ie, the evaluation value with a rank of 1). At this time, the user may rearrange the evaluation values of the evaluation value data 224 in ascending order by performing an input operation on the input unit 101a. That is, the evaluation value output unit 13 sorts the evaluation value data 224 so that each evaluation value is in descending order and each rank is in ascending order. This makes it easier to find the largest evaluation value.
  • state estimation using a Kalman filter is conceptually based on the observed value at the previous time at each time (t-4, t-3, t-2, t-1, t). This corresponds to sequentially estimating the mean and variance of the predicted distribution of the state. An estimated value can be obtained based on the mean and variance of this predicted distribution. Note that sequentially estimating the mean and variance of the predicted distribution of the state can also be interpreted as being corrected using the observed value at the previous time.
  • the prediction distribution is calculated using the Kalman filter. Find (the mean and variance of). Then, the next control point to be set is selected based on the evaluation function determined by Bayesian optimization.
  • the optimization method is the same as control method 1, and if the position of the immediately preceding characteristic point (observed value) is above the center of the standard, the optimization method (optimization objective) is minimized; If the position of the characteristic point is below the center of the specification, control may be performed to switch the optimization objective to maximization.
  • FIGS. 35A and 35B are diagrams for explaining a control image of control method 2 in which the optimization target value according to the present embodiment is set to the upper and lower limits of the specification. Note that in FIGS. 35A and 35B as well, control images are shown where the number of product characteristics is one-dimensional to simplify the explanation.
  • time-series data control method 2 when the optimization objective is within the standard range and the previous characteristic point yt is above the center of the standard as in the example shown in FIG. 35A, the lower limit of the standard is set to the optimization target value. Set to control.
  • FIG. 35B when the immediately preceding characteristic point yt +1 is below the standard center, control is performed by setting the standard upper limit to the optimization target value.
  • Control may be performed by setting the local maximum value as the optimization target value. Note that the optimization method (optimization objective) is controlled as maximization because the immediately preceding characteristic point y t+1 is below the standard center.
  • EHVIC which is an acquisition function of Bayesian optimization
  • time series data control method 3 3.
  • control will be described as an example where the number of product characteristics is two or more dimensions.
  • EHVIC The idea of EHVIC is to divide the area in the characteristic space to be evaluated so that when maximizing or minimizing multiple product characteristics simultaneously, the search for the optimal solution proceeds efficiently based on the observed Pareto points. become For example, when there are two product characteristics, even if one product characteristic obtains a characteristic point extremely close to the optimization target value, that value does not mean that the entire characteristic space is partitioned, and all Pareto points are Partition in stages by coordinates. Utilizing this idea of EHVIC, in the time-series data control method 3, it may be possible to suppress the optimization target value from rapidly moving from the initial value, which is the upper and lower limits of the specification, to the center of the specification.
  • EHVIC has a problem in that the calculation cost increases exponentially depending on the number of product characteristics. Therefore, by using the area reduction rule, which is a rule for reducing the active area, the calculation cost can be reduced to the order of a polynomial function while maintaining search accuracy.
  • the optimization target value cannot be uniquely determined if the above-mentioned area reduction rule is applied as is. Therefore, in the following, a method will be described in which when the number of product characteristics is two or more dimensions, the optimization target value can be uniquely determined even when the above-mentioned area reduction rule is applied.
  • FIG. 37 is a diagram for explaining that there are multiple optimization target value candidates when the number of product characteristics is two or more dimensions.
  • the entire characteristic space can be divided into 2D divided areas by dividing each product characteristic into areas above and below the center of the standard. In each divided region, the region reduction rule is applied, and the Pareto boundary on the opposite side is set as the optimization target value for each product characteristic.
  • FIG. 37 shows four divided areas when the number of product characteristics is two-dimensional. Further, FIG. 37 shows Pareto boundaries and active regions calculated from observed characteristic points existing in each of the four divided regions. Note that the Pareto boundary shown in FIG. 37 is a temporary Pareto boundary whose optimization target value cannot be uniquely determined, so it will be referred to as a temporary Pareto boundary below. Further, the active area under the temporary Pareto boundary will be referred to as a temporary active area below. Furthermore, we will refer to the non-temporary Pareto boundary as a combined Pareto boundary, and the active area under the combined Pareto boundary as a combined active area.
  • the optimization method (optimization purpose) in the first product characteristic (Y 1 ) is determined to be minimized, and the optimization target value may be determined from the Pareto boundary that exists in the region below the center of the specification.
  • the optimization target value may be determined from the Pareto boundary that exists in the region below the center of the specification.
  • FIG. 37 for example, if there is a provisional Pareto boundary shown in two ways, A and B, in the first product characteristic (Y 1 ), there are two candidates for the optimization target value. It turns out.
  • FIG. 38A is a diagram showing a combined active area when the center does not coincide with the standard center.
  • FIG. 38B is a diagram showing the combined active area when the center coincides with the standard center.
  • the combined Pareto boundaries are individually set above and below the center of the specification for the first product characteristic (Y 1 ), and above and below the center of the specification for the second product characteristic (Y 2 ). It is better not to specify it. This is because although the combined active area under the individually determined combined Pareto boundary becomes a single rectangle as shown in FIG. 38A, the center is likely to be at a position different from the standard center. Additionally, when controlling time series data using control method 3 when the center of the combined active area is at a different position from the standard center, the control proceeds so that the observed values are driven toward the center of the combined active area. . For this reason, there is a high possibility that the control will be deviated from the standard center, which is the original target (optimization target value).
  • the combined Pareto boundaries for each product characteristic are uniquely defined in the upper region and the lower region from the center of the specification.
  • the first product characteristic (Y 1 ) there are two provisional Pareto boundaries that exist in the area above the center of the specification, and two provisional Pareto boundaries that exist in the area below the center of the specification, for a total of 4
  • the Y1 coordinate that is larger than the standard center by an intermediate value such as the average value of the distance from the standard center may be determined as the joint Pareto boundary of the upper region.
  • the combined active area under the determined combined Pareto boundary can be made into a single area as shown in FIG. 38B. can be made into a rectangle, and its center can be made to coincide with the standard center.
  • the method of determining the combined Pareto boundary in both the upper and lower regions from the center of the standard is not limited to using the intermediate position of the Y 1 coordinate of multiple temporary Pareto boundaries;
  • a joint Pareto boundary may be defined that passes through the positions where the values are closest or farthest. If you define a joint Pareto boundary that passes through the position where the value is closest (Y 1 coordinate), an overshoot phenomenon is likely to occur, and the joint Pareto boundary that passes through the position where the value is the farthest (Y 1 coordinate) is determined. If this is determined, the hunting phenomenon is likely to occur. Therefore, unless there is a particular reason, it is sufficient to define a combined Pareto boundary that passes through the position where the above-mentioned intermediate value is obtained.
  • time-series data control method 3 As described above, simply combining the Kalman filter and Bayesian optimization will cause an overshoot phenomenon. Therefore, by performing time-series data control method 3, the hunting phenomenon can also be suppressed by setting the initial value of the optimization target value to the upper and lower limits of the specification and gradually moving it toward the center of the specification.
  • time-series data control method 3 a Bayesian optimization acquisition function is used for quantitative evaluation of candidate control points, and a high-speed algorithm that suppresses calculation costs by applying area reduction rules is implemented. realizable.
  • the analyst (user) can inform the analyst (user) of the candidate control points to be set next in order of ranking, etc. offered and selected. Therefore, it is possible to quantitatively evaluate candidate control points that can suppress the occurrence of overshoot or hunting phenomena, regardless of the analyst, thereby suppressing the occurrence of overshoot or hunting phenomena and generating stable and efficient time series data. real-time control can be achieved.
  • the evaluation device selects a plurality of candidate control points at the second time following the first time based on the known characteristic points corresponding to the controlled control points at the first time.
  • An evaluation device that evaluates a plurality of corresponding unknown characteristic points by Bayesian optimization, wherein the evaluation device acquires control result data indicating a controlled control point at the first time and a known characteristic point at the first time.
  • the calculation means can calculate the evaluation values of the plurality of unknown characteristic points at the second time following the first time based on the control result data, objective data, constraint condition data, and area reduction rule data.
  • weighting is given to the evaluation value for at least one product characteristic according to the degree of compliance with the constraint conditions, and this at least one product characteristic has an optimization purpose. Therefore, the Kalman filter and Bayesian optimization can be applied to an optimization problem in which constraints are imposed on time-series product characteristics that have the objective of the optimization problem. In this way, evaluation values of a plurality of unknown characteristic points can be calculated from the known characteristic points corresponding to the controlled control point at the first time, which is past control result information. Control points can be quantitatively analyzed.
  • the constraint condition is the standard range
  • the optimization objectives include a first objective of keeping the product characteristics within the standard range and a second objective of minimizing or maximizing the product characteristics.
  • the evaluation value calculation unit 12 determines whether (i) the section of the product characteristic used to calculate the evaluation value is outside the standard range, and (ii) the section is outside the standard range. (iii) when the interval is within the specification range and the optimization objective is the first objective; and (iii) when the interval is within the specification range and the optimization objective is the second objective. Evaluation values are calculated by performing different weighting processes. That is, the evaluation value is calculated based on the above (Formula 6) and (Formula 7).
  • the evaluation value of the candidate control point can be appropriately calculated based on Bayesian optimization. That is, even if the purpose of optimizing product characteristics is within the standard range, maximizing or minimizing, the evaluation value of the candidate control point can be appropriately calculated based on Bayesian optimization.
  • the interval of product characteristics is within the standard range and the optimization objective is the second objective, so unlike the method of Non-Patent Document 2, the optimization objective is to maximize Alternatively, even if the product characteristic to be minimized has a standard range as a constraint, the evaluation value can be quantitatively and appropriately calculated.
  • the calculation means may calculate a predicted distribution at the plurality of candidate control points using a Kalman filter, and use the calculated predicted distribution to calculate the evaluation value.
  • a candidate control point to be set as the next control point can be selected based on the evaluation value.
  • the user of the evaluation device can select the candidate control point as the next control point based on those evaluation values, and use that control point.
  • the characteristic points obtained by the controlled control can be used to calculate the evaluation value of each candidate control point.
  • the evaluation value calculation unit 12 may calculate the evaluation value of each candidate control point using the Monte Carlo method.
  • the Monte Carlo method is an approximation method, even if it is difficult to calculate the evaluation value analytically, it can be calculated approximately. In other words, it is possible to calculate approximately by using the Monte Carlo method without strictly performing the calculations of (Formula 6) and (Formula 7). Note that as long as it is an approximation method, other methods may be used instead of the Monte Carlo method.
  • Example 2 As an example of the processing when searching for the optimal solution of product characteristics by applying the area reduction rule to the above control method 3, the process of determining the active area from the area reduction rule when a set of control results is added. An example will be explained.
  • FIG. 39 is a diagram showing an example of a control result data sheet obtained when searching for an optimal solution according to an example of this embodiment.
  • FIG. 40A is a diagram showing the combined active area at the time when the control results up to the 26th in FIG. 39 are obtained.
  • FIG. 40B is a diagram showing the combined active region at the time when the 27th control result in FIG. 39 is obtained. Note that the prediction distribution calculation process using the Kalman filter and the acquisition function calculation process using EHVI and EHVIC are the same as described above, and therefore their explanation will be omitted. Furthermore, non-standard characteristic points are omitted in FIGS. 40A and 40B.
  • FIG. 39 shows, for each control number, the process conditions of the control point used in the control identified by that control number, and the product characteristics of the characteristic point that is the control result obtained by that control. . Furthermore, FIG. 39 shows a case where the number of process conditions is two (first process condition, second process condition) and the number of product characteristics is two (first product characteristic, second product characteristic 2). has been done.
  • the optimization purpose is to minimize both the first product characteristic (Y 1 ) and the second product characteristic (Y 2 ), and the standard range of the first product characteristic (Y 1 ) is 1.1 to 1.7, and the second product characteristic is The standard range of (Y 2 ) is 0.5 to 0.7.
  • the initial value of the combined active area is the standard range, that is, the initial value of the combined Pareto boundary is the standard upper and lower limits of each dimension.
  • the standard range in this embodiment can be expressed as (Equation 11).
  • the regions (regions R 1-1 to R 4-2 ) to which the region reduction rules in FIGS. 40A and 40B are applied are determined as in (Equation 12).
  • the minimum value of the coordinates of the first product characteristic (Y 1 ) among the characteristic points in each region is set as a temporary Pareto boundary.
  • the maximum value of the coordinates of the first product characteristic (Y 1 ) among the characteristic points in each region is set as a temporary Pareto boundary.
  • the minimum value of the coordinates of the second product characteristic (Y 2 ) among the characteristic points in each region is set as a temporary Pareto boundary.
  • the maximum value of the coordinates of the second product characteristic (Y 2 ) among the characteristic points in each region is set as a temporary Pareto boundary.
  • a boundary passing through coordinates that are separated from the standard center by the average value of the distance from the standard center of the temporary Pareto boundary is defined as a combined Pareto boundary, and a rectangular area defined thereby is determined as a combined active area.
  • the coordinates in the first product characteristic (Y 1 ) will be referred to as Y 1 coordinates
  • the coordinates in the second product characteristic (Y 2 ) will be referred to as Y 2 coordinates.
  • the characteristic points A, B, C, D, E, F, and G shown in FIG. 40A are the region R 1-1 , region R 1-2 , region R 2-1 , region R 2-2 , and region shown in FIG. 40A. These are Pareto points in R 3-1 , region R 3-2 , and region R 4-1 . Note that in the example shown in FIG. 40A, there is no Pareto point in region R 4-2 . Furthermore, in the example shown in FIG. 40A, the joint Pareto boundary of the first product characteristic (Y 1 ) is between the Y 1 coordinates of characteristic points A, C, E, and G and the standard center of the first product characteristic (Y 1 ).
  • the position passes through the Y1 coordinate, which is away from the center of the standard by the average value of the distance.
  • the joint Pareto boundary of the second product characteristic (Y 2 ) is the Y 2 coordinates of characteristic points B, D, F, the standard upper limit value of the second product characteristic (Y 2 ), and the standard center of the second product characteristic (Y 2 ).
  • the position passes through the Y2 coordinate, which is away from the center of the standard by the average value of the distance between.
  • the upper and lower limits of the first product characteristic (Y 1 ) of the combined active region shown in FIG. 40A are as follows. That is, the upper limit of the first product characteristic (Y 1 ) is the average value (0.017975) of the distance between the Y 1 coordinates (1.4256, 1.4181, 1.3791, 1.3927) of characteristic points A, C, E, G and the center of the specification (1.4000). ) is added from the standard center to the Y 1 coordinate (1.417975).
  • the lower limit of the first product characteristic (Y 1 ) is the average value (0.017975) of the distance between the Y 1 coordinates (1.4256, 1.4181, 1.3791, 1.3927) of characteristic points A, C, E, G and the center of the specification (1.4000). ) is subtracted from the standard center to give the Y 1 coordinate (1.382025).
  • the upper and lower limits of the second product characteristic (Y 2 ) of the combined active region shown in FIG. 40A are as follows. In other words, the upper limit of the second product characteristic (Y 2 ) is the average distance between the Y 2 coordinates (0.6173, 0.5790, 0.5239) of characteristic points B, D, and F and the standard upper limit value (0.7) and the standard center (0.6000).
  • the value (0.053600) is added to the Y2 coordinate (0.653600) from the center of the standard.
  • the lower limit of the second product characteristic (Y 2 ) is the Y 2 coordinate of characteristic points B, D, and F (0.6173, 0.5790, 0.5239) and the average distance between the lower limit value of the specification (0.5) and the center of the specification (0.6000).
  • the Y2 coordinate (0.546400) is obtained by subtracting the value (0.053600) from the standard center.
  • the evaluation value calculation unit 12 predicts the 27th characteristic point for each candidate control point using a Kalman filter update formula such as (Formula 5) based on the control results with control numbers up to the 26th. Calculate the distribution.
  • the evaluation value calculation unit 12 calculates an acquisition function for each candidate control point using the calculated prediction distribution and a Bayesian optimization update formula such as (Formula 7).
  • the candidate control point with the maximum evaluation value obtained from the calculated acquisition function is adopted as the 27th control point.
  • the characteristic point with the 27th most recent control number is within the standard range and belongs to region R 4-2 .
  • the characteristic point with the 27th control number corresponds to characteristic point H in FIG. 40B, and becomes a new Pareto point belonging to region R 4-2 . Therefore, even if the active area is removed in region R 4-2 , the upper and lower limits of the first product characteristic (Y 1 ) of the combined active area do not change and are not updated, and the second product characteristic (Y 2 ) of the combined active area The upper and lower limits of are updated.
  • the upper limit of the second product characteristic (Y 2 ) is only the average value (0.030450) of the distance between the Y 2 coordinates (0.6173, 0.5790, 0.5239, 0.6074) of characteristic points B, D, F, and H and the center of the specification (0.6000). This is the Y2 coordinate (0.630450) added from the center of the standard.
  • the lower limit of the second product characteristic (Y 2 ) is only the average value (0.030450) of the distance between the Y 2 coordinates (0.6173, 0.5790, 0.5239, 0.6074) of characteristic points B, D, F, and H and the center of the specification (0.6000). This is the Y2 coordinate (0.569550) subtracted from the center of the standard.
  • characteristic point with the 27th control number corresponds to characteristic point H in FIG. Located higher up.
  • the optimization method (optimization objective) for the first product characteristic (Y 1 ) and the second product characteristic (Y 2 ) is set to maximization and minimization. Further, the optimization target value for the first product characteristic (Y 1 ) is set to the coordinates (1.417975) of the upper limit of the first product characteristic (Y 1 ) in the combined active region. The optimization target value in the second product characteristic (Y 2 ) is set to the coordinate (0.569550) of the lower limit of the second product characteristic (Y 2 ) in the combined active region.
  • each component may be configured with dedicated hardware, or may be realized by executing a software program suitable for each component.
  • Each component may be realized by a program execution unit such as a CPU or a processor reading and executing a software program recorded on a recording medium such as a hard disk or a semiconductor memory.
  • the software that implements the evaluation apparatus and the like of the above embodiments is a program that causes a computer to execute each step of the flowchart shown in FIG. 10, for example.
  • At least one of the above devices is specifically a computer system consisting of a microprocessor, ROM (Read Only Memory), RAM (Random Access Memory), hard disk unit, display unit, keyboard, mouse, etc. be.
  • a computer program is stored in the RAM or hard disk unit.
  • the at least one device described above achieves its functions by the microprocessor operating according to a computer program.
  • a computer program is configured by combining a plurality of instruction codes indicating instructions to a computer in order to achieve a predetermined function.
  • a part or all of the components constituting at least one of the above devices may be composed of one system LSI (Large Scale Integration).
  • a system LSI is a super-multifunctional LSI manufactured by integrating multiple components onto a single chip, and specifically, it is a computer system that includes a microprocessor, ROM, RAM, etc. .
  • a computer program is stored in the RAM. The system LSI achieves its functions by the microprocessor operating according to a computer program.
  • An IC card or module is a computer system composed of a microprocessor, ROM, RAM, etc.
  • the IC card or module may include the above-mentioned super multifunctional LSI.
  • An IC card or module achieves its functions by a microprocessor operating according to a computer program. This IC card or this module may be tamper resistant.
  • the present disclosure may be the method described above. Furthermore, it may be a computer program that implements these methods using a computer, or it may be a digital signal formed from a computer program.
  • the present disclosure also provides a method for storing computer programs or digital signals on computer-readable recording media, such as flexible disks, hard disks, CD (Compact Disc)-ROMs, DVDs, DVD-ROMs, DVD-RAMs, and BDs (Blu-ray). (Registered Trademark) Disc), semiconductor memory, etc. Further, it may be a digital signal recorded on these recording media.
  • computer-readable recording media such as flexible disks, hard disks, CD (Compact Disc)-ROMs, DVDs, DVD-ROMs, DVD-RAMs, and BDs (Blu-ray). (Registered Trademark) Disc), semiconductor memory, etc. Further, it may be a digital signal recorded on these recording media.
  • the present disclosure may transmit a computer program or a digital signal via a telecommunication line, a wireless or wired communication line, a network typified by the Internet, data broadcasting, or the like.
  • the program or digital signal may be implemented by another independent computer system by recording it on a recording medium and transferring it, or by transferring the program or digital signal via a network or the like.
  • candidate control points that can suppress overshoot or hunting phenomena can be quantitatively evaluated for control problems of time-series data.
  • the evaluation device of the present disclosure has the effect of being able to quantitatively and efficiently search for optimal process conditions even when there are constraints on product characteristics due to standard ranges, and can be used not only in mass production sites of industrial products but also in spatial simulation and dynamic It can also be applied to optimal control applications such as object trajectory control.
  • Evaluation value data 100 Evaluation device 101a Input unit 101b Communication unit 102 Arithmetic circuit 103 Memory 104 Display unit 105 Storage unit 200 Program 201 Characteristic point data 210 Setting information 211 Process condition data 212 Objective data 213 Constraint condition data 214 Area reduction rule data 221 Candidate control point data 222 Control result data 223 Predicted distribution data 224 Evaluation value data 300 Reception image 310 Process condition area 311 to 314, 321 to 328 Input field 320 Product Characteristic area 330 Area reduction rule area

Landscapes

  • Engineering & Computer Science (AREA)
  • Physics & Mathematics (AREA)
  • Theoretical Computer Science (AREA)
  • General Physics & Mathematics (AREA)
  • Software Systems (AREA)
  • Evolutionary Computation (AREA)
  • General Engineering & Computer Science (AREA)
  • Artificial Intelligence (AREA)
  • Data Mining & Analysis (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Mathematical Physics (AREA)
  • Computing Systems (AREA)
  • Business, Economics & Management (AREA)
  • Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Probability & Statistics with Applications (AREA)
  • General Business, Economics & Management (AREA)
  • Economics (AREA)
  • Human Resources & Organizations (AREA)
  • Marketing (AREA)
  • Automation & Control Theory (AREA)
  • Strategic Management (AREA)
  • Algebra (AREA)
  • Computational Mathematics (AREA)
  • General Health & Medical Sciences (AREA)
  • Mathematical Analysis (AREA)
  • Mathematical Optimization (AREA)
  • Pure & Applied Mathematics (AREA)
  • Primary Health Care (AREA)
  • Manufacturing & Machinery (AREA)
  • Tourism & Hospitality (AREA)
  • Computer Hardware Design (AREA)
  • Geometry (AREA)
  • Feedback Control In General (AREA)

Abstract

評価装置（１００）は、第１時刻における制御済みの制御点に対応する既知の特性点に基づき、第１時刻に続く第２時刻における複数の候補制御点に対応する未知の特性点をベイズ最適化によって評価する装置であって、第１時刻における制御済みの制御点および第１時刻における既知の特性点を示す制御結果データ（２２２）と、最適化目的を示す目的データ（２１２）と、制約条件を示す制約条件データ（２１３）と、領域縮小規則データ（２１４）を取得する受信制御部（１０）と、制御結果データ（２２２）、目的データ（２１２）、制約条件データ（２１３）および領域縮小規則データ（２１４）に基づき、複数の未知の特性点それぞれの評価値を算出する評価値算出部（１２）と、評価値を出力する評価値出力部（１３）とを備える。

Description

評価装置、評価方法、およびプログラム
　本開示は、一般の工業製品の量産生産において、製品特性値を規格内へ効率的に制御するための技術に関する。
　工業製品の量産現場において、工業製品の良否を判断する規格など、要求される製品特性の要件を満たすように、設定したプロセス条件を最適な条件で制御する必要がある。例えば、電池の極板の量産現場（生産工程）では、プロセス条件としてスラリーの温度、粘度及び流量、ポンプ回転数、室温等設定され、製品特性として塗布重量等が設定される。
　プロセス条件の最適解は、プロセス条件と製品特性との間の関係性が物理式で表すことができる場合、数学的最適化手法で探索することができることが多い。しかし、その関係性が未知である場合、プロセス条件の値の組み合わせ（すなわち制御点）を一組選択して、実際の制御（すなわち生産）を行う。そして、制御結果として、その制御点に対応する製品特性の値の組み合わせ（すなわち特性点）を獲得する。このような制御を繰り返すことでプロセス条件の最適解を探索することができる。
　一般に、複雑なモノづくりまたは多様な外乱の影響を受ける量産現場では、変化への対応が遅れたり、変化への対処の精度が低かったりする場合、不良品を大量に生産してしまい、多大な金銭的または時間的コストを費やすといった生産ロスが発生する。よって、生産ロスの少ない量産業務遂行のためには、迅速かつ高精度に最適解（すなわち最適なプロセス条件）を探索することが重要になる。
　ところで、従来、量産現場の時系列データの制御では、現場作業者による勘または経験のアプローチ、古典制御または現代制御理論を用いたアプローチがその最適解の探索に用いられてきた。しかし、勘または経験のアプローチは、現場作業者の力量に依存する。また、古典制御または現代制御理論を用いたアプローチは、最適解の探索段階であるパラメータ調整等の段階において、解析者の試行錯誤を必要とする。なお、このアプローチの場合、パラメータ調整次第では、特性値が目標値を行き過ぎてしまうオーバーシュート現象、または、特性値が目標値付近で振動するハンチング現象が発生してしまう。例えば特許文献１には、現代制御手法の一つであるＰＩＤ制御を用いて、ゲイン値の調整をシステマティックに実行する方法が開示されている。
　また、近年、機械学習の分野において、ベイズ最適化を用いたデータ駆動型アプローチが注目されている（例えば非特許文献１および非特許文献２参照）。ベイズ最適化は、入出力の間の対応関係を表現する数理モデルとしてガウス過程を仮定した最適化手法である。制御を行うたびに入出力の関係が一組得られる系においてベイズ最適化を用いる場合、制御結果が得られるごとに、既知の制御結果から算出することができる対応関係の同時分布をもとに、既知の制御結果で周辺化した予測分布を算出する。さらに、獲得関数と呼ばれる評価基準とを用いて、最適な次の制御条件を選択する。これにより、解析者の力量に依らない定量的な評価が可能になり、最適解探索作業の自動化にも貢献できる。
特開２０２１－１１１０１７号公報
M.Emmerich, A.Deutz, J.W.Klinkenberg, "The computation of the expected improvement in dominated hypervolume of Pareto front approximations,"Repport Technique, Leiden University, Vol.34, 2008． M.Abdolshah, A.Shilton, S.Rana, S.Gupta, S.Venkatesh, "Expected Hypervolume Improvement with Constraints," International Conference on Pattern Recognition (ICPR), 2018．
　本開示の一態様に係る評価装置は、第１時刻における制御済みの制御点に対応する既知の特性点に基づいて、前記第１時刻に続く第２時刻における複数の候補制御点に対応する複数の未知の特性点をベイズ最適化によって評価する評価装置であって、前記第１時刻における制御済みの制御点および前記第１時刻における既知の特性点を示す制御結果データを取得する第１受信手段と、前記未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有し、当該最適化目的を示す目的データを取得する第２受信手段と、前記少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データを取得する第３受信手段と、少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、前記分割方法により分割された特性空間の領域ごとに示す領域縮小規則データを取得する第４受信手段と、前記制御結果データ、前記目的データ、前記制約条件データおよび前記領域縮小規則データに基づいて、前記複数の未知の特性点それぞれの評価値を算出する算出手段と、前記評価値を出力する出力手段と、を備え、前記算出手段は、前記制約条件の適合度合いに応じた重み付けを、前記少なくとも１つの製品特性に対する評価値に付与する。
図１は、実施の形態に係る評価装置の概略的な動作を説明するための図である。 図２は、実施の形態に係る各候補制御点および各特性点のそれぞれをグラフで表した一例を示す図である。 図３Ａは、実施の形態に係る部分行列を示す図である。 図３Ｂは、実施の形態に係る部分行列を示す図である。 図４は、実施の形態に係る評価装置の構成を示す図である。 図５は、実施の形態に係る演算回路の機能構成を示すブロック図である。 図６は、実施の形態に係る設定情報の入力を受け付けるために表示部に表示される受付画像の一例を示す図である。 図７は、実施の形態に係る領域縮小規則データの入力を受け付けるために表示部に表示される受付画像の一例を示す図である。 図８は、実施の形態に係る領域縮小規則データの一例を示す図である。 図９Ａは、実施の形態に係る規格範囲の一例を示す図である。 図９Ｂは、実施の形態に係る規格範囲の他の例を示す図である。 図１０は、実施の形態に係る評価装置の処理動作を示すフローチャートである。 図１１Ａは、実施の形態に係る候補制御点データの一例を示す図である。 図１１Ｂは、実施の形態に係る候補制御点データの他の例を示す図である。 図１２は、実施の形態に係る制御結果データの一例を示す図である。 図１３は、実施の形態に係る評価値算出部による処理を説明するための図である。 図１４は、実施の形態に係る予測分布データの一例を示す図である。 図１５Ａは、実施の形態に係る改善領域の一例を示す図である。 図１５Ｂは、実施の形態に係る改善領域の他の例を示す図である。 図１６Ａは、実施の形態に係る、改善領域の体積の算出方法を説明するための図である。 図１６Ｂは、実施の形態に係る、特性空間の全体を複数の小領域に分割する例を示す図である。 図１６Ｃは、実施の形態に係る小領域の下端点および上端点の一例を示す図である。 図１７は、実施の形態に係る領域縮小規則が適用される場合に領域分割される特性空間の一例を示す図である。 図１８Ａは、実施の形態に係る領域縮小規則が適用される場合に領域分割される特性空間の別の例を示す図である。 図１８Ｂは、実施の形態に係る領域縮小規則が適用される場合に領域分割される特性空間の別の例を示す図である。 図１８Ｃは、実施の形態に係る領域縮小規則が適用される場合に領域分割される特性空間の別の例を示す図である。 図１９は、実施の形態に係る領域縮小規則が適用されるパレート境界の算出方法を説明するための図である。 図２０は、実施の形態に係る領域縮小規則が適用されるパレート境界の算出方法を説明するための図である。 図２１は、実施の形態に係る領域縮小規則が適用されるパレート境界の算出方法を説明するための図である。 図２２は、実施の形態に係る領域縮小規則が適用されるパレート境界の算出方法を説明するための図である。 図２３は、実施の形態に係る制約条件がない場合のパレート境界の一例を示す図である。 図２４は、図２３に示されるパレート境界のもとでの改善領域の一例を示す図である。 図２５は、実施の形態に係る制約条件がある場合のパレート境界の一例を示す図である。 図２６は、図２５に示されるパレート境界のもと制約条件がある場合に定められる改善領域の一例を示す図である。 図２７は、実施の形態に係る規格範囲および管理範囲の一例を示す図である。 図２８は、実施の形態に係る評価値データの一例を示す図である。 図２９は、実施の形態に係る表示部に表示される変更後の評価値データの一例を示す図である。 図３０は、実施の形態に係るカルマンフィルタによる予測分布の導出を概念的に示す図である。 図３１は、実施の形態に係る複数の候補制御点の予測分布と製品特性の規格範囲との関係を概念的に示す図である。 図３２Ａは、実施の形態に係る最適化目標値を規格中央に設定する制御方法１の制御イメージを説明するための図である。 図３２Ｂは、実施の形態に係る最適化目標値を規格中央に設定する制御方法１の制御イメージを説明するための図である。 図３３は、オーバーシュート現象およびハンチング現象を概念的に示す図である。 図３４は、カルマンフィルタの特性であるノイズキャンセリング効果を概念的に示す図である。 図３５Ａは、実施の形態に係る最適化目標値を規格上下限に設定する制御方法２の制御イメージを説明するための図である。 図３５Ｂは、実施の形態に係る最適化目標値を規格上下限に設定する制御方法２の制御イメージを説明するための図である。 図３６Ａは、実施の形態に係る最適化目標値を徐々に規格中央に移動させていく制御方法３の制御イメージを説明するための図である。 図３６Ｂは、実施の形態に係る最適化目標値を徐々に規格中央に移動させていく制御方法３の制御イメージを説明するための図である。 図３７は、実施の形態に係る製品特性数が２次元以上の場合に最適化目標値の候補が複数存在することを説明するための図である。 図３８Ａは、実施の形態に係る中心が規格中央と一致しない場合の結合アクティブ領域を示す図である。 図３８Ｂは、実施の形態に係る中心が規格中央と一致する場合の結合アクティブ領域を示す図である。 図３９は、実施の形態の実施例に係る最適解の探索を行う際に得られた制御結果データシートの一例を示す図である。 図４０Ａは、図３９の２６番目までの制御結果を得た時点での結合アクティブ領域を示す図である。 図４０Ｂは、図３９の２７番目の制御結果を得た時点での結合アクティブ領域を示す図である。
　（本発明の基礎となった知見）
　本発明者は、「背景技術」の欄において記載した上記非特許文献１および非特許文献２に関し、以下の課題が生じることを見出した。
　複数の製品特性を同時に最適化する多目的ベイズ最適化に関する手法がいくつか提案されている。例えば、上記非特許文献１では、その一種であるＥＨＶＩ（Expected Hypervolume Improvement）の最適解探索原理および具体的な算出方法が開示されている。これにより、最適化したい製品特性が複数存在していても、最適解探索の定量評価が可能になる。
　また、工業製品開発または製造プロセス開発において、製品特性の値に関する制約条件として規格範囲が設けられている場合がある。例えば、「電池容量は１８５０～１９５０［ｍＡｈ］の規格範囲に収まってほしい」、「寿命は３年を最小値として長いほどよい（すなわち、寿命の規格範囲の最小値は３年で、最大値は＋∞）」等の規格範囲である。規格範囲がある場合に、従来のベイズ最適化を適用すると、計算効率の悪い探索をしたり、最適解ではない、別の領域に探索が進んだりする恐れがある。
　そこで、制約条件付きベイズ最適化に関する手法もいくつか提案されている。例えば、上記特許文献１の手法では、各候補制御点について、ガウス過程回帰により求まった予測分布から、規格範囲に収まる確率を算出し、その確率が或るしきい値を上回る候補制御点のみを抽出して、獲得関数を評価する。これにより、制約条件付き最適化問題に対応している。
　また、例えば、上記非特許文献２では、ＥＨＶＩを制約条件がある場合に拡張したＥＨＶＩＣ（Expected Hypervolume Improvement with Constraints）が開示されている。上記非特許文献２に記載の獲得関数の設計手法は、規格範囲内に入る確率と改善量とを総合的に指標化し、すべての候補制御点について評価する。そのため、探索効率が向上する（すなわち、真の最適解が求まる）可能性が高い。
　しかしながら、上記非特許文献２の獲得関数の設計手法では、その獲得関数を評価する上で、最大化または最小化したい製品特性の数と暫定の最適解（非劣解）であるパレート点の数に関して、計算量が指数関数オーダーで増大してしまうという問題もある。
　また、特許文献１には、上述したように、現代制御手法の一つであるＰＩＤ制御を用いて、ゲイン値の調整をシステマティックに実行する方法が開示されている。しかし、特許文献１に開示される方法はルールベースでのゲイン調整方法であるため定量評価が難しく、量産現場の時系列データの制御においてオーバーシュート現象またはハンチング現象を抑制できる良好な制御結果が得られるとは限らない。このため、量産現場の時系列データの制御において、パラメータの評価が人依存になってしまう。
　一方、時系列で変化する状態を推定するための計算手法として、カルマンフィルタが広く知られている。カルマンフィルタは、状態空間モデルと呼ばれる数理モデルにより、システム内部の見えない状態を推定する計算手法である。このため、カルマンフィルタは、推定した状態の情報をもとに量産現場の時系列データの制御に応用され得る。
　以上を鑑みて、本開示では、カルマンフィルタによる状態推定と、ベイズ最適化による獲得関数評価のアイデアとを利用することを考える。
　しかしながら、単にカルマンフィルタとベイズ最適化とを組み合わせるだけでは、オーバーシュート現象またはハンチング現象が発生してしまうことを抑制できない。
　そこで、本開示の評価装置は、時系列データの制御問題に対して、オーバーシュート現象またはハンチング現象の発生を抑制できる候補制御点を定量的に評価できることを目的とする。
　本開示の一態様に係る評価装置は、第１時刻における制御済みの制御点に対応する既知の特性点に基づいて、前記第１時刻に続く第２時刻における複数の候補制御点に対応する複数の未知の特性点をベイズ最適化によって評価する評価装置であって、前記第１時刻における制御済みの制御点および前記第１時刻における既知の特性点を示す制御結果データを取得する第１受信手段と、前記未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有し、当該最適化目的を示す目的データを取得する第２受信手段と、前記少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データを取得する第３受信手段と、少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、前記分割方法により分割された特性空間の領域ごとに示す領域縮小規則データを取得する第４受信手段と、前記制御結果データ、前記目的データ、前記制約条件データおよび前記領域縮小規則データに基づいて、前記複数の未知の特性点それぞれの評価値を算出する算出手段と、前記評価値を出力する出力手段と、を備え、前記算出手段は、前記制約条件の適合度合いに応じた重み付けを、前記少なくとも１つの製品特性に対する評価値に付与する。
　これにより、算出手段が、制御結果データ、目的データ、制約条件データおよび領域縮小規則データに基づいて、第１時刻に続く第２時刻における複数の未知の特性点の評価値を算出する。複数の未知の特性点の評価値を算出するときには、制約条件の適合度合いに応じた重み付けを、少なくとも１つの製品特性に対する評価値に付与し、この少なくとも１つの製品特性は最適化目的を有する。したがって、最適化問題の目的を有する時系列の製品特性に対して制約条件が付与されている最適化問題に対して、カルマンフィルタとベイズ最適化を適用することができる。このようにして、過去の制御結果情報である第１時刻における制御済みの制御点に対応する既知の特性点から、複数の未知の特性点の評価値を算出できるので、次に設定すべき候補制御点を定量的に解析することができる。
　以上のように、本発明の解析システムによれば、次に設定すべき最適な候補制御点を定量的に算出し、解析者の力量に依らず、高精度な制御の実現が期待できる。
　また、前記制約条件は、少なくとも１つの制約範囲であってもよい。前記最適化目的には、製品特性を前記少なくとも１つの制約範囲のうちの何れかの制約範囲内に収める第１目的と、製品特性を最小化または最大化する第２目的とがあってもよい。前記算出手段は、少なくとも１つの製品特性のそれぞれについて、以下の（ｉ）～（ｉｉｉ）の場合で互いに異なる重み付け処理を行うことによって、前記評価値を算出してもよい。
　（ｉ）評価値を算出するために用いられる当該製品特性の区間が前記少なくとも１つの制約範囲のそれぞれの外にある場合。
　（ｉｉ）前記区間が前記少なくとも１つの制約範囲のうちの何れかの制約範囲内であって、かつ、前記最適化目的が前記第１目的である場合。
　（ｉｉｉ）前記区間が前記少なくとも１つの制約範囲のうちの何れかの制約範囲内であって、かつ、前記最適化目的が前記第２目的である場合。
　例えば、複数の制約範囲は、規格範囲と、その規格範囲に含まれる管理範囲である。そして、（ｉｉ）の場合は、区間が規格範囲内および管理範囲外であって、かつ、最適化目的が第１目的である第１の場合と、区間が管理範囲内であって、かつ、最適化目的が第１目的である第２の場合とに、場合分けされる。また、例えば、第１の場合よりも第２の場合の方が大きい重みを用いた重み付け処理が行われる。このように、複数の制約範囲があり、（ｉｉ）の場合をさらに複数の場合に場合分けすることによって、複数の制約範囲のそれぞれに段階的に重み付けを行うことができる。したがって、製品特性の値が規格範囲に収まり、可能な限り管理範囲に収まって欲しいような場合であっても、評価値を適切に算出することができる。その結果、最適化問題に対する適用場面をさらに拡張することができる。
　また、複数のプロセス条件のそれぞれの所定の条件を満たす値を組み合わせることによって、前記複数の候補制御点を作成する候補制御点作成手段をさらに備えてもよい。
　例えば、所定の条件は、複数のプロセス条件のそれぞれの比率変数の値の和が１であるという条件である。より具体的な一例では、その比率変数は、プロセス条件に対応する化合物などの材料の配合比である。したがって、複数種の化合物の配合比の組み合わせごとに、その組み合わせに対する評価値を算出することができる。その結果、それらの化合物の配合によって得られる合成材料の少なくとも１つの製品特性に対する最適解を適切に探索することができる。
　また、前記算出手段は、前記少なくとも１つの制約範囲のうち、矩形と異なる形状の制約範囲に基づいて、前記評価値を算出してもよい。
　これにより、例えば２つの製品特性によって表現される特性空間において、円形、楕円形、星形などの制約範囲に基づいて評価値が算出されるため、制約範囲の形状が矩形の場合に限定されることなく、適用場面をさらに拡張することができる。
　また、前記算出手段は、カルマンフィルタを用いて前記複数の候補制御点における予測分布を算出し、算出した前記予測分布を用いて、前記評価値を算出してもよい。
　これにより、カルマンフィルタで求まる予測分布と予測分布からベイズ最適化を用いて求まる評価値とを逐次的に算出することができる。よって、次に制御点として設定する候補制御点を評価値に基づいて選択できる。つまり、各候補制御点に対して算出された評価値が出力されるため、評価装置のユーザは、それらの評価値に基づいて候補制御点を次の制御点として選択し、その制御点を用いた制御によって得られる特性点を、次の各候補制御点の評価値の算出に利用することができる。このような制御と評価値の算出および出力との繰り返しによって、各製品特性の最適化目的を満たす候補制御点の解、すなわち最適解を得ることができる。
　また、前記算出手段は、モンテカルロ法を用いて前記評価値を算出してもよい。
　これにより、モンテカルロ法が近似手法であるため、評価値が解析的に算出困難な場合でも、近似的に算出することができる。
　以下、実施の形態について、図面を参照しながら具体的に説明する。
　なお、以下で説明する実施の形態は、いずれも包括的または具体的な例を示すものである。以下の実施の形態で示される数値、形状、材料、構成要素、構成要素の配置位置および接続形態、ステップ、ステップの順序などは、一例であり、本開示を限定する主旨ではない。また、以下の実施の形態における構成要素のうち、独立請求項に記載されていない構成要素については、任意の構成要素として説明される。また、各図は、模式図であり、必ずしも厳密に図示されたものではない。また、各図において、同じ構成部材については同じ符号を付している。
　（実施の形態１）
　［概要］
　図１は、本実施の形態に係る評価装置の概略的な動作を説明するための図である。
　本実施の形態における評価装置１００は、複数の候補制御点のそれぞれに対する評価値を算出し、それらの評価値を示す評価値データ２２４を表示する。候補制御点は、制御点の候補とされる点である。制御点は、制御条件（制御空間上における各プロセス条件の値の組み合わせ）を示す制御空間上における点である。評価値は、その候補制御点にしたがった制御によって得られると予測される製品特性の評価結果を示す値である。例えば、評価値は、制御によって得られると予測される製品特性が最適化目的に合致している度合いを示し、評価値が大きいほど、その度合いは大きい。
　ユーザは、その評価値データ２２４によって示される各候補制御点の評価値を参照し、それらの候補制御点のうちの１つを次の制御点として選択する。ユーザは、その選択された制御点にしたがった制御を、制御設備（量産設備）を用いて行う。制御によって、その制御点に対応する特性点が得られる。特性点は、例えば製品特性の値を示し、複数の製品特性があれば、複数の製品特性の値の組み合わせとして示される。ユーザは、その得られた特性点を制御点に対応付けて評価装置１００に入力する。その結果、評価装置１００は、その制御によって得られた特性点を用いて、各候補制御点に対する評価値を再び算出し、それらの評価値を示す評価値データ２２４を再表示する。つまり、評価値データ２２４が更新される。評価装置１００は、このような評価値データ２２４の更新を繰り返すことによって、製品特性の最適解を探索する。
　図２は、各候補制御点および各特性点のそれぞれをグラフで表した一例を示す図である。具体的には、図２の（ａ）のグラフは、制御空間に配置される各候補制御点を示し、図２の（ｂ）のグラフは、特性空間に配置される各特性点を示す。
　制御空間における候補制御点は、図２の（ａ）に示すように、第１プロセス条件および第２プロセス条件の値の組み合わせに対応する格子点上に配置される。図２の（ａ）に示す各候補制御点に対応する特性点は、図２の（ｂ）に示すように、特性空間に配置される。具体的には、候補制御点が制御点として選択され、その制御点にしたがった制御によって第１製品特性および第２製品特性のそれぞれの値が得られる場合、その制御点に対応する特性点は、その第１製品特性の値と第２製品特性の値との組み合わせによって表現される位置に配置される。ここで、候補制御点と特性点との間には１対１の対応関係があるが、その対応関係（すなわち、図２中の関数ｆ）は未知である。
　制御を一回実行するとは、一つの候補制御点を選択し、その選択された候補制御点に対応する特性点との対応関係を一組獲得すること、として換言できる。
　また、図２の（ｂ）に示すように、設定された規格範囲により、特性空間は規格範囲内領域と規格範囲外領域とに分割される。また、本実施の形態では、最適化目的を有する製品特性に対して、規格範囲である制約条件が付与されていてもよい。制約条件は、製品特性に付与される条件であって、例えば、製品特性の値の範囲を条件として指定する制約範囲がある。制約範囲として、例えば、製品特性の規格によって定められる規格範囲や、ユーザが適宜設定可能な管理範囲などがある。
　なお、本実施の形態では、第１プロセス条件および第２プロセス条件のように、プロセス条件の数が２つであり、かつ、第１製品特性および第２製品特性のように、製品特性の数が２つである例について主に説明するが、プロセス条件の数および製品特性の数は２つに限らない。プロセス条件の数は１つであっても、３つ以上であってもよく、製品特性の数は１つであってもよく、３つ以上であってもよい。また、プロセス条件の数と製品特性の数とは等しくても異なっていてもよい。
　本実施の形態では、その対応関係は普遍的なものである必要はなく、一時点前の制御結果に依存し、時系列で変化すると仮定する。
　評価装置１００では、候補制御点と特性点との対応関係をカルマンフィルタで記述する。
　以下、評価装置１００で利用されるカルマンフィルタについて説明する。
　＜カルマンフィルタ＞
　カルマンフィルタは、状態空間モデルと呼ばれる数理モデルにより、システム内部の見えない状態を推定する計算手法である。換言すると、カルマンフィルタは、誤差を含んだ観測値から、時間的に変化するシステム内部の状態量を推定する計算手法であり、ベイズ統計に基づく枠組みに含まれる。
　カルマンフィルタの数理モデルは、時刻tにおけるシステムの内部状態量をX(t)、時刻tにおける観測量をY(t)としたとき、時間的な変化するシステムの内部状態量を（式１）で表現することができる。
　ここで、Ａ，Ｂ，Ｃは変換を定める行列を表す。v(t),w(t)は時刻tにおけるガウシアンノイズを表す。ガウシアンノイズの平均および分散は、例えば０と１とに設定する等、解析者が適宜設定することができる。
　なお、（式１）の上式は、状態方程式とも称され、観測されていないシステムの内部状態の時間的発展を記述している。（式１）の下式は、観測方程式とも称され、システムの内部状態から我々が観測する観測量への変換を記述している。システムの内部状態量がすべて観測されている場合、Ｃを単位行列とすればよい。Ａ，Ｂ，Ｃの各要素の値は通常、未知であるので、自己回帰モデル等の時系列解析手法で逐次的に推定しながら進めることもできる。
　上記の（式１）をもとに定式化されるカルマンフィルタは、線形ガウスカルマンフィルタと称される。線形ガウスカルマンフィルタでは、観測済み状態量から予測される、次時点における内部状態量は、（式２）に示すような予測分布を用いて導出される。
　ここで、Y(1:t-1)は、時刻1から時刻t-1までのYをまとめたベクトルを表す。また、（式２）の右辺における正規分布の平均および分散は、下記の（式３）の５つの更新式を各時刻で解くことで算出される。
　なお、カルマンフィルタには、例えばＥｎＫＦ（Ensemble Kalman Filter）、ＥＫＦ（Extended Kalman Filter）、ＵＫＦ (Unscented Kalman Filter)、粒子フィルタなど、変換が非線形の場合、ノイズが非ガウス分布の場合、外部入力が影響する場合などによる様々なパターンが存在する。つまり、本実施の形態におけるカルマンフィルタの形式は、上記のいずれのパターンであってもよく特に指定されない。また、候補制御点と特性点との対応関係が何らかの確率分布で推定できるのであれば、カルマンフィルタ以外の手法を用いてもよい。
　以下では、説明を簡単にするため、制御点と特性点とをまとめて内部状態として扱う線形ガウス形式のカルマンフィルタを例に挙げて説明する。
　候補制御点ベクトルX=(X1, ….,XDx)と特性点ベクトルY=(Y1,….,YDy)とをまとめてZ=(X1,….,XDx,Y1,….,YDy)と表す場合、（式１）に示される状態方程式と観測方程式とは以下の（式４）のように記述される。
　この（式４）のように定式化されるカルマンフィルタでは、（式３）の第４更新式におけるY(t)に、Z観測(t)を代入すれば、Zの内部状態量の予測分布が逐次的に算出される。
　ここで、多次元正規分布は、一部の次元で条件付け操作をしても正規性が保存される性質を有している。この性質を用いることで、プロセス条件としての制御因子Xが与えられたもとでの製品特性としての目的特性Ｙの予測分布が、正規分布として導出されることになる。その正規分布の平均および分散は、具体的には、（式５）で与えられる。
　ここで、
および
は、（式３）を用いて算出される量であるとし、AY,P－ XX,P－ XY,P－ YY,P－ YXは、図３Ａ及び図３Ｂのように、行列の一部の要素を抽出した部分行列を表すとする。図３Ａ及び図３Ｂは、部分行列を示す図である。図３Ａには、変換を定める行列Ａの要素とその部分行列AYとが示されている。図３Ｂには、行列P－の要素とその部分行列P－ XX,P－ XY,P－ YY,P－ YXとが示されている。
　［ハードウェア構成］
　図４は、本実施の形態に係る評価装置１００の構成を示す図である。
　評価装置１００は、入力部１０１ａ、通信部１０１ｂ、演算回路１０２、メモリ１０３、表示部１０４、および、記憶部１０５を備える。
　入力部１０１ａは、ユーザによる入力操作を受け付けるＨＭＩ（Human Machine Interface）である。入力部１０１ａは、例えばキーボード、マウス、タッチセンサ、タッチパッドなどである。
　例えば、入力部１０１ａは、ユーザからの入力として、設定情報２１０を受け付ける。設定情報２１０は、プロセス条件データ２１１、目的データ２１２、制約条件データ２１３および領域縮小規則データ２１４を含む。プロセス条件データ２１１は、例えば、図２の（ａ）に示すように、プロセス条件の取り得る値を示すデータである。プロセス条件の値は、連続値でも離散値でもよい。目的データ２１２は、例えば、最小化／最大化などの製品特性の最適化目的を示すデータである。制約条件データ２１３は、例えば、制約範囲などの制約条件を示すデータである。領域縮小規則データ２１４は、パレート境界を算出する規則を示すデータであり、改善量の算出方法を変更させる。より具体的には、領域縮小規則データ２１４は、少なくとも２つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、分割方法により分割された特性空間の領域ごとに示す。詳細については後述する。
　通信部１０１ｂは、他の機器と有線または無線で接続し、他の機器とデータを送受信する。例えば、通信部１０１ｂは、他の装置（例えば、制御装置）から特性点データ２０１を受信する。
　表示部１０４は、画像または文字などを表示する。表示部１０４は、例えば液晶ディスプレイ、プラズマディスプレイ、有機ＥＬ（Electro-Luminescence）ディスプレイなどである。なお、表示部１０４は、入力部１０１ａと一体となっているタッチパネルでもよい。
　記憶部１０５は、演算回路１０２への各命令が記述されたプログラム（すなわちコンピュータプログラム）２００および各種データを格納している。記憶部１０５は、不揮発性の記録媒体であって、例えば、ハードディスクなどの磁気記憶装置、ＳＳＤ（Solid State Drive）などの半導体メモリ、光ディスクなどである。
　なお、プログラム２００および各種データは、例えば、上述の他の機器から通信部１０１ｂを介して評価装置１００に提供され、記憶部１０５に格納されてもよい。記憶部１０５は、各種データとして、候補制御点データ２２１、制御結果データ２２２、予測分布データ２２３、評価値データ２２４を格納する。
　候補制御点データ２２１は、各候補制御点を示すデータである。図２の（ａ）に示す例では、各候補制御点は、第１プロセス条件および第２プロセス条件の値の組み合わせによって表現される。候補制御点データ２２１は、第１プロセス条件および第２プロセス条件の値の組み合わせが列挙されたテーブル形式のデータであってもよい。このような候補制御点データ２２１の具体例については、図１１Ａおよび図１１Ｂを用いて詳細に説明する。
　制御結果データ２２２は、制御に用いられた１以上の制御点と、その１以上の制御点のそれぞれに対応する特性点とを示すデータである。例えば、制御結果データ２２２は、図２の（ａ）の制御空間上の制御点と、その制御点を用いた制御によって得られた図２の（ｂ）の特性空間上の特性点との組み合わせを示す。制御点は、第１プロセス条件および第２プロセス条件の値の組み合わせによって表現される。特性点は、第１製品特性および第２製品特性の値の組み合わせによって表現される。制御結果データ２２２は、その制御点および特性点の組み合わせが列挙されたテーブル形式のデータであってもよい。この制御結果データ２２２の具体例については、図１２を用いて詳細に説明する。
　予測分布データ２２３は、候補制御点データ２２１によって示される全候補制御点の予測分布を示すデータである。予測分布は、上述したようにカルマンフィルタに基づく（式３）によって求まる分布であって、例えば、平均および分散によって表現される。例えば、予測分布データ２２３は、各候補制御点に対して、第１製品特性の予測分布と、第２製品特性の予測分布とを対応付けて示すテーブル形式のデータであってもよい。この予測分布データ２２３の具体例については、図１４を用いて詳細に説明する。
　評価値データ２２４は、例えば図１に示すように、複数の候補制御点のそれぞれに対する評価値を示すデータである。例えば、評価値データ２２４は、複数の候補制御点のそれぞれに評価値を関連付けて示すテーブル形式のデータであってもよい。この評価値データ２２４の他の具体例については、図２８等を用いて詳細に説明する。
　演算回路１０２は、記憶部１０５からプログラム２００をメモリ１０３に読み出し、展開されたプログラム２００を実行する回路である。演算回路１０２は、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphics Processing Unit）などである。
　［機能構成］
　図５は、演算回路１０２の機能構成を示すブロック図である。
　演算回路１０２は、プログラム２００を実行することによって、評価値データ２２４を生成するための複数の機能を実現する。具体的には、演算回路１０２は、受信制御部（第１受信手段、第２受信手段、第３受信手段および第４受信手段）１０、候補制御点作成部（候補制御点作成手段）１１、評価値算出部（算出手段）１２、および評価値出力部（出力手段）１３を備える。
　＜受信制御部１０＞
　受信制御部１０は、入力部１０１ａまたは通信部１０１ｂを介して、特性点データ２０１、プロセス条件データ２１１、目的データ２１２、制約条件データ２１３および領域縮小規則データ２１４を受信する。例えば、受信制御部１０は、ユーザによる入力部１０１ａへの入力操作によって特性点データ２０１が入力されると、その特性点データ２０１に示される特性点を制御点に対応付けて記憶部１０５の制御結果データ２２２に書き込む。これにより、制御結果データ２２２が更新される。この制御結果データ２２２が更新されると、受信制御部１０は、評価値算出部１２に対して、その更新後の制御結果データ２２２を用いた処理を実行させる。つまり、受信制御部１０は、評価値算出部１２に対して評価値の算出を実行させる。なお、このときには、評価値算出部１２は、記憶部１０５に既に格納されている候補制御点データ２２１を用いて評価値の算出を実行する。このように、受信制御部１０は、特性点データ２０１の入力をトリガーに、評価値算出部１２に対して評価値の算出を開始させる。
　また、受信制御部１０は、他のトリガーに応じて、評価値算出部１２に対して評価値の算出を開始させてもよい。例えば、制御結果データ２２２が記憶部１０５に既に格納されていれば、受信制御部１０は、ユーザによる制御点の水準の入力をトリガーに、評価値算出部１２に対して評価値の算出を開始させてもよい。なお、制御点の水準は、例えば、プロセス条件が取り得る値の最小値、最大値、および離散幅などである。つまり、受信制御部１０は、ユーザによって制御点の水準が入力され、その水準に基づいて候補制御点データ２２１が生成されると、その候補制御点データ２２１と制御結果データ２２２と領域縮小規則データ２１４とに基づく評価値の算出を評価値算出部１２に対して開始させる。
　あるいは、候補制御点データ２２１が記憶部１０５に既に格納されていれば、受信制御部１０は、ユーザによる制御結果データ２２２の入力をトリガーに、評価値算出部１２に対して評価値の算出を開始させてもよい。受信制御部１０は、ユーザによって制御結果データ２２２が入力されると、その制御結果データ２２２と候補制御点データ２２１と領域縮小規則データ２１４とに基づく評価値の算出を評価値算出部１２に対して開始させる。
　あるいは、候補制御点データ２２１が記憶部１０５に既に格納されていれば、受信制御部１０は、通信部１０１ｂによる制御結果データ２２２の受信をトリガーに、評価値算出部１２に対して評価値の算出を開始させてもよい。例えば、制御設備、制御装置、または製造装置などが制御結果データ２２２を評価装置１００に送信し、通信部１０１ｂがその制御結果データ２２２を受信する。受信制御部１０は、通信部１０１ｂによって制御結果データ２２２が受信されると、その制御結果データ２２２と候補制御点データ２２１と領域縮小規則データ２１４とに基づく評価値の算出を評価値算出部１２に対して開始させる。
　このように、受信制御部１０は、候補制御点データ２２１および制御結果データ２２２があれば、それらに基づく評価値の算出を評価値算出部１２に対して開始させる。なお、制御結果データ２２２が記憶部１０５に既に格納されていれば、受信制御部１０は、ユーザによる候補制御点データ２２１の入力をトリガーに、評価値算出部１２に対して評価値の算出を開始させてもよい。また、候補制御点データ２２１および制御結果データ２２２が記憶部１０５に既に格納されていれば、受信制御部１０は、ユーザによる開始指示の入力をトリガーに、評価値算出部１２に対して評価値の算出を開始させてもよい。
　＜候補制御点作成部１１＞
　候補制御点作成部１１は、受信制御部１０によって取得されたプロセス条件データ２１１に基づいて、候補制御点データ２２１を生成する。換言すると、候補制御点作成部１１は、複数のプロセス条件のそれぞれの所定の条件を満たす値を組み合わせることによって、複数の候補制御点を作成する。本実施の形態では、候補制御点作成部１１は、複数の候補制御点のそれぞれを、１以上のプロセス条件のそれぞれの値を用いて作成する。そして、候補制御点作成部１１は、その生成された候補制御点データ２２１を記憶部１０５に格納する。
　＜評価値算出部１２＞
　評価値算出部１２は、候補制御点データ２２１および制御結果データ２２２を記憶部１０５から読み出し、それらのデータに基づいて予測分布データ２２３を生成し、その予測分布データ２２３を記憶部１０５に格納する。さらに、評価値算出部１２は、その予測分布データ２２３と、受信制御部１０によって取得された目的データ２１２、制約条件データ２１３および領域縮小規則データ２１４とに基づいて、評価値データ２２４を生成し、その評価値データ２２４を記憶部１０５に格納する。
　本実施の形態では、評価値算出部１２は、カルマンフィルタを用いて複数の候補制御点における予測分布を算出し、算出した予測分布を用いて、評価値を算出する。なお、評価値算出部１２は、後述するように、少なくとも１つの制約範囲のうち、矩形と異なる形状の制約範囲に基づいて、評価値を算出してもよい。
　＜評価値出力部１３＞
　評価値出力部１３は、記憶部１０５から評価値データ２２４を読み出して、その評価値データ２２４を表示部１０４に出力する。あるいは、評価値出力部１３は、通信部１０１ｂを介して、評価値データ２２４を外部の装置に出力してもよい。つまり、評価値出力部１３は、各候補制御点の評価値を出力する。なお、評価値出力部１３は、評価値算出部１２から評価値データ２２４を直接取得して、その評価値データ２２４を表示部１０４に出力してもよい。同様に、評価値出力部１３は、記憶部１０５から予測分布データ２２３を読み出して、その予測分布データ２２３を表示部１０４に出力する。なお、評価値出力部１３は、評価値算出部１２から予測分布データ２２３を直接取得して、その予測分布データ２２３を表示部１０４に出力してもよい。
　［入力］
　図６は、設定情報２１０の入力を受け付けるために表示部１０４に表示される受付画像の一例を示す図である。
　受付画像３００は、プロセス条件領域３１０と、製品特性領域３２０とを含む。プロセス条件領域３１０は、プロセス条件データ２１１の入力を受け付けるための領域である。製品特性領域３２０は目的データ２１２および制約条件データ２１３の入力を受け付けるための領域である。
　＜プロセス条件領域３１０＞
　プロセス条件領域３１０は、入力フィールド３１１～３１４を有する。入力フィールド３１１は、第１プロセス条件の名称を入力するためのフィールドである。例えば、入力フィールド３１１には、第１プロセス条件の名称として「Ｘ１」が入力される。入力フィールド３１２は、第１プロセス条件の値を入力するためのフィールドである。例えば、この入力フィールド３１２には、第１プロセス条件の値として「－５，－４，－３，－２，－１，０，１，２，３，４，５」が入力される。同様に、入力フィールド３１３は、第２プロセス条件の名称を入力するためのフィールドである。例えば、入力フィールド３１３には、第２プロセス条件の名称として「Ｘ２」が入力される。入力フィールド３１４は、第２プロセス条件の値を入力するためのフィールドである。例えば、この入力フィールド３１４には、第２プロセス条件の値として「－５，－４，－３，－２，－１，０，１，２，３，４，５」が入力される。
　このような入力フィールド３１１～３１４への入力によって、その入力結果に応じたプロセス条件データ２１１が評価装置１００に入力される。
　＜製品特性領域３２０＞
　製品特性領域３２０は、入力フィールド３２１～３２８を有する。入力フィールド３２１および３２５は、第１製品特性の名称および第２製品特性の名称を入力するためのフィールドである。例えば、入力フィールド３２１には、第１製品特性の名称として「Ｙ１」が入力され、入力フィールド３２５には、第２製品特性の名称として「Ｙ２」が入力される。入力フィールド３２２および３２６は、第１製品特性および第２製品特性の最適化目的を選択するためのフィールドである。具体的には、その入力フィールド３２２および３２６のそれぞれは、「最大化」、「最小化」、および「規格範囲内」のうちの何れか１つを目的として選択するための３つのラジオボタンを有する。目的としての「最大化」は、第１製品特性または第２製品特性の値の最大化を目的とし、「最小化」は、第１製品特性または第２製品特性の値の最小化を目的とする。「規格範囲内」は、第１製品特性または第２製品特性の値が規格範囲内に収まることを目的とする。例えば、ユーザによる入力部１０１ａへの入力操作によって、「規格範囲内」を示すラジオボタンが選択された場合、評価装置１００は、第１製品特性または第２製品特性の最適化目的として規格範囲内を選択する。入力フィールド３２３および３２４は、第１製品特性の規格範囲を示す最小値および最大値をそれぞれ入力するためのフィールドである。例えば、入力フィールド３２３には、最小値として「３０」が入力され、入力フィールド３２４に、最大値として「４０」が入力された場合、評価装置１００は、規格範囲を３０～４０とする。入力フィールド３２７および３２８は、第２製品特性の規格範囲における最小値および最大値をそれぞれ入力するためのフィールドである。例えば、入力フィールド３２７に、規格範囲の最小値として「１０」が入力され、入力フィールド３２８が未入力である場合、評価装置１００は、規格範囲を１０～＋∞とする。なお、入力フィールド３２７が未入力の場合には、評価装置１００は、規格範囲の最小値を－∞とする。
　このような入力フィールド３２１～３２８への入力によって、その入力結果に応じた目的データ２１２および制約条件データ２１３が評価装置１００に入力される。つまり、受信制御部１０は、入力フィールド３２２および３２６への入力に応じた目的データ２１２を取得し、入力フィールド３２３、３２４、３２７、および３２８への入力に応じた制約条件データ２１３を取得する。図６の例では、その目的データ２１２は、第１製品特性の最適化目的として、その第１製品特性の値が規格範囲内に収まることを示し、かつ、第２製品特性の最適化目的として、その第２製品特性の値の最小化を示す。さらに、制約条件データ２１３は、第１製品特性の規格範囲が３０～４０であることを示し、かつ、第２製品特性の規格範囲が１０～＋∞であることを示す。
　＜領域縮小規則領域３３０＞
　図７は、領域縮小規則データ２１４の入力を受け付けるために表示部１０４に表示される受付画像の一例を示す図である。
　受付画像３００は、図６に示すプロセス条件領域３１０と、製品特性領域３２０と以外に、図７に示す領域縮小規則領域３３０を含む。領域縮小規則領域３３０は、領域縮小規則データ２１４の入力を受け付けるための領域である。
　領域縮小規則領域３３０は、入力フィールド３３１および３３２を有する。入力フィールド３３１は、領域縮小規則を適用するかどうかを入力するためのフィールドである。入力フィールド３３２は、領域縮小規則を適用する場合、領域分割に空集合を含むかどうかを入力するためのフィールドである。具体的には、入力フィールド３３１には、領域縮小規則を「適用する」または「適用しない」を選択するための２つのラジオボタンを有する。また、入力フィールド３３２には、領域分割において空集合設定を行うかどうかを示す「有り」または「なし」を選択するための２つのラジオボタンを有する。
　例えば、ユーザによる入力部１０１ａへの入力操作によって、「適用する」を示すラジオボタンが選択された場合、評価装置１００は、評価値の算出処理において領域縮小規則を適用する。さらに、ユーザによる入力部１０１ａへの入力操作によって、「適用する」を示すラジオボタンが選択された場合、評価装置１００は、評価値の算出処理において、空集合を含む領域分割を行った上で領域縮小規則を適用して改善領域の体積を算出する。
　＜領域縮小規則データ２１４の一例＞
　図８は、領域縮小規則データ２１４の一例を示す図である。
　図７の受付画像３００の領域縮小規則領域３３０によって入力された領域縮小規則データ２１４は、例えば図８に示すように、パレート境界の定義が変更されることと、アクティブ領域の定義が変更されることと、領域分割が適用されることとを示す。これにより改善量の算出方法が変更される。なお、変更された定義および適用される領域分割の具体例は後述するため、ここでの説明は省略する。
　なお、プロセス条件データ２１１に含まれるプロセス条件の値は、絶対値であってもよいし、他のプロセス条件の値または全プロセス条件の値の総和に対する比率などの相対値であってもよい。また、第１プロセス条件の連続変数の値と、第２プロセス条件の比率変数の値と、第３プロセス条件の比率変数の値とを示してもよい。変数には、連続変数とは異なる離散変数がある。プロセス条件が離散変数である場合、その離散変数は「リンゴ、ミカン、バナナ」または「触媒あり、触媒なし」のように、大小関係および数値的な大きさを持たない。比率変数は、例えば、第１プロセス条件の材料と、第２プロセス条件の材料とが配合されることによって生成される合成材料において、第１プロセス条件または第２プロセス条件の材料の配合比を示す。
　また、目的データ２１２は、例えば、第１製品特性の最適化目的と、第２製品特性の最適化目的とを示す。また、図６の受付画像３００の製品特性領域３２０によって入力された制約条件データ２１３は、第１製品特性の規格範囲と、第２製品特性の規格範囲とを示していてもよい。具体的には、目的データ２１２には、第１製品特性の最適化目的として「規格範囲内」が示され、第２製品特性の最適化目的として「最小化」が示されていてもよい。また、制約条件データ２１３には、第１製品特性の規格範囲として、最小値「３０」～最大値「４０」の範囲が示され、第２製品特性の規格範囲として、最小値「１０」～最大値「＋∞」の範囲が示されている。したがって、第２製品特性の最適化目的は、最小値「１０」以上の規格範囲内での第２製品特性の値の最小化である。
　＜規格範囲の一例＞
　図９Ａは、規格範囲の一例を示す図である。
　制約条件データ２１３によって示される規格範囲は、例えば図９Ａのように、特性空間上において矩形の範囲で表現される。なお、図９Ａに示す例では、規格範囲の形状は矩形であるが、他の形状であってもよい。つまり、規格範囲の形状は、後述の評価値の算出が実装可能であれば、任意の形状であってもよい。
　図９Ｂは、規格範囲の他の例を示す図である。
　規格範囲は、図９Ｂに示すように、例えば円形であってもよい。具体的な一例では、第１製品特性および第２製品特性の特性空間における規格範囲は、円の中心（２０，２０）と半径１０とによって表現される。なお、その規格範囲の形状は、円形以外の形状であってもよく、楕円形、星形などであってもよい。
　このように、本実施の形態では、評価値算出部１２は、矩形と異なる形状の規格範囲に基づいて、各候補制御点の評価値を算出してもよい。これにより、特性空間において、円形、楕円形、星形などの規格範囲に基づいて評価値が算出されるため、規格範囲の形状が矩形の場合に限定されることなく、適用場面をさらに拡張することができる。
　［処理動作］
　評価装置１００は、上述のように入力された各データを用いて評価値の算出および出力に関する処理を行う。
　図１０は、本実施の形態に係る評価装置１００の処理動作を示すフローチャートである。
　まず、候補制御点作成部１１は、プロセス条件データ２１１を用いて候補制御点データ２２１を生成する（ステップＳ２１）。
　次に、受信制御部１０は、目的データ２１２を取得する（ステップＳ２２）。つまり、受信制御部１０は、最適化目的を示す目的データ２１２を取得する第２受信ステップを実行する。ここで、未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有する。さらに、受信制御部１０は、制約条件データ２１３を取得する（ステップＳ２３）。つまり、受信制御部１０は、少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データ２１３を取得する第３受信ステップを実行する。さらに、受信制御部１０は、パレート境界を算出する規則を示す領域縮小規則データ２１４を取得する（ステップＳ２４）。つまり、受信制御部１０は、少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、当該分割方法により分割された特性空間の領域ごとに示す領域縮小規則データ２１４を取得する第４受信ステップを実行する。さらに、受信制御部１０は、制御結果データ２２２を記憶部１０５から読み出す（ステップＳ２５）。つまり、受信制御部１０は、第１時刻における制御済みの制御点および第１時刻における既知の特性点を示す制御結果データ２２２を取得する第１受信ステップを実行する。なお、制御結果データ２２４に、何れの特性点も示されていない場合には、このステップＳ２５を含むステップＳ２５～Ｓ２７の処理はスキップされる。
　そして、評価値算出部１２は、目的データ２１２、制約条件データ２１３、領域縮小規則データ２１４、候補制御点データ２２１、および制御結果データ２２２に基づいて、各候補制御点の評価値を算出する（ステップＳ２６）。つまり、評価値算出部１２は、それらのデータに基づいて、複数の未知の特性点それぞれの評価値を算出する算出ステップを実行する。具体的には、評価値算出部１２は、候補制御点データ２２１に示される各候補制御点の評価値を、領域縮小規則に基づいて変更した改善量の算出方法を用いて算出する。また、この算出ステップにおいて、評価値算出部１２は、制約条件の適合度合いに応じた重み付けを、少なくとも１つの製品特性に対する評価値に付与してもよい。そして、評価値算出部１２は、その算出された各候補制御点の評価値を示す評価値データ２２４を生成する。
　次に、評価値出力部１３は、ステップＳ５で算出された評価値、すなわち評価値データ２２４を表示部１０４に出力する（ステップＳ２７）。つまり、評価値出力部１３は、評価値を出力する出力ステップを実行する。これにより、評価値データ２２４が、例えば表示部１０４に表示される。
　そして、受信制御部１０は、ユーザによる入力部１０１ａへの入力操作に応じて、その入力部１０１ａから操作信号を取得する。操作信号は、最適解の探索の終了、または、最適解の探索の続行を示す。なお、最適解の探索は、新規制御結果に基づく各候補制御点の評価値の算出および出力を行う処理である。受信制御部１０は、その操作信号が最適解の探索の終了を示すか、続行を示すかを判定する（ステップＳ２８）。
　ここで、受信制御部１０は、操作信号が最適解の探索の終了を示すと判定すると（ステップＳ２８の「終了」）、全ての処理を終了する。一方、受信制御部１０は、操作信号が最適解の探索の続行を示すと判定すると（ステップＳ２８の「続行」）、次の制御点として選択された候補制御点を、記憶部１０５の制御結果データ２２２に書き込む。量産現場（生産工程）においては、生産が継続している間は続行を選択するとよい。例えば、ユーザが入力部１０１ａに対する入力操作を行うことによって、受信制御部１０は、評価値データ２２４から候補制御点を次の制御点として選択する。受信制御部１０は、このように選択された候補制御点を制御結果データ２２２に書き込む。そして、ユーザは、次の制御点に対応する特性点が制御によって得られると、入力部１０１ａに対する入力操作を行うことによって、その特性点を示す特性点データ２０１を評価装置１００に入力する。受信制御部１０は、その入力された特性点データ２０１を取得し、その特性点データ２０１によって示される特性点を、記憶部１０５の制御結果データ２２２に書き込む。このとき、その特性点は、直近に選択されて書き込まれた制御点に対応付けられる。これにより、新規制御結果が制御結果データ２２２に記録される（ステップＳ２９）。つまり、制御結果データ２２２が更新される。制御結果データ２２２が更新されると、評価値算出部１２は、ステップＳ２５からの処理を繰り返し実行する。
　以上のようなフローを経る過程で、過去の制御結果から、次に行うべき最適な制御条件（すなわち候補制御点）を定量的に解析することができる。その結果、ユーザなどの解析者の力量に依らず開発サイクルの短縮が期待できる。
　＜候補制御点データ２２１の一例＞
　図１１Ａは、候補制御点データ２２１の一例を示す図である。
　候補制御点作成部１１は、プロセス条件データ２１１に基づいて、図１１Ａに示す候補制御点データ２２１を生成する。例えば、候補制御点作成部１１は、プロセス条件データ２１１によって示される全てのプロセス条件のそれぞれの値が連続変数の値であって、その値に関して制約がない場合には、各プロセス条件の値の組み合わせ全通りのそれぞれを候補制御点として作成する。例えば、プロセス条件データ２１１が、第１プロセス条件の連続変数の値「１０，２０，３０，４０，５０」と、第２プロセス条件の連続変数の値「１００，２００，３００，４００，５００」とを示すとする。この場合、候補制御点作成部１１は、第１プロセス条件の値「１０」と第２プロセス条件の値「１００」との組み合わせ、第１プロセス条件の値「１０」と第２プロセス条件の値「２００」との組み合わせなど、全ての組み合わせのそれぞれを候補制御点として作成する。候補制御点作成部１１は、その作成された候補制御点に対して制御点番号を関連付け、その制御点番号が関連付けられた候補制御点を示す候補制御点データ２２１を生成する。
　具体的な一例では、候補制御点データ２２１は、図１１Ａに示すように、制御点番号「１」に関連付けられた候補制御点（１０，１００）、制御点番号「２」に関連付けられた候補制御点（１０，２００）、制御点番号「３」に関連付けられた候補制御点（１０，３００）などを示す。なお、これらの候補制御点の第１成分は、第１プロセス条件の値を示し、第２成分は、第２プロセス条件の値を示す。
　ここで、値の組み合わせ全通りのうち、ある制約を満たす値の組み合わせのみが候補制御点として作成されてもよい。例えば、材料開発において、第１プロセス条件と第２プロセス条件として第１化合物と第２化合物がそれぞれ設定され、値としてそれらの配合比が設定される場合、候補制御点作成部１１は、和が１を満たす値の組み合わせのみを候補制御点として採用する。図１１Ｂの候補制御点データ２２１には、その一例が示されている。
　図１１Ｂは、候補制御点データ２２１の他の例を示す図である。
　候補制御点作成部１１は、プロセス条件データ２１１に基づいて、図１１Ｂに示す候補制御点データ２２１を生成する。例えば、プロセス条件データ２１１が、第２プロセス条件の比率変数の値として「０．０，　０．２，　０．４，　０．６，　０．８，　１．０」を示し、第３プロセス条件の比率変数の値として「０．０，　０．２，０．４，　０．６，　０．８，　１．０」を示すとする。この場合、これらの比率変数の値の組み合わせは、上述の第１化合物と第２化合物との配合比に相当する。したがって、候補制御点作成部１１は、第２プロセス条件の比率変数の値と、第３プロセス条件の比率変数の値との和が１を満たすように、第１プロセス条件の値と、第２プロセス条件の値と、第３プロセス条件の値との組み合わせを、候補制御点として作成すればよい。例えば、候補制御点作成部１１は、第１プロセス条件の値「１０」と、第２プロセス条件の値「０．２」と、第３プロセス条件の値「０．８」との組み合わせなど、比率変数の値の和が１を満たす値の組み合わせを候補制御点として作成すればよい。候補制御点作成部１１は、その作成された候補制御点に対して制御点番号を関連付け、その制御点番号が関連付けられた候補制御点を示す候補制御点データ２２１を生成する。
　具体的な一例では、候補制御点データ２２１は、図１１Ｂに示すように、制御点番号「１」に関連付けられた候補制御点（１０，　０．０，　１．０）、制御点番号「２」に関連付けられた候補制御点（１０，　０．２，　０．８）、制御点番号「３」に関連付けられた候補制御点（１０，　０．４，　０．６）などを示す。なお、これらの候補制御点の第１成分は、第１プロセス条件の値を示し、第２成分は、第２プロセス条件の値を示し、第３成分は、第３プロセス条件の値を示す。
　このように、本実施の形態では、複数のプロセス条件がある場合、候補制御点作成部１１は、複数の候補制御点のそれぞれを作成するときには、その複数のプロセス条件のそれぞれの、所定の条件を満たす値を組み合わせることによって、当該候補制御点を作成する。例えば、所定の条件は、図１１Ｂに示すように、複数のプロセス条件のそれぞれの比率変数の値の和が１であるという条件である。より具体的な一例では、その比率変数は、プロセス条件に対応する化合物などの材料の配合比である。したがって、複数種の化合物の配合比の組み合わせごとに、その組み合わせに対する評価値を算出することができる。その結果、それらの化合物の配合によって得られる合成材料の１以上の製品特性に対する最適解を適切に探索することができる。
　＜制御結果データ２２２の一例＞
　図１２は、制御結果データ２２２の一例を示す図である。
　評価値算出部１２は、評価値を算出するために、記憶部１０５に格納されている制御結果データ２２２を読み出す。この制御結果データ２２２は、図１２に示すように、制御番号ごとに、その制御番号によって識別される制御で用いられた制御点と、その制御によって得られた制御結果である特性点とを示す。制御点は、各プロセス条件の値の組み合わせによって表現される。例えば、制御点は、第１プロセス条件の値「１０」と、第２プロセス条件の値「１００」との組み合わせである値の組み合わせによって表現される。特性点は、制御で得られた各製品特性の値の組み合わせによって表現される。なお、製品特性の値は、以下、製品特性値とも呼ばれる。例えば、特性点は、第１製品特性の値「８」と、第２製品特性の値「０．０」との組み合わせによって表現される。
　具体的な一例では、制御結果データ２２２は、図１２に示すように、制御番号「１」に関連付けられた制御点（１０，１００）および特性点（８，　０．０）、制御番号「２」に関連付けられた制御点（１０，５００）および特性点（４０，　１．６）、制御番号「３」に関連付けられた制御点（５０，１００）および特性点（４０，　１．６）などを示す。
　［評価値の算出処理］
　図１３は、評価値算出部１２による処理を説明するための図である。評価値算出部１２は、候補制御点作成部１１によって生成された候補制御点データ２２１と、記憶部１０５にある制御結果データ２２２とに基づいて、予測分布データ２２３を生成する。そして、評価値算出部１２は、各製品特性の最適化目的を示す目的データ２１２と、各製品特性の規格範囲を示す制約条件データ２１３と、パレート境界を算出する規則を示す領域縮小規則データ２１４と、予測分布データ２２３とに基づいて、評価値データ２２４を生成する。
　ここで、制御結果データ２２２は、複数の候補制御点のうちの、既に制御に用いられた１以上の候補制御点である１以上の制御点と、その１以上の制御点のそれぞれに対応する特性点であって、その制御点を用いた１以上の製品特性の制御結果とを示す。したがって、本実施の形態における評価値算出部１２は、（ａ）１以上の製品特性のそれぞれの最適化目的および規格範囲と、（ｂ）複数の候補制御点のうちの、既に制御に用いられた１以上の候補制御点である１以上の制御点と、（ｃ）１以上の制御点のそれぞれに対応する特性点であって、当該制御点を用いた１以上の製品特性の制御結果を示す特性点とに基づいて、各候補制御点の評価値をベイズ最適化に基づいて算出する。
　評価値算出部１２は、その生成された評価値データ２２４を評価値出力部１３に出力する。なお、評価値算出部１２は、予測分布データ２２３も評価値出力部１３に出力してもよい。または、評価値算出部１２は、予測分布データ２２３を記憶部１０５に格納し、評価値出力部１３は、ユーザによる入力部１０１ａへの入力操作に応じて、その記憶部１０５から予測分布データ２２３を読み出してもよい。
　評価値算出部１２は、候補制御点と特性点との対応関係を上述したカルマンフィルタを用いて記述する。
　評価値算出部１２は、上記ステップＳ２５において記憶部１０５から読み出された制御結果データ２２２に示される既知の制御結果に対して、上記（式３）を用いた演算を行うことにより、予測分布データ２２３を生成する。
　図１４は、予測分布データ２２３の一例を示す図である。予測分布データ２２３は、各候補制御点における予測分布の平均と分散とを示す。この予測分布は、各製品特性についてカルマンフィルタに基づく（式３）によって算出された分布である。例えば、予測分布データ２２３は、図１４に示すように、制御点番号ごとに、その制御点番号に対応する、第１製品特性の予測分布の平均および分散と、第２製品特性の予測分布の平均および分散とを示す。
　具体的な一例では、予測分布データ２２３は、図１４に示すように、制御点番号「１」に対応する第１製品特性の平均「２３．５３２２」および分散「１９．４０１２」と、第２製品特性の平均「０．７７６６１」および分散「０．９７００６」とを示す。さらに、予測分布データ２２３は、制御点番号「２」に対応する第１製品特性の平均「３０．２５３６」および分散「２１．５５２１」と、第２製品特性の平均「１．１１２６８」および分散「１．０７７６１」とを示す。なお、制御点番号は、図１１Ａまたは図１１Ｂに示すように、候補制御点に対応付けられている。
　評価値算出部１２は、ベイズ最適化における獲得関数と呼ばれる評価基準に基づき評価値を算出する。この評価値の算出には、上述の予測分布が用いられる。また、本実施の形態における獲得関数は、制約条件があるベイズ最適化における獲得関数である。
　＜制約条件のないベイズ最適化の獲得関数＞
　まず、本実施の形態における獲得関数の説明の前に、制約条件のないベイズ最適化の獲得関数（すなわち非特許文献１のＥＨＶＩ）について説明する。ただし、最大化と最小化については、それらのうちの一方の符号を反転させると、他方と等価になるため、最小化に統一して説明する。ＥＨＶＩでは、改善領域の体積（改善量とも呼ばれる）が大きいほど、暫定の制御結果から大きく改善された特性点が得られたと考える。その改善領域は、行った制御から既に得られている少なくとも１つの特性点の中のパレート点（すなわち非劣解）の座標から定まるパレート境界と、新規特性点が観測された際に新規特性点により新たに定まるパレート境界とで囲まれた領域である。なお、パレート点は、現時点で暫定的にパレート解とされている特性点である。例えば、第１製品特性および第２製品特性のそれぞれの最適化目的が最小化である場合には、パレート点よりも第１製品特性および第２製品特性のいずれの値も小さい他の特性点は存在しない。パレート境界は、パレート点の座標のそれぞれを、第１製品特性および第２製品特性の方向に沿って結ぶことにより定まる境界線である。また、以降では、パレート境界で特性空間全体が区分けされたうち、各製品特性に関して値が小さい側をアクティブ領域、値が大きい側を非アクティブ領域と呼ぶ。新規特性点が非アクティブ領域に入ったときの改善量は０とする。
　図１５Ａは、改善領域の一例を示す図である。
　例えば、図１５Ａに示すように、４つのパレート点２１～２４から定まるパレート境界３１と、１つの新規特性点ｙnewが得られた際に新たに定まるパレート境界３２とで囲まれた領域が、改善領域として特定される。
　ここで、カルマンフィルタにより、各候補制御点を選択した場合の各製品特性値の振る舞いは正規分布の形で表現されており、観測された特性点の位置によって改善量も変動する。ＥＨＶＩは、以下の（式６）のように、各候補制御点について、予測分布で改善量の期待値を取った量として定義される。ＥＨＶＩによって得られる値が大きい候補制御点ほど改善量の期待値が大きく、次に実行すべき制御点を表している。
　（式６）において、Ｄは製品特性の数（すなわち次元数）を表し、
はＤ次元ユークリッド空間を表し、Ｉ（ｙnew）は改善量を表す。また、ｐ（ｙnew｜ｘnew）は、少なくとも１つの候補制御点の中から１つの候補制御点を新規制御点ｘnewとして選択した際の、その新規制御点ｘnewに対応する特性点ｙnewの予測分布を表す。その特性点ｙnewの各次元の予測分布、すなわち平均および分散は、上記（式３）により求まっている。
　＜制約条件があるベイズ最適化の獲得関数＞
　次に、本実施の形態における獲得関数について説明する。本実施の形態における獲得関数は、制約条件がある場合のベイズ最適化の獲得関数である。なお、制約条件がない場合には、評価値算出部１２は、上記（式６）で表される獲得関数を用いればよい。Ｄ個の製品特性のうち、ｙ１～ｙDminimizeのＤDminimize個の製品特性の最適化目的が最小化であり、残りのｙDminimize+1～ｙDのＤrange（＝Ｄ-Ｄminimize）個の製品特性の最適化目的が規格範囲内であるとする。このとき、本実施の形態における獲得関数、すなわち制約条件付きＥＨＶＩＣは、以下の（式７）のように定義される。
　（式７）において、Ｒminimizeは、最適化目的が最小化である製品特性ｙ1～ｙDminimizeについて、すべて規格範囲内である領域を表す。Ｒrangeは、最適化目的が規格範囲内である製品特性ｙDminimize+1～ｙDについて、すべて規格範囲内である領域を表す。なお、ＲminimizeおよびＲrangeのそれぞれの領域は、その領域に対応する規格範囲の形状を示す関数によって表現される。図９Ｂに示すように規格範囲の形状が円であれば、ＲminimizeおよびＲrangeのそれぞれの領域は、その円を示す関数によって表現される。また、規格範囲の形状が星形であれば、ＲminimizeおよびＲrangeのそれぞれの領域は、その星形を示す関数によって表現される。ｙnew，minimizeは、特性点ｙnewの全ての次元から、最適化目的が最小化である製品特性の各次元を抽出することによって得られるベクトルを表す。ｙnew，rangeは、特性点ｙnewの全ての次元から、最適化目的が規格範囲内である製品特性の各次元を抽出することによって得られるベクトルを表す。ＩＣ（ｙnew）は、制約条件がある場合の改善量であり、既存のパレート境界と、新たに定まるパレート境界とで囲まれた領域の体積を表す。その既存のパレート境界は、規格範囲内に存在する少なくとも１つのパレート点およびその規格範囲のそれぞれの座標から定まる境界である。新たに定まるパレート境界は、新規特性点が観測された際に、その新規特性点であるパレート点および規格範囲のそれぞれの座標から定まる境界である。Pr{A}は、事象Ａが成立する確率を表し、例えば（式３）で算出される平均および分散を用いて表現される。
　図１５Ｂは、本実施の形態に係る改善領域の他の例を示す図である。本実施の形態と非特許文献２との大きな違いは、最適化目的が最小化である製品特性について、本実施の形態では積分範囲が特性空間全体から規格範囲内に制限されていて、改善量の測り方が規格範囲に応じて変わっていることである。規格範囲における最大値および最小値が指定されない場合は、最大値は＋∞として設定され、最小値は－∞として設定される。最適化目的が最小化であるすべての製品特性の規格範囲における最大値および最小値がそれぞれ＋∞および－∞で、かつ、Ｄrange＝０のとき、本実施の形態における獲得関数であるＥＨＶＩＣは、非特許文献１のＥＨＶＩに帰着される。また、最適化目的が最小化であるすべての製品特性の規格範囲における最大値および最小値がそれぞれ＋∞および－∞で、かつ、Ｄrange＞＝１のとき、本実施の形態における獲得関数であるＥＨＶＩＣは、非特許文献２のＥＨＶＩＣに帰着される。したがって、本実施の形態における評価装置１００は、従来手法でも評価値を算出することができる。
　また、非特許文献２では、最適化目的が最小化である製品特性が一つ以上存在する最適化問題、つまり、Ｄminimize＞＝１を想定しているが、本実施の形態における獲得関数では、Ｄminimize＝０（すなわちＤrange＝Ｄ）の場合でも不都合なく定式化が可能である。したがって、本実施の形態における獲得関数は、すべての製品特性の最適化目的が規格範囲内である最適化問題にも自然に拡張される。
　次に、本実施の形態における獲得関数であるＥＨＶＩＣの具体的な算出方法に関して説明する。
　図１６Ａは、改善領域の体積の算出方法を説明するための図である。なお、図１６Ａの（ａ）は、特性空間における改善領域を示し、図１６Ａの（ｂ）は、分割の対象とされるその改善領域を示し、図１６Ａの（ｃ）は、改善領域の分割によって得られる複数の小領域を示す。
　評価値算出部１２は、最適化目的が最小化である製品特性の次元については、図１６Ａのように、改善領域の体積である改善量（すなわちＩＣ（ｙnew））を算出する。つまり、評価値算出部１２は、パレート点および新規特性点のそれぞれの座標で改善領域を複数の小領域に分割し、各小領域の体積の期待値を算出したのち、それらの期待値の和を取ることで改善量（すなわちＩＣ（ｙnew））を算出する。また、評価値算出部１２は、最適化目的が規格範囲内である製品特性の次元については、各製品特性値が規格範囲に入る確率を算出する。
　図１６Ｂは、特性空間の全体を複数の小領域に分割する例を示す図である。
　評価値算出部１２は、最適化目的が最小化である製品特性の次元と、最適化目的が規格範囲内である製品特性の次元とについて、図１６Ｂのように特性空間全体を複数の小領域に分割し、以下の（式８）を用いることによって、獲得関数を統一的に算出する。つまり、評価値算出部１２は、パレート点、新規特性点、および規格値のそれぞれの座標で特性空間全体を複数の小領域に分割し、各小領域の体積の算出を以下の（式８）のような場合分け計算で実行する。なお、上述の規格値は、規格範囲の最大値および最小値である。そして、評価値算出部１２は、それらの小領域の体積を期待値処理したものの和を取ることで、制約条件がある場合の獲得関数を統一的に算出する。なお、その体積は、Ｄ次元超体積とも呼ばれる。
　（式８）において、ｙdは、小領域の下端点(y1,…,yD)の第d成分を表し、y’dは、小領域の上端点(y’1,…,y’D)の第d成分を表す。
　図１６Ｃは、小領域の下端点および上端点の一例を示す図である。
　Ｄ＝２の場合、図１６Ｃに示すように、（y1,y2）は、小領域の下端点を表し、（y’1,y’2)は、小領域の上端点を表す。
　また、（式８）における（ｉ）は、区間[yd,y’d]が次元ｄに関して規格範囲外であるときに適用される。（ｉｉ）は、区間[yd,y’d]が次元ｄに関して規格範囲内であり、かつ、次元ｄの製品特性の最適化目的が規格範囲内であるときに適用される。（ｉｉｉ）は、区間[yd,y’d]が次元ｄに関して規格範囲内であり、かつ、次元ｄの製品特性の最適化目的が最小化であるときに適用される。ｃdは、重み係数であり、製品特性の次元ｄごとに探索の優先度を付与する際等に適宜設定される。例えば、次元ｄの優先度が高いほど小さい重み係数ｃdが用いられ、逆に、次元ｄの優先度が低いほど大きい重み係数ｃdが用いられる。重み係数ｃｄの逆数が優先度であってもよい。特に指定がなければ、すなわち、各次元ｄの優先度が等しい場合には、各次元ｄのｃdは、例えばすべて１に設定される。
　以上が本実施の形態における獲得関数と、その獲得関数の具体的な算出方法の説明である。
　これにより、制約条件のある多目的最適化問題においても、一貫した手続きでの定量的な評価が可能になる。より具体的には、制約条件がある場合のベイズ最適化の獲得関数を利用して、各候補制御点について、特性空間における制約範囲内の体積を、最適化の改善量として算出し、その改善量から評価値を適切に算出することできる。
　しかしながら、制約条件がある場合のベイズ最適化の獲得関数を算出する際、例えば図１６Ａのように改善領域を複数の小領域に分割し、各小領域の体積の期待値を算出したのち、それらの期待値の和を取る必要がある。より一般化すると、製品特性の数すなわち次元数をＤとしたとき、改善領域は、複数のＤ次元超直方体の和領域で表すことができる。このため、制約条件がある場合のベイズ最適化の獲得関数を算出する際、改善領域を複数のＤ次元超直方体に分割し、各超直方体の体積の期待値を算出したのち、それらの期待値の和を取ることで算出する必要がある。したがって、獲得関数の計算量は、改善領域を構成する超直方体の個数に大きく依存することになる。超直方体の個数は、製品特性の数（次元数）であるＤと、観測された特性点のうちの中のパレート点の数であるパレート点数Nparetoを用いて、
オーダーで指数関数的に増加する。なお、実用的には、概ね3次元（Ｄ＝３）が計算できる限界となる。
　そこで、以下では、探索効率を保ったまま計算量を削減し、定量的な評価を高速に実行するために改良した改善領域の体積の算出方法について説明する。
　＜改良した改善領域の体積の算出方法＞
　改良した改善領域の体積の算出方法について説明する前に、パレート境界を算出する規則を示し改善量の算出方法を変更させる領域縮小規則について説明する。
　領域縮小規則では、特性空間を所定数の領域に分割する方法と、パレート境界の算出方法とが示される。以下では、説明を簡単にするため、規格範囲が設定されていないとして説明する。
　領域縮小規則が適用される場合、製品特性の数Ｄ個すなわちＤ次元の製品特性に対して、特性空間全体をＤ＋１個の領域に分割（領域分割）される。領域分割する方法は、任意の方法でよい。領域分割された特性空間の領域は、順に領域１，領域２，・・・，領域Ｄ，領域Ｄ＋１と名付けられるとする。
　図１７は、領域縮小規則が適用される場合に領域分割される特性空間の一例を示す図である。図１７に示される例では、製品特性の数が２個（Ｄ＝２）であることから、特性空間が３個の領域すなわち領域１、領域２及び領域３に領域分割されている。
　図１７に示される例では、特性空間において「－∞」～「１０」の第１製品特性の範囲と「－∞」～「１０」の第２製品特性の範囲とからなる第３領域が示されている。また、特性空間において第３領域を除く領域で４５度の傾きで規定される直線よりも下の領域である第１領域と、第３領域を除く領域で４５度の傾きで規定される直線よりも上の領域である第２領域とが示されている。
　なお、領域縮小規則が適用される場合、領域分割された特性空間のＤ＋１個の領域には、空集合が設定されてもよいが、Ｄ＋１個の領域すべて空集合はナンセンスであるため少なくとも１つは空でない集合とする。
　図１８Ａ～図１８Ｃは、本実施の形態に係る領域縮小規則が適用される場合に領域分割される特性空間の別の例を示す図である。
　図１８Ａに示される例では、２次元の製品特性に対して、領域３が空集合に設定されることで特性空間が２個の領域に分割される場合が示されている。より具体的には、図１８Ａに示されるように、第３領域が空集合であることから、特性空間において４５度の傾きで規定される直線よりも下の領域である第１領域と、特性空間において４５度の傾きで規定される直線よりも上の領域である第２領域とに領域分割されている。
　図１８Ｂに示される例では、２次元の製品特性に対して、領域２及び領域３が空集合に設定されることで特性空間が１個の領域１に分割される場合が示されている。より具体的には、図１８Ｂに示されるように、第２領域及び第３領域が空集合であることから、特性空間全体が第１領域のみに領域分割されている。
　図１８Ｃに示される例では、２次元の製品特性に対して、領域３が空集合に設定されることで特性空間が２個の領域に分割される場合が示されている。つまり、図１８Ｃに示される例では、第３領域が空集合であることから、特性空間全体が第１領域及び第２領域に領域分割されている。より具体的には、図１８Ｃに示されるように、特性空間において、円の中心（１０，１０）と半径５とによって表現される領域である第１領域と、特性空間において第１領域以外の領域である第２領域とに領域分割されている。
　このように、領域縮小規則が適用される場合、領域分割された特性空間のＤ＋１個の領域には、空集合が設定されてもよい。ただし、この場合、特性空間上の任意の点は、空集合以外のいずれか一つの領域に振り分けられる。
　続いて、領域縮小規則が適用される場合におけるパレート境界の算出方法について説明する。
　ここで、パレート点は、上述したように、現時点で暫定的にパレート解とされている特性点であり、非劣解とも称される。例えば２次元の製品特性から構成される特性空間において、第１製品特性および第２製品特性のそれぞれの最適化目的が最小化であるとする。この場合、パレート点は、観測されたすべての他の特性点と比較して、その点よりも第１製品特性および第２製品特性のいずれの値も小さい他の特性点は存在しない特性点となる。
　また、パレート境界は、少なくとも１つのパレート点の座標から定まる境界である。例えば図１６Ａに示されるように、上述したパレート境界は、パレート点の座標のそれぞれを、第１製品特性および第２製品特性の値が大きい方向へ延長して結ぶことにより定まる境界線であった。一方、領域縮小規則が適用されるとパレート境界の算出方法が変更される。すなわち、領域縮小規則が適用される場合、アクティブ領域を縮小する次元が領域ごとに定められることでパレート境界の定義が変更される。
　ここで、領域縮小規則が適用されてパレート境界の定義が変更された場合におけるパレート境界の算出方法について説明する。
　図１９～図２２は、領域縮小規則が適用されるパレート境界の算出方法を説明するための図である。図１９～図２２の（ａ）には、定義変更前のパレート境界の算出方法すなわち領域縮小規則が適用されないパレート境界の算出方法の例が示されている。図１９～図２２の（ｂ）には、定義変更後のパレート境界の算出方法すなわち領域縮小規則が適用されるパレート境界の算出方法の例が示されている。なお、図１９～図２２では、例えば２次元の製品特性から構成される特性空間において、第１製品特性および第２製品特性のそれぞれの最適化目的が最小化であるとしている。また、図１９～図２２の（ｂ）では、特性空間において原点を通り４５度の傾きで規定される直線で領域１と領域２とに領域分割されているとしている。
　例えば図１９において、１つ目の新規特性点ｙnew(1)が得られた際、新規特性点ｙnew(1)がパレート点となる。この場合、図１９の（ａ）に示すように、領域縮小規則が適用されず定義変更されないときには、新規特性点ｙnew(1)の座標それぞれを、第１製品特性および第２製品特性の方向に沿って結ぶことにより定まる境界線がパレート境界となる。一方、図１９の（ｂ）に示すように、領域縮小規則が適用されて定義変更されるときには、新規特性点ｙnew(1)は、領域２に位置するため、新規特性点ｙnew(1)の第１製品特性の座標を通り、第２製品特性の軸に平行により定まる境界線がパレート境界となる。換言すると、新規特性点ｙnew(1)の第１製品特性の座標ｙnew1(1)は、新規特性点ｙnew(1)の第２製品特性の座標ｙnew2(1)よりも大きい。このため、新規特性点ｙnew1(1)より右側の領域を、非アクティブ領域とする。また、このことは、新規特性点ｙnew(1)の座標ｙnew1(1)でアクティブ領域を縮小（削減）すると表現することもできる。
　次に、例えば図２０において、２つ目の新規特性点ｙnew(2)が得られた際、新規特性点ｙnew(2)はパレート点となる。この場合、図２０の（ａ）に示すように、領域縮小規則が適用されず定義変更されないときには、特性点ｙnew(1)の座標と新規特性点ｙnew(2)の座標とを、第１製品特性および第２製品特性の値が大きい方向へ延長して結ぶことにより定まる境界線がパレート境界となる。一方、図２０の（ｂ）に示すように、領域縮小規則が適用され定義変更されるときには、新規特性点ｙnew(2)は、領域１に位置する。このため、新規特性点ｙnew(2)の第２製品特性の座標を通り、第１製品特性の軸に平行により定まる線と特性点ｙnew(1)の第１製品特性の座標を通り、第２製品特性の軸に平行により定まる線とからなる境界線がパレート境界となる。換言すると、新規特性点ｙnew(2)の第２製品特性の座標ｙnew2(2)は、新規特性点ｙnew(2)の第１製品特性の座標ｙnew1(2)よりも大きい。このため、新規特性点ｙnew1(1)より右側の領域、または、新規特性点ｙnew2(2)より上側の領域を、非アクティブ領域とする。また、このことは、新規特性点ｙnew(2)の座標ｙnew2(2)でアクティブ領域をさらに縮小（削減）していると表現することもできる。
　なお、例えば図２１において、３つ目の新規特性点ｙnew(3)が得られた際、新規特性点ｙnew(3)はパレート点とならない。この場合、図２１の（ａ）および（ｂ）に示すように、パレート境界は変更されないことになる。換言すると、新規特性点ｙnew(3)がアクティブ領域に含まれず、非アクティブ領域に含まれる場合には、パレート点とならないため、パレート境界は変更されない。
　次に、例えば図２２の（ａ）において、４つ目の新規特性点ｙnew(4)が得られた際、新規特性点ｙnew(4)はパレート点となる。この場合、図２２の（ａ）では、新規特性点ｙnew(4)がアクティブ領域に含まれるため、パレート境界は変更される。一方、図２２の（ｂ）では、新規特性点ｙnew(4)が非アクティブ領域に含まれるため、パレート境界は変更されないことになる。
　このように、領域縮小規則が適用される場合、アクティブ領域を縮小する次元が領域ごとに定められることでパレート境界の定義が変更されることになる。
　図２３は、制約条件がない場合のパレート境界の一例を示す図である。図２３には、図１７に示されるように領域分割された場合に算出されたパレート境界の一例が示されている。なお、図２３では、例えば２次元の製品特性から構成される特性空間において第１製品特性および第２製品特性のそれぞれの最適化目的は最小化であるとしている。また、図２３に示す例では、領域１および領域２は、パレート境界の算出に用いられる一方で、領域３ではパレート境界の算出には用いられないとしている。
　そして、このようなパレート境界は、（式９）のように定式化できる。
　ここで、Ｄ次元の製品特性に対して、特性空間全体はＤ＋１個の領域に分割されており、各ｄ＝１，…Ｄに対して、観測された特性点のうち、領域ｄに含まれる特性点の中で、ｙd座標が一番小さい座標をy’dとおく。なお、領域Ｄ＋１に含まれる特性点に対しては何もしない。また、y’dの初期値は各規格上限値とする。このとき、（式９）で表される領域が、領域縮小規則が適用される場合のパレート境界となる。（式９）において、Ｄは製品特性の数（次元数）を表し、
は、ｙがＤ次元ユークリッド空間の要素であることを表す。＼は、バックスラッシュの左側の集合からバックスラッシュの右側の集合に含まれる要素を取り除いた集合（差集合）であることを表す。
すなわちターンエーは、集合における「任意」の要素を取ることを表す。
　したがって、評価値算出部１２は、（式９）により、領域ｄに含まれる特性点の中で、ｙd座標が一番小さい座標y’dを有する特性点における座標y’dで定められる境界をパレート境界として算出することになる。
　図２４は、図２３に示されるパレート境界のもとでの改善領域の一例を示す図である。
　このように、評価値算出部１２は、最適化目的が最小化である製品特性の次元と、最適化目的が規格範囲内である製品特性の次元とについて、図２４のように、改善領域の体積である改善量（すなわちＩ（ｙnew））を算出することができる。つまり、領域縮小規則が適用される場合、評価値算出部１２は、既存のパレート境界と、新たに定まるパレート点とで定まる１つの小領域の体積の期待値を算出することで改善量（すなわちＩ（ｙnew））を算出することができる。これを直観的に説明すると、評価値算出部１２は、１つの小領域の体積の期待値で表すことのできる非アクティブ領域の増加量により、改善量を算出することができる。
　このような改善量を算出する方法は、新規特性点ｙnewのｙd座標をｙnew,dとおくと、ｙnewが観測された際の改善領域の体積（改善量と呼ばれる）として（式１０）のように定義できる。
　（式１０）において、ｙnew,dは、新規特性点の次元ｄの座標（第ｄ成分）を表し、ｙ’dは、パレート境界を定める観測された特性点（パレート点）の次元ｄの座標（第ｄ成分）を表す。また、（式１０）において、ある次元ｄに対して、ｙnew,dがｙ’d未満であれば、改善量（すなわちＩ（ｙnew））は非負実数となる一方で、ある次元ｄに対して、ｙnew,dがｙ’d以上なら改善量（すなわちＩ（ｙnew））は、負実数ではなく０となる。
　そして、このように算出される改善量の体積で表される改善量（すなわちＩ（ｙnew））が大きいほど、暫定の実験結果から大きく改善された特性点が得られたと考えることができる。
　また、本実施の形態では、領域縮小規則が適用される場合における獲得関数（すなわち制約条件付きＥＨＶＩＣ）は、ＥＨＶＩと同様に、各候補制御点について、カルマンフィルタに基づき算出された予測分布を用いて定義される。より具体的には、領域縮小規則が適用される場合における獲得関数は、（式７）のように、改善量の期待値を取った量で定義できる。そして、改善量の期待値を取った量の値の大小で、次に実行すべき制御点の良し悪しを評価する。
　以上のように、領域縮小規則が適用される場合における獲得関数の算出方法によれば、領域縮小規則が適用されない場合のような小領域への分割および和計算を必要とせずに、単一のＤ次元超直方体の体積の期待値さえ算出すれば、獲得関数を算出することができる。これにより、領域縮小規則が適用される場合における獲得関数の計算量は、パレート点数Nparetoには非依存で、製品特性数Dの増加に関して多項式オーダーでの増加に抑えられるので、探索効率を保ったまま高速な解析処理を実現できる。
　なお、上記では、規格範囲が設定されていない場合に適用される領域縮小規則について説明したがこれに限らない。規格範囲が設定されていてもよく、同様に領域縮小規則が適用される。
　図２５は、制約条件がある場合のパレート境界の一例を示す図である。図２６は、図２５に示されるパレート境界のもと制約条件がある場合に定められる改善領域の一例を示す図である。
　図２６は、図２３と比較して規格範囲が設定されている点が異なり、その他は同じである。すなわち、図２５には、図１７に示されように領域分割され、かつ、２次元の製品特性から構成される特性空間において規格範囲が設定されているときに算出されたパレート境界の一例が示されている。なお、図２５に示される例でも、第１製品特性および第２製品特性のそれぞれの最適化目的は最小化であるとしている。
　このような場合でも、規格範囲内において領域縮小規則が定められるので、新規特性点が観測されるたびに、上述したのと同様の手続きでアクティブ領域の縮小が実行される。すなわち、評価値算出部１２は、（式９）により、領域ｄに含まれる特性点の中で、ｙd座標が一番小さい座標ｙ’dを有し、かつ規格範囲内にある特性点の座標ｙ’dで定められる境界をアクティブ境界として算出する。
　すると、改善領域が図２６のように定まるので、評価値算出部１２は、（式７）および（式８）を用いて、獲得関数の評価を算出することができる。
　以上の獲得関数の算出方法は、厳密解を求めるための方法であり、予測分布がカルマンフィルタなどを用いて正規分布で算出でき、かつ改善領域が単一の矩形のように、獲得関数が解析的に算出可能であれば、上記の方法で算出可能である。しかし、予測分布が正規分布でない、または改善領域が複雑な形状である場合には、獲得関数を解析的に算出できない可能性がある。そのような場合には、モンテカルロ法等を用いて獲得関数を近似的に算出してもよい。その場合でも、特性空間の小領域への分割および改善領域等は、上述の説明と変わらない。
　なお、上記では、制約条件として規格範囲が設けられている場合についても説明したが、これに限らない。規格範囲だけでなく、その規格範囲とは別に範囲が設けられてもよい。例えば、特性点が最低限収まってほしい規格範囲の中に、特性点が可能な限り収まってほしい管理範囲が設定されているケースも、実務においてしばしば要請される。なお、規格範囲および管理範囲は、それぞれ制約条件である制約範囲の一例である。
　図２７は、規格範囲および管理範囲の一例を示す図である。
　図２７に示す例では、第１製品特性の規格範囲は、最小値「１０」～最大値「５０」であり、第２製品特性の規格範囲は、最小値「１０」～最大値「５０」である。また、第１製品特性の管理範囲は、その規格範囲よりも狭い範囲、すなわち最小値「２０」～最大値「４０」であり、第２製品特性の管理範囲は、その規格範囲よりも狭い範囲、すなわち最小値「２０」～最大値「４０」である。このように、管理範囲は規格範囲に含まれてもよい。
　このような場合、評価値算出部１２は、（式８）において、０と１との間にある例えば０．５等の中間的な値を、ｌ（ｙｄ，ｙｄ’）に対してさらに設定することで評価値を算出することができる。なお、０．５は一例であって、他の数値であってもよい。
　［出力］
　評価値出力部１３は、評価値算出部１２によって上述のように算出された各候補制御点の評価値を示す評価値データ２２４を取得し、その評価値データ２２４を表示部１０４に表示させる。なお、評価値出力部１３は、評価値算出部１２からその評価値データ２２４を直接取得してもよく、評価値算出部１２によって記憶部１０５に格納された評価値データ２２４を読み出すことによって、その評価値データ２２４を取得してもよい。
　図２８は、評価値データ２２４の一例を示す図である。評価値データ２２４は、例えば図２８に示すように、各候補制御点における評価値およびその順位を示す。具体的には、評価値データ２２４は、制御点番号ごとに、その制御点番号に対応する評価値と、その評価値の順位とを示す。各制御点番号は、図１１Ａおよび図１１Ｂに示すように、候補制御点に対応付けられている。したがって、評価値データ２２４は、候補制御点ごとに、その候補制御点に対応する評価値と、その評価値の順位とを示していると言える。また、順位は、評価値が大きいほど小さい数値を示し、逆に、評価値が小さいほど大きい数値を示す。
　具体的な一例では、評価値データ２２４は、図２８に示すように、制御点番号「１」に対応する評価値「０．０００００」および順位「２３」、制御点番号「２」に対応する評価値「０．８７６８２」および順位「１」、制御点番号「３」に対応する評価値「０．６２３４２」および順位「４」などを示す。
　図２９は、表示部１０４に表示される変更後の評価値データ２２４の一例を示す図である。
　評価値出力部１３は、評価値データ２２４を評価値順位でソートするなどにより、その評価値データ２２４を変更し、その変更後の評価値データ２２４を表示部１０４に表示してもよい。
　変更後の評価値データ２２４は、例えば図２９に示すように、評価値の順位ごとに、その順位に該当する評価値と、その評価値に対応する候補制御点とを示す。また、評価値の各順位は昇順に配列されている。つまり、評価値の大きい順に各候補制御点が配列されている。例えば、評価値データ２２４は、順位「１」に該当する評価値「０．８７６８２」と、その評価値に対応する候補制御点（１０，２００）とを示す。さらに、評価値データ２２４は、順位「２」に該当する評価値「０．８７６８２」と、その評価値に対応する候補制御点（２０，１００）とを示す。
　このような評価値データ２２４が表示部１０４に表示されることによって、ユーザは、最適解の探索を続行するか終了するかを判断することができる。さらに、ユーザは、最適解の探索を続行する場合には、表示されている各評価値および各順位に基づいて、表示されている全ての制御点番号、すなわち全ての候補制御点から、次の制御点とされる候補制御点を選択することができる。例えば、ユーザは、最も大きい評価値（すなわち、順位が１である評価値）に対応する候補制御点を選択する。このとき、ユーザは、入力部１０１ａに対する入力操作を行うことによって、評価値データ２２４の各評価値を大きい順に並び替えさせてもよい。つまり、評価値出力部１３は、評価値データ２２４の各評価値が降順となり、各順位が昇順となるように、それらをソートする。これにより、最も大きい評価値を見つけやすくすることができる。
　［時系列データの制御への応用例］
　次に、時系列データの制御への応用例について説明する。
　量産現場では、例えば量産ラインの稼働再開時に製品特性を早期に安定させるため、リアルタイムで最適なプロセス条件の設定値（制御点）を探索したいというニーズがある。これは、時系列性を考慮した最適解の探索問題に該当し、各時刻で過去の制御情報（プロセス条件設定値などで定まる制御点についての情報）をもとに、次時点の最適な対応関係（制御点と特性点）を獲得する問題と解することができる。
　そこで、本実施の形態では、カルマンフィルタで求まる予測分布とベイズ最適化で求まる評価関数とを用いて、あらゆる関係性（対応関係）を確率分布で考慮しながら評価することができる評価値を用いて、最適な対応関係を獲得する。
　図３０は、カルマンフィルタによる予測分布の導出を概念的に示す図である。
　図３０に示すように、カルマンフィルタによる状態推定は、概念的には、毎時刻(t-4、t-3、t-2、t-1、t)において、１つ前の時刻の観測値から状態の予測分布の平均と分散とを逐次推定していることに該当する。この予測分布の平均と分散とに基づき、推定値を得ることができる。なお、状態の予測分布の平均と分散とを逐次推定することは、１つ前の時刻の観測値で補正されていると解釈することもできる。
　換言すると、本実施の形態では、毎時刻で制御結果すなわち制御に用いられた１以上の制御点と当該１以上の制御点のそれぞれに対応する特性点とが得られるごとに、カルマンフィルタで予測分布（の平均と分散）を求める。そして、ベイズ最適化で求まる評価関数に基づき、次に設定する制御点を選択する。
　図３１は、複数の候補制御点の予測分布と製品特性の規格範囲との関係を概念的に示す図である。図３１の（ａ）には候補制御点Ａの予測分布、図３１の（ｂ）には候補制御点Ｂの予測分布、図３１の（ｃ）には候補制御点Ｃの予測分布が概念的に示されている。
　図３１に示すように、候補制御点Ａ、候補制御点Ｂおよび候補制御点Ｃの予測分布（の平均および分散）が同じ場合、製品特性を早期に安定させるという観点で見ると予測分布の平均は規格中央に近い方がよい。したがって、図３１に示される候補制御点Ａ、候補制御点Ｂおよび候補制御点Ｃのうちでは、候補制御点Ｂが次に設定する制御点として選択されるとよい。
　つまり、ベイズ最適化において、最適化目的が規格範囲内である場合、規格中央を最適化目標値に設定し、直前の特性点（観測値）の位置が規格中央より上側なら最適化方法（最適化目的）を最小化、直前の特性点の位置が規格中央より下側なら最適化目的（最適化目的）を最大化に切り替える制御とほぼ等価となる。
　＜時系列データの制御方法１＞
　図３２Ａおよび図３２Ｂは、本実施の形態に係る最適化目標値を規格中央に設定する制御方法１の制御イメージを説明するための図である。なお、図３２Ａおよび図３２Ｂでは、説明を簡潔にするため製品特性数が１次元である場合の制御イメージが示されている。
　最適化目的が規格範囲内である場合、規格中央を最適化目標値に設定し、かつ、直前の特性点（観測値）の位置が規格中央より上側なら最適化方法（最適化目的）を最小化、直前の特性点の位置が規格中央より下側なら最適化目的（最適化目的）を最大化に切り替える制御であると解釈できる。このような制御方法を、以下では時系列データの制御方法１と称する。
　つまり、図３２Ａおよび図３２Ｂに示されるように、時系列データの制御方法１では、最適化目的が規格範囲内である場合、規格中央が最適化目標値に設定される。具体的には、図３２Ａに示す例では、直前の特性点ｙｔが規格中央よりも上側であるので、最適化方法（最適化目的）を最小化にして制御を行えばよい。一方、図３２Ｂに示す例では、直前の特性点ｙｔ＋１が規格中央よりも下側であるので、最適化方法（最適化目的）を最大化にして制御を行えばよい。
　このように時系列データの制御方法１では、最適化目的が規格範囲内であることを、最適化目標値が規格中央であり、最適化方法が直前の特性点（観測点）の位置に応じて最大化または最小化に切り替わる制御方法であると解釈する。
　図３３は、オーバーシュート現象およびハンチング現象を概念的に示す図である。図３４は、カルマンフィルタの特性であるノイズキャンセリング効果を概念的に示す図である。
　カルマンフィルタで求まった予測分布を用いて、上記の制御方法１の解釈で最適化目的が規格範囲内である場合のベイズ最適化で制御を行えば、基本的にはうまく動作する。しかし、制御方法１での制御を行うと、例えば図３３に示すような、観測値が規格中央を行き過ぎてしまうオーバーシュート現象が発生してしまう。これは、カルマンフィルタには、図３４に示されるように、観測値の変動に過度に影響されずに予測（予測分布を推定）する特性すなわちノイズキャンセリング効果があるためである。つまり、カルマンフィルタを用いると、予測分布の挙動が観測値の挙動よりも小さく算出されやすい。このため、制御方法１のように規格中央を最適化目標値にして制御する場合、観測値が規格中央を跨いだ後でも、予測値はまだ跨いでないと判断され、予測値はそのまましばらく上昇または下降を続けるように算出される。そして、予測値が規格中央を跨いだときに初めて予測値の軌道修正がなされるように算出される。
　＜時系列データの制御方法２＞
　時系列データの制御方法１で発生してしまうオーバーシュート現象を抑制するために、最適化目的が規格範囲内である場合の解釈を変えた次善策である制御方法２を用いる。すなわち、時系列データの制御方法２では、最適化目的が規格範囲内である場合、直前の特性点（観測値）が規格中央より上側のときは最適化目標値を規格下限値に、直前の特性点が（観測値）が規格中央より下側のときは最適化目標値を規格上限値に設定する制御を行ってもよい。なお、最適化方法（最適化目的）については制御方法１と同様であり、直前の特性点（観測値）の位置が規格中央より上側なら最適化方法（最適化目的）を最小化、直前の特性点の位置が規格中央より下側なら最適化目的（最適化目的）を最大化に切り替える制御を行えばよい。
　図３５Ａおよび図３５Ｂは、本実施の形態に係る最適化目標値を規格上下限に設定する制御方法２の制御イメージを説明するための図である。なお、図３５Ａおよび図３５Ｂでも、説明を簡潔にするため製品特性数が１次元である場合の制御イメージが示されている。
　時系列データの制御方法２では、最適化目的が規格範囲内であり、図３５Ａに示す例のように直前の特性点ｙｔが規格中央よりも上側である場合、規格下限を最適化目標値に設定して制御を行う。一方、図３５Ｂに示すように直前の特性点ｙｔ＋１が規格中央よりも下側である場合、規格上限を最適化目標値に設定して制御を行う。
　このように、時系列データの制御方法２では、最適化目的が規格範囲内であることを、最適化目標値が直前の特性点（観測点）の位置に応じて規格上限値または規格下限値に切り替える制御方法であると解釈する。さらに、時系列データの制御方法２では、最適化目的が規格範囲内である場合、最適化方法が直前の特性点（観測点）の位置に応じて最大化または最小化に切り替わる制御方法であると解釈する。
　これにより、制御方法２では、予測値が規格中央を跨いだときに初めて予測値の軌道修正がなさるように算出される制御方法１と比較して、予測値の軌道修正のタイミングが早くなるので、オーバーシュート現象を抑制しうる。
　しかしながら、制御方法２での制御を行うと、例えば図３３に示すような、観測値が規格中央付近で振動するハンチング現象が発生しやすくなってしまう。なぜなら、制御方法２での制御を行うと、最適化目標値が規格中央から遠い規格下限値または規格上限値の位置に設定されるため、複数の候補制御点のうちから特性点が急激に上昇または下降するような候補制御点が次に設定する制御点として常に選択されやすくなるためである。
　＜時系列データの制御方法３＞
　最適化目的が規格範囲内である場合において、制御方法１のように最適化目標値を規格中央に設定するとオーバーシュート現象が生じ、制御方法２のように最適化目標値を規格上下限に設定するとハンチング現象が発生してしまう。そこで、最適化目標値の初期値を規格上下限に設定し、以下で説明する規則に従って、最適化目標値を徐々に規格中央に移動させていく制御方法３を行うことで、オーバーシュート現象およびハンチング現象を抑制できることが期待される。つまり、時系列データの制御方法３では、最適化目的が規格範囲内である場合、最適化目標値を規格上下限値から規格中央に寄せながら、最適化方法（最適化目的）を直前の特性値の位置に応じて最大化または最小化に切り替える制御を行う。
　図３６Ａおよび図３６Ｂは、本実施の形態に係る最適化目標値を徐々に規格中央に移動させていく制御方法３の制御イメージを説明するための図である。なお、図３６Ａおよび図３６Ｂでも、説明を簡潔にするため製品特性数が１次元である場合の制御イメージが示されている。
　時系列データの制御方法３では、最適化目的が規格範囲内である場合、最適化目標値を徐々に規格中央に移動させながら制御を行う。例えば図３６Ａに示すように直前の特性点ｙｔが規格中央よりも上側であり、かつ、過去の制御において特性点ｙｔ－４が規格範囲内で極小値である場合、その極小値を最適化目標値に設定して制御を行う。なお、最適化方法（最適化目的）は、直前の特性点ｙｔが規格中央よりも上側であるので最小化であるとして制御される。一方、図３６Ｂに示す例のように、直前の特性点ｙｔ＋１が規格中央よりも下側であり、かつ、過去の制御において特性点ｙｔ－１が規格範囲内で極大値である場合、その極大値を最適化目標値に設定して制御を行えばよい。なお、最適化方法（最適化目的）は、直前の特性点ｙｔ＋１が規格中央よりも下側であるので最大化であるとして制御される。
　以下では、最適化目標値を規格中央へ移動させる方法について説明する。
　一番平易な発想としては、規格上下限から規格中央への距離を等分割し、時刻ごとに一定の速度で移動させる方法がある。しかし、この方法では、分割数つまり移動させる速度を解析者が決定する必要がある。最適化目標値を移動させる速度が早すぎると、初めから最適化目標値を規格中央に設定した場合とほぼ等価になるため、オーバーシュート現象が発生しやすくなる。逆に、最適化目標値を移動させる速度が遅すぎると、初めから最適化目標値を規格上下限に設定した場合とほぼ等価になるため、ハンチング現象が発生しやすくなる。したがって、オーバーシュート現象およびハンチング現象をともに抑制する分割数を解析者が選択しなければならないという課題が生じる。また、オーバーシュート現象およびハンチング現象をともに抑制する分割数は、対象のタスクによって様々であり、解析者による試行錯誤が発生するため、現実的でない。
　それに対して、上記の制御方法３では、特性点の推移が上昇から下降に転じたら、極大点を次に上昇する際の最適化目標値に設定し、逆に特性点の推移が下降から上昇に転じたら、極小点を次に下降する際の最適化目標値に設定する。このように、制御方法３での制御を行うと、解析者に依存する作業は発生させずに、オーバーシュート現象およびハンチング現象を抑制できる。
　＜ＥＨＶＩＣを活用した制御方法３＞
　上述した制御方法３では、初期段階の極大点または極小点の位置が規格中央に近い場合、初めから最適化目標値を規格中央に設定した場合とほぼ等価になる。このため、初期段階の極大点または極小点の位置が規格中央に近い場合、オーバーシュート現象が発生する可能性がある。
　そこで、時系列データの制御方法３にベイズ最適化の獲得関数であるＥＨＶＩＣのアイデアを活用する。なお、以下では、製品特性数が２次元以上の場合の制御を例に挙げて説明する。
　ＥＨＶＩＣのアイデアでは、複数の製品特性を同時に最大化または最小化する際、観測済みのパレート点をもとに最適解の探索が効率的に進むように、評価すべき特性空間内の領域を区画化する。例えば２つの製品特性があるとき一方の製品特性が極端に最適化目標値に近い特性点を獲得しても、その値で特性空間全域に亘って区画化しているわけではなく、全パレート点の座標で段階的に区画化する。このＥＨＶＩＣのアイデアを活用して、時系列データの制御方法３において、最適化目標値が初期値である規格上下限から規格中央へ急激に移動してしまうことを抑制してもよい。
　しかし、上述したようにＥＨＶＩＣは、製品特性数に応じて、計算コストが指数関数的に増大するという問題がある。そこで、アクティブ領域を縮小させる規則である領域縮小規則を利用することで、探索精度を保ったまま、計算コストを多項式関数オーダーまで削減できる。
　なお、製品特性数が２次元以上の場合、上述した領域縮小規則をそのまま適用すれば最適化目標値を一意に決定できない。したがって、以下では、製品特性数が２次元以上の場合に、上述した領域縮小規則を適用しても最適化目標値を一意に決定できる方法について説明する。
　図３７は、製品特性数が２次元以上の場合に最適化目標値の候補が複数存在することを説明するための図である。
　製品特性数がＤ次元の場合、各製品特性に対して、規格中央より上側および下側に分けることで、特性空間全体は２Ｄ個の分割領域に分割することができる。分割領域それぞれにおいて、領域縮小規則を適用し、各製品特性について反対側のパレート境界を最適化目標値として設定する。図３７には、製品特性数が２次元の場合の４個の分割領域が示されている。また、図３７には、４個の分割領域それぞれに存在する観測済み特性点から算出されたパレート境界とアクティブ領域とが示されている。なお、図３７に示されるパレート境界は、最適化目標値を一意に定められない暫定的なパレート境界であるため、以下では暫定パレート境界と称することにする。また、暫定パレート境界のもとでのアクティブ領域を、以下では暫定アクティブ領域と称することにする。さらに、暫定的でないパレート境界を結合パレート境界と称し、結合パレート境界のもとでのアクティブ領域を、結合アクティブ領域と称することにする。
　例えば、第１製品特性（Ｙ１）における直前の特性点が規格中央より上側に位置する場合に領域縮小規則を適用するときには、第１製品特性（Ｙ１）において最適化方法（最適化目的）は最小化に決定し、最適化目標値は規格中央より下側領域に存在するパレート境界から定めればよい。しかし、例えば図３７に示すように、第１製品特性（Ｙ１）においてＡとＢとの２通りで示される暫定パレート境界がある場合、最適化目標値となる候補が２つ存在してしまうことになる。この場合には、ＡとＢとで示される暫定パレート境界のうち、規格中央から近い方（Ａの暫定パレート境界）、規格中央から遠い方（Ｂの暫定パレート境界）または、これらの中間的な位置の境界などを結合パレート境界として定めるとよい。なお、第２製品特性（Ｙ２）においても同様のことが言えるため説明を省略する。
　図３８Ａは、中心が規格中央と一致しない場合の結合アクティブ領域を示す図である。図３８Ｂは、中心が規格中央と一致する場合の結合アクティブ領域を示す図である。
　また、図３７に示す例において、第１製品特性（Ｙ１）の規格中央より上側および下側、第２製品特性（Ｙ２）の規格中央より上側および下側で、個別に結合パレート境界を定めない方がよい。個別に定めた結合パレート境界のもとでの結合アクティブ領域は、図３８Ａに示すように、単一の矩形となるものの、その中心は規格中央と異なる位置となる可能性が高いからである。また、結合アクティブ領域の中心が規格中央と異なる位置となる場合に制御方法３で時系列データの制御を実行したときには、結合アクティブ領域の中央に向かって観測値が追い込まれるように制御が進行する。このため、本来のターゲット（最適化目標値）である規格中央とずれて制御されてしまう可能性が高くなる。
　このことを鑑みて、本実施の形態におけるＥＨＶＩＣを活用した制御方法３では、各製品特性における結合パレート境界を、規格中央より上側領域と下側領域とにおいて合わせて一意に定める。例えば図３７に示す例では、第１製品特性（Ｙ１）において、規格中央より上側領域に存在する２つの暫定パレート境界、および規格中央より下側領域に存在する２つの暫定パレート境界、計4つのY1座標に対して、規格中央からの距離の平均値などの中間値だけ規格中央より大きいY1座標を上側領域の結合パレート境界として定めればよい。同様にその中間値だけ規格中央より小さいY1座標を下側領域の結合パレート境界として定めればよい。このようにして、上側領域の結合パレート境界とした下側領域の結合パレート境界とは、第１製品特性（Ｙ１）の規格中央から等距離になるように定められる。第２製品特性（Ｙ２）においても同様のため説明を省略する。
　このように、規格中央より上側領域と下側領域とにおいて合わせて結合パレート境界を定めることで、定められた結合パレート境界のもとでの結合アクティブ領域を、図３８Ｂに示すように、単一の矩形にさせ、かつ、その中心を規格中央と一致させることができる。
　なお、規格中央より上側領域と下側領域とにおいて合わせて結合パレート境界を定める方法としては、複数の暫定パレート境界のY1座標の中間位置を用いる場合に限らず、規格中央から距離が一番近い値または一番遠い値となる位置を通る結合パレート境界を定めてもよい。一番近い値となる位置（Y1座標）を通る結合パレート境界を定めた場合には、オーバーシュート現象が発生しやすくなり、一番遠い値となる位置（Y1座標）を通る結合パレート境界を定めた場合には、ハンチング現象が発生しやすくなる。このため、特に理由がなければ、上述した中間的な値となる位置を通る結合パレート境界を定めればよい。
　以上のように、単にカルマンフィルタとベイズ最適化とを組み合わせるだけでは、オーバーシュート現象が発生してしまう。そこで、時系列データの制御方法３を行うことで、最適化目標値の初期値を規格上下限に設定し、かつ徐々に規格中央へ寄せることでハンチング現象も抑制できる。また、時系列データの制御方法３によれば、候補制御点の定量評価のためにベイズ最適化の獲得関数を利用し、さらに領域縮小規則を適用することで計算コストを抑制した高速化アルゴリズムを実現できる。
　このように、カルマンフィルタで求まる予測分布と予測分布からベイズ最適化を用いて求まる評価値とを逐次的に算出することで、次に設定する候補制御点をランキング順などで解析者（ユーザ）に提供され、選択される。よって、解析者に依らず、オーバーシュート現象またはハンチング現象の発生を抑制できる候補制御点を定量的に評価できるので、オーバーシュート現象またはハンチング現象の発生を抑制し、安定した効率のよい時系列データのリアルタイム制御を実現できる。
　［効果等］
　以上のように、本実施の形態における評価装置は、第１時刻における制御済みの制御点に対応する既知の特性点に基づいて、前記第１時刻に続く第２時刻における複数の候補制御点に対応する複数の未知の特性点をベイズ最適化によって評価する評価装置であって、前記第１時刻における制御済みの制御点および前記第１時刻における既知の特性点を示す制御結果データを取得する第１受信手段と、前記未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有し、当該最適化目的を示す目的データを取得する第２受信手段と、前記少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データを取得する第３受信手段と、少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、前記分割方法により分割された特性空間の領域ごとに示す領域縮小規則データを取得する第４受信手段と、前記制御結果データ、前記目的データ、前記制約条件データおよび前記領域縮小規則データに基づいて、前記複数の未知の特性点それぞれの評価値を算出する算出手段と、前記評価値を出力する出力手段と、を備え、前記算出手段は、前記制約条件の適合度合いに応じた重み付けを、前記少なくとも１つの製品特性に対する評価値に付与する。
　これにより、算出手段が、制御結果データ、目的データ、制約条件データおよび領域縮小規則データに基づいて、第１時刻に続く第２時刻における複数の未知の特性点の評価値を算出することができる。複数の未知の特性点の評価値を算出するときには、制約条件の適合度合いに応じた重み付けを、少なくとも１つの製品特性に対する評価値に付与し、この少なくとも１つの製品特性は最適化目的を有する。したがって、最適化問題の目的を有する時系列の製品特性に対して制約条件が付与されている最適化問題に対して、カルマンフィルタとベイズ最適化を適用することができる。このようにして、過去の制御結果情報である第１時刻における制御済みの制御点に対応する既知の特性点から、複数の未知の特性点の評価値を算出できるので、次に設定すべき候補制御点を定量的に解析することができる。
　また、本実施の形態では、制約条件は規格範囲であり、最適化目的には、製品特性を規格範囲内に収める第１目的と、製品特性を最小化または最大化する第２目的とがある。そして、評価値算出部１２は、少なくとも１つの製品特性のそれぞれについて、（ｉ）評価値を算出するために用いられる当該製品特性の区間が規格範囲外にある場合と、（ｉｉ）その区間が規格範囲内であって、かつ、最適化目的が第１目的である場合と、（ｉｉｉ）その区間が規格範囲内であって、かつ、最適化目的が第２目的である場合とで、互いに異なる重み付け処理を行うことによって、評価値を算出する。つまり、上記（式６）および（式７）に基づいて評価値が算出される。これにより、製品特性の最適化目的が第１目的であっても、第２目的であっても、候補制御点の評価値をベイズ最適化に基づいて適切に算出することができる。つまり、製品特性の最適化目的が規格範囲内であっても、最大化または最小化であっても、候補制御点の評価値をベイズ最適化に基づいて適切に算出することができる。また、（ｉｉｉ）の場合は、製品特性の区間が規格範囲内であって、かつ、最適化目的が第２目的であるため、非特許文献２の手法とは異なり、最適化目的が最大化または最小化である製品特性に、制約条件として規格範囲がある場合であっても、評価値を定量的に適切に算出することができる。
　その結果、規格範囲などの制約条件がある最適化問題に対しても適用することができる。つまり、適用場面を拡張し、最適解の探索効率向上のための定量評価を行うことができる。
　また、前記算出手段は、カルマンフィルタを用いて前記複数の候補制御点における予測分布を算出し、算出した前記予測分布を用いて、前記評価値を算出してもよい。
　これにより、カルマンフィルタで求まる予測分布と予測分布からベイズ最適化を用いて求まる評価値とを逐次的に算出することができる。よって、次に制御点として設定する候補制御点を評価値に基づいて選択できる。つまり、各候補制御点に対して算出された評価値が出力されるため、評価装置のユーザは、それらの評価値に基づいて候補制御点を次の制御点として選択し、その制御点を用いた制御によって得られる特性点を、次の各候補制御点の評価値の算出に利用することができる。このような制御と評価値の算出および出力との繰り返しによって、各製品特性の最適化目的を満たす候補制御点の解、すなわち最適解を得ることができる。
　また、獲得関数の算出において、パレート点からある規則によって定まるパレート境界と、そのパレート境界で特性空間全体がアクティブ領域および非アクティブ領域に区分けされる。したがって、算出手段は、さらに、領域縮小規則データに基づき、パレート境界の定義を変更できるので、単一の超直方体の体積の期待値の算出のみで、ベイズ最適化における獲得関数を算出することができるようになる。これにより、算出手段は、探索効率を保ったまま獲得関数の計算量を抑制できるので、未知の特性点の評価を高速で行えるようになる。
　また、評価値算出部１２は、モンテカルロ法を用いて各候補制御点の評価値を算出してもよい。これにより、モンテカルロ法が近似手法であるため、評価値が解析的に算出困難な場合でも、近似的に算出することができる。つまり、（式６）および（式７）の演算を厳密に行わずに、モンテカルロ法を用いることによって、近似的に算出することができる。なお、近似手法であれば、モンテカルロ法に限らず、他の手法が用いられてもよい。
　（実施例）
　上記の制御方法３に領域縮小規則を適用して製品特性の最適解の探索を行う際の処理の実施例として、制御結果が一組追加された際に、領域縮小規則からアクティブ領域を定める処理の一例を説明する。
　図３９は、本実施の形態の実施例に係る最適解の探索を行う際に得られた制御結果データシートの一例を示す図である。図４０Ａは、図３９の２６番目までの制御結果を得た時点での結合アクティブ領域を示す図である。図４０Ｂは、図３９の２７番目の制御結果を得た時点での結合アクティブ領域を示す図である。なお、カルマンフィルタを用いた予測分布の算出処理およびＥＨＶＩおよびＥＨＶＩＣを用いた獲得関数の算出処理は、上述した通りであるので説明を省略する。また、図４０Ａおよび図４０Ｂにおいて規格外の特性点は省略している。
　図３９には、制御番号ごとに、その制御番号によって識別される制御で用いられた制御点のプロセス条件と、その制御によって得られた制御結果である特性点の製品特性とが示されている。また、図３９には、プロセス条件の数が２つ（第１プロセス条件、第２プロセス条件）、製品特性の数が２つ（第１製品特性、第２製品特性２）である場合が示されている。最適化目的は第１製品特性（Ｙ１）第２製品特性（Ｙ２）とも最小化であり、第１製品特性（Ｙ１）の規格範囲を１．１～１．７、第２製品特性（Ｙ２）の規格範囲を０．５～０．７とする。また、結合アクティブ領域の初期値は規格範囲、つまり結合パレート境界の初期値は各次元の規格上下限とする。本実施例における規格範囲は、（式１１）のように表記できる。また、図４０Ａおよび図４０Ｂにおける領域縮小規則が適用される領域（領域Ｒ１－１～領域Ｒ４－２）は、（式１２）のように定めている。
　領域Ｒ１－１および領域Ｒ２－１においては、各領域内の特性点のうち、第１製品特性（Ｙ１）の座標の最小値を暫定パレート境界とする。領域Ｒ３－１および領域Ｒ４－１においては、各領域内の特性点のうち、第１製品特性（Ｙ１）の座標の最大値を暫定パレート境界とする。領域Ｒ１－２および領域Ｒ４－２においては、各領域内の特性点のうち、第２製品特性（Ｙ２）の座標の最小値を暫定パレート境界とする。領域Ｒ２－２および領域Ｒ３－２においては、各領域内の特性点のうち、第２製品特性（Ｙ２）の座標の最大値を暫定パレート境界とする。
　また、各次元について暫定パレート境界の規格中央からの距離の平均値だけ規格中央から離れた座標を通る境界を結合パレート境界とし、それにより定まる矩形領域を結合アクティブ領域として定める。また、以下では、第１製品特性（Ｙ１）における座標をY1座標と称し、第２製品特性（Ｙ２）における座標をY2座標と称する。
　まず、図３９に示す制御結果データにおいて、制御番号が２６番目までの制御結果が与えられているとする。この時点での結合アクティブ領域は、図４０Ａに示す結合アクティブ領域となる。
　図４０Ａに示す特性点Ａ，Ｂ，Ｃ，Ｄ，Ｅ，Ｆ，Ｇは、図４０Ａに示す領域Ｒ１－１，領域Ｒ１－２，領域Ｒ２－１，領域Ｒ２－２，領域Ｒ３－１，領域Ｒ３－２，領域Ｒ４－１におけるパレート点である。なお、図４０Ａに示す例では、領域Ｒ４－２にはパレート点はない。また、図４０Ａに示す例では、第１製品特性（Ｙ１）の結合パレート境界は、特性点Ａ，Ｃ，Ｅ，ＧのY1座標と第１製品特性（Ｙ１）における規格中央との距離の平均値だけ規格中央から離れたY1座標を通る位置となっている。第２製品特性（Ｙ２）の結合パレート境界は、特性点B,D,FのY2座標と第２製品特性（Ｙ２）における規格上限値と第２製品特性（Ｙ２）における規格中央との距離の平均値だけ規格中央から離れたY2座標を通る位置となっている。
　具体的には、図４０Ａに示す結合アクティブ領域の第１製品特性（Ｙ１）における上限と下限とは、次の通りである。すなわち、第１製品特性（Ｙ１）の上限は、特性点Ａ，Ｃ，Ｅ，ＧのY1座標（1.4256, 1.4181, 1.3791, 1.3927）と規格中央（1.4000）との距離の平均値（0.017975）だけ規格中央から加算したY1座標（1.417975）となる。また、第１製品特性（Ｙ１）の下限は、特性点Ａ，Ｃ，Ｅ，ＧのY1座標（1.4256, 1.4181, 1.3791, 1.3927）と規格中央（1.4000）との距離の平均値（0.017975）だけ規格中央から減算したY1座標（1.382025）となる。また、図４０Ａに示す結合アクティブ領域の第２製品特性（Ｙ２）における上限と下限とは、次の通りである。すなわち、第２製品特性（Ｙ２）の上限は、特性点Ｂ，Ｄ，ＦのY2座標（0.6173, 0.5790, 0.5239）および規格上限値（0.7）と規格中央（0.6000）との距離の平均値（0.053600）だけ、規格中央から加算したY2座標（0.653600）となる。また、第２製品特性（Ｙ２）の下限は、特性点Ｂ，Ｄ，ＦのY2座標（0.6173, 0.5790, 0.5239）および規格下限値（0.5）と規格中央（0.6000）との距離の平均値（0.053600）だけ規格中央から減算したY2座標（0.546400）となる。
　図３９に示すように、直近となる制御番号が２６番目の特性点は、規格範囲内であり、領域Ｒ２－１に属する。なお、制御番号が２６番目の特性点は、図４０Ａでは特性点Ｃに該当し、第１製品特性（Ｙ１）においては規格中央より上側、第２製品特性（Ｙ２）においては規格中央より下側に位置する。したがって、第１製品特性（Ｙ１）および第２製品特性（Ｙ２）における最適化方法（最適化目的）は、最小化および最大化に設定される。また、第１製品特性（Ｙ１）における最適化目標値は、結合アクティブ領域の第１製品特性（Ｙ１）の下限の座標（1.382025）に設定される。第２製品特性（Ｙ２）における最適化目標値は、結合アクティブ領域の第２製品特性（Ｙ２）の上限の座標（0.653600）に設定される。
　次に、評価値算出部１２は、制御番号が２６番目までの制御結果をもとに、（式５）などのカルマンフィルタの更新式を用いて、候補制御点ごとに２７番目の特性点の予測分布を算出する。評価値算出部１２は、算出した予測分布と（式７）などのベイズ最適化の更新式とを用いて、候補制御点ごとに獲得関数を算出する。算出した獲得関数から得た評価値が最大となる候補制御点が２７番目の制御点として採用される。
　次に、図３９に示す制御結果データにおいて、制御番号が２７番目の制御結果が得られると、同様の処理を実行する。この時点での結合アクティブ領域は、図４０Ｂに示す結合アクティブ領域となる。
　図３９に示すように、直近となった制御番号が２７番目の特性点は、規格範囲内であり、領域Ｒ４－２に属する。制御番号が２７番目の特性点は、図４０Ｂでは特性点Ｈに該当し、領域Ｒ４－２に属する新規パレート点となる。したがって、領域Ｒ４－２においてアクティブ領域を削ったとしても結合アクティブ領域の第１製品特性（Ｙ１）における上下限は変化しないので更新されず、結合アクティブ領域の第２製品特性（Ｙ２）における上下限が更新される。第２製品特性（Ｙ２）における上限は、特性点Ｂ，Ｄ，Ｆ，ＨのY2座標（0.6173, 0.5790, 0.5239, 0.6074）と規格中央（0.6000）との距離の平均値（0.030450）だけ規格中央から加算したY2座標（0.630450）となる。第２製品特性（Ｙ２）における下限は、特性点Ｂ，Ｄ，Ｆ，ＨのY2座標（0.6173, 0.5790, 0.5239, 0.6074）と規格中央（0.6000）との距離の平均値（0.030450）だけ規格中央から減算したY2座標（0.569550）となる。
　なお、制御番号が２７番目の特性点は、図４０Ｂでは特性点Ｈに該当し、第１製品特性（Ｙ１）においては規格中央より下側、第２製品特性（Ｙ２）においては規格中央より上側に位置する。
　したがって、第１製品特性（Ｙ１）および第２製品特性（Ｙ２）における最適化方法（最適化目的）は、最大化および最小化に設定される。また、第１製品特性（Ｙ１）における最適化目標値は、結合アクティブ領域の第１製品特性（Ｙ１）の上限の座標（1.417975）に設定される。第２製品特性（Ｙ２）における最適化目標値は、結合アクティブ領域の第２製品特性（Ｙ２）の下限の座標（0.569550）に設定される。
　評価値算出部１２は、制御番号が２７番目までの制御結果をもとに、（式５）などのカルマンフィルタの更新式を用いて、候補制御点ごとに２８番目の特性点の予測分布を算出する。評価値算出部１２は、算出した予測分布と（式７）などのベイズ最適化の更新式とを用いて、候補制御点ごとに獲得関数を算出する。算出した獲得関数から得た評価値が最大となる候補制御点が２８番目の制御点として採用される。
　以上のような処理を、制御結果が一組追加されるごとに行うことで、最適化目標値は初期値である規格上下限から徐々に規格中央へ移動させていく制御方法３を行う。これにより、時系列データの制御において人に依存しない高精度かつ高効率な制御が可能となる。
　以上、本開示の一態様に係る評価装置１００について、上記実施の形態および各変形例に基づいて説明したが、本開示は、その実施の形態および各変形例に限定されるものではない。本開示の趣旨を逸脱しない限り、当業者が思いつく各種変形を上記実施の形態または各変形例に施したものも本開示に含まれてもよい。
　なお、上記実施の形態等において、各構成要素は、専用のハードウェアで構成されるか、各構成要素に適したソフトウェアプログラムを実行することによって実現されてもよい。各構成要素は、ＣＰＵまたはプロセッサなどのプログラム実行部が、ハードディスクまたは半導体メモリなどの記録媒体に記録されたソフトウェアプログラムを読み出して実行することによって実現されてもよい。ここで、上記実施の形態等の評価装置などを実現するソフトウェアは、例えば図１０に示すフローチャートの各ステップをコンピュータに実行させるプログラムである。
　なお、以下のような場合も本開示に含まれる。
　（１）上記の少なくとも１つの装置は、具体的には、マイクロプロセッサ、ＲＯＭ（Read Only Memory）、ＲＡＭ（Random Access Memory）、ハードディスクユニット、ディスプレイユニット、キーボード、マウスなどから構成されるコンピュータシステムである。そのＲＡＭまたはハードディスクユニットには、コンピュータプログラムが記憶されている。マイクロプロセッサが、コンピュータプログラムにしたがって動作することにより、上記の少なくとも１つの装置は、その機能を達成する。ここでコンピュータプログラムは、所定の機能を達成するために、コンピュータに対する指令を示す命令コードが複数個組み合わされて構成されたものである。
　（２）上記の少なくとも１つの装置を構成する構成要素の一部または全部は、１個のシステムＬＳＩ（Large Scale Integration：大規模集積回路）から構成されているとしてもよい。システムＬＳＩは、複数の構成部を１個のチップ上に集積して製造された超多機能ＬＳＩであり、具体的には、マイクロプロセッサ、ＲＯＭ、ＲＡＭなどを含んで構成されるコンピュータシステムである。前記ＲＡＭには、コンピュータプログラムが記憶されている。マイクロプロセッサが、コンピュータプログラムにしたがって動作することにより、システムＬＳＩは、その機能を達成する。
　（３）上記の少なくとも１つの装置を構成する構成要素の一部または全部は、その装置に脱着可能なＩＣカードまたは単体のモジュールから構成されているとしてもよい。ＩＣカードまたはモジュールは、マイクロプロセッサ、ＲＯＭ、ＲＡＭなどから構成されるコンピュータシステムである。ＩＣカードまたはモジュールは、上記の超多機能ＬＳＩを含むとしてもよい。マイクロプロセッサが、コンピュータプログラムにしたがって動作することにより、ＩＣカードまたはモジュールは、その機能を達成する。このＩＣカードまたはこのモジュールは、耐タンパ性を有するとしてもよい。
　（４）本開示は、上記に示す方法であるとしてもよい。また、これらの方法をコンピュータにより実現するコンピュータプログラムであるとしてもよいし、コンピュータプログラムからなるデジタル信号であるとしてもよい。
　また、本開示は、コンピュータプログラムまたはデジタル信号をコンピュータで読み取り可能な記録媒体、例えば、フレキシブルディスク、ハードディスク、ＣＤ（Compact Disc）－ＲＯＭ、ＤＶＤ、ＤＶＤ－ＲＯＭ、ＤＶＤ－ＲＡＭ、ＢＤ（Blu-ray（登録商標） Disc）、半導体メモリなどに記録したものとしてもよい。また、これらの記録媒体に記録されているデジタル信号であるとしてもよい。
　また、本開示は、コンピュータプログラムまたはデジタル信号を、電気通信回線、無線または有線通信回線、インターネットを代表とするネットワーク、データ放送等を経由して伝送するものとしてもよい。
　また、プログラムまたはデジタル信号を記録媒体に記録して移送することにより、またはプログラムまたはデジタル信号を、ネットワーク等を経由して移送することにより、独立した他のコンピュータシステムにより実施するとしてもよい。
　本開示の評価装置によれば、時系列データの制御問題に対して、オーバーシュート現象またはハンチング現象を抑制できる候補制御点を定量的に評価できる。
　本開示の評価装置は、製品特性に規格範囲による制約条件がある場合でも、最適プロセス条件を定量的に効率良く探索できるという効果を奏し、工業製品の量産現場に限らず、空間シミュレーションや動的物体の軌道制御等の最適制御の用途にも適用できる。
　１０　　受信制御部
　１１　　候補制御点作成部
　１２　　評価値算出部
　１３　　評価値出力部
　１００　　評価装置
　１０１ａ　　入力部
　１０１ｂ　　通信部
　１０２　　演算回路
　１０３　　メモリ
　１０４　　表示部
　１０５　　記憶部
　２００　　プログラム
　２０１　　特性点データ
　２１０　　設定情報
　２１１　　プロセス条件データ
　２１２　　目的データ
　２１３　　制約条件データ
　２１４　　領域縮小規則データ
　２２１　　候補制御点データ
　２２２　　制御結果データ
　２２３　　予測分布データ
　２２４　　評価値データ
　３００　　受付画像
　３１０　　プロセス条件領域
　３１１～３１４、３２１～３２８　　入力フィールド
　３２０　　製品特性領域
　３３０　　領域縮小規則領域

Claims (8)

  1. 　第１時刻における制御済みの制御点に対応する既知の特性点に基づいて、前記第１時刻に続く第２時刻における複数の候補制御点に対応する複数の未知の特性点をベイズ最適化によって評価する評価装置であって、
    　前記第１時刻における制御済みの制御点および前記第１時刻における既知の特性点を示す制御結果データを取得する第１受信手段と、
    　前記未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有し、当該最適化目的を示す目的データを取得する第２受信手段と、
    　前記少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データを取得する第３受信手段と、
    　少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、前記分割方法により分割された特性空間の領域ごとに示す領域縮小規則データを取得する第４受信手段と、
    　前記制御結果データ、前記目的データ、前記制約条件データおよび前記領域縮小規則データに基づいて、前記複数の未知の特性点それぞれの評価値を算出する算出手段と、
    　前記評価値を出力する出力手段と、を備え、
    　前記算出手段は、前記制約条件の適合度合いに応じた重み付けを、前記少なくとも１つの製品特性に対する評価値に付与する、
    　評価装置。
  2. 　前記制約条件は、少なくとも１つの制約範囲であり、
    　前記最適化目的には、製品特性を前記少なくとも１つの制約範囲のうちの何れかの制約範囲内に収める第１目的と、製品特性を最小化または最大化する第２目的とがあり、
    　前記算出手段は、
    　前記少なくとも１つの製品特性のそれぞれについて、
    　（ｉ）前記評価値を算出するために用いられる当該製品特性の区間が前記少なくとも１つの制約範囲のそれぞれの外にある場合と、
    　（ｉｉ）前記区間が前記少なくとも１つの制約範囲のうちの何れかの制約範囲内であって、かつ、前記最適化目的が前記第１目的である場合と、
    　（ｉｉｉ）前記区間が前記少なくとも１つの制約範囲のうちの何れかの制約範囲内であって、かつ、前記最適化目的が前記第２目的である場合とで、
    　互いに異なる重み付け処理を行うことによって、前記評価値を算出する、
    　請求項１に記載の評価装置。
  3. 　複数のプロセス条件のそれぞれの所定の条件を満たす値を組み合わせることによって、前記複数の候補制御点を作成する候補制御点作成手段をさらに備える、
    　請求項２に記載の評価装置。
  4. 　前記算出手段は、
    　前記少なくとも１つの制約範囲のうち、矩形と異なる形状の制約範囲に基づいて、前記評価値を算出する、
    　請求項２に記載の評価装置。
  5. 　前記算出手段は、
    　カルマンフィルタを用いて前記複数の候補制御点における予測分布を算出し、算出した前記予測分布を用いて、前記評価値を算出する、
    　請求項１～４のいずれか１項に記載の評価装置。
  6. 　前記算出手段は、
    　モンテカルロ法を用いて前記評価値を算出する、
    　請求項１～４のいずれか１項に記載の評価装置。
  7. 　評価装置が、第１時刻における制御済みの制御点に対応する既知の特性点に基づいて、前記第１時刻に続く第２時刻における複数の候補制御点に対応する複数の未知の特性点をベイズ最適化によって評価する評価方法であって、
    　前記第１時刻における制御済みの制御点および前記第１時刻における既知の特性点を示す制御結果データを取得する第１受信ステップと、
    　前記未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有し、当該最適化目的を示す目的データを取得する第２受信ステップと、
    　前記少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データを取得する第３受信ステップと、
    　少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、前記分割方法により分割された特性空間の領域ごとに示す領域縮小規則データを取得する第４受信ステップと、
    　前記制御結果データ、前記目的データ、前記制約条件データおよび前記領域縮小規則データに基づいて、前記複数の未知の特性点それぞれの評価値を算出する算出ステップと、
    　前記評価値を出力する出力ステップと、を含み、
    　前記算出ステップにおいて、前記制約条件の適合度合いに応じた重み付けを、前記少なくとも１つの製品特性に対する評価値に付与する、
    　評価方法。
  8. 　コンピュータが、第１時刻における制御済みの制御点に対応する既知の特性点に基づいて、前記第１時刻に続く第２時刻における複数の候補制御点に対応する複数の未知の特性点をベイズ最適化によって評価するためのプログラムであって、
    　前記第１時刻における制御済みの制御点および前記第１時刻における既知の特性点を示す制御結果データを取得する第１受信ステップと、
    　前記未知の特性点は、１または複数の製品特性の値を示し、少なくとも１つの製品特性が最適化目的を有し、当該最適化目的を示す目的データを取得する第２受信ステップと、
    　前記少なくとも１つの製品特性に対して付与された制約条件を示す制約条件データを取得する第３受信ステップと、
    　少なくとも１つの製品特性によって表現される特性空間の分割方法を示し、かつ、アクティブ領域を縮小する次元を、前記分割方法により分割された特性空間の領域ごとに示す領域縮小規則データを取得する第４受信ステップと、
    　前記制御結果データ、前記目的データ、前記制約条件データおよび前記領域縮小規則データに基づいて、前記複数の未知の特性点それぞれの評価値を算出する算出ステップと、
    　前記評価値を出力する出力ステップと、を前記コンピュータに実行させ、
    　前記算出ステップにおいて、前記制約条件の適合度合いに応じた重み付けを、前記少なくとも１つの製品特性に対する評価値に付与する、
    　プログラム。
PCT/JP2023/018114 2022-06-08 2023-05-15 評価装置、評価方法、およびプログラム WO2023238606A1 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2022-093250 2022-06-08
JP2022093250 2022-06-08

Publications (1)

Publication Number Publication Date
WO2023238606A1 true WO2023238606A1 (ja) 2023-12-14

Family

ID=89118230

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2023/018114 WO2023238606A1 (ja) 2022-06-08 2023-05-15 評価装置、評価方法、およびプログラム

Country Status (1)

Country Link
WO (1) WO2023238606A1 (ja)

Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005285090A (ja) * 2003-12-24 2005-10-13 Yamaha Motor Co Ltd 多目的最適化装置、多目的最適化方法および多目的最適化プログラム
JP2011048768A (ja) * 2009-08-28 2011-03-10 Hitachi Ltd 最適設計装置
JP2012123592A (ja) * 2010-12-08 2012-06-28 Fujitsu Ltd 最適化プログラム、装置及びプログラム
JP2018045266A (ja) * 2016-09-12 2018-03-22 株式会社日立製作所 設計支援装置
JP2022077760A (ja) * 2020-11-12 2022-05-24 富士通株式会社 情報処理プログラム、装置、及び方法

Patent Citations (5)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2005285090A (ja) * 2003-12-24 2005-10-13 Yamaha Motor Co Ltd 多目的最適化装置、多目的最適化方法および多目的最適化プログラム
JP2011048768A (ja) * 2009-08-28 2011-03-10 Hitachi Ltd 最適設計装置
JP2012123592A (ja) * 2010-12-08 2012-06-28 Fujitsu Ltd 最適化プログラム、装置及びプログラム
JP2018045266A (ja) * 2016-09-12 2018-03-22 株式会社日立製作所 設計支援装置
JP2022077760A (ja) * 2020-11-12 2022-05-24 富士通株式会社 情報処理プログラム、装置、及び方法

Similar Documents

Publication Publication Date Title
Janga Reddy et al. An efficient multi-objective optimization algorithm based on swarm intelligence for engineering design
Onken et al. Discretize-optimize vs. optimize-discretize for time-series regression and continuous normalizing flows
Örkcü et al. Estimating the parameters of 3-p Weibull distribution using particle swarm optimization: A comprehensive experimental comparison
Blanquero et al. Variable selection in classification for multivariate functional data
CN113574475A (zh) 确定用于控制环境的因果模型
JP7232122B2 (ja) 物性予測装置及び物性予測方法
US20240028795A1 (en) Design assitance device, design assitance method, and design assitance program
CN113597582A (zh) 使用因果模型调谐pid参数
CN102129242A (zh) 基于两层混合智能优化的批处理生产过程产品质量控制方法
Phong et al. PSO-convolutional neural networks with heterogeneous learning rate
Goswami et al. Variants of genetic algorithms and their applications
WO2023238606A1 (ja) 評価装置、評価方法、およびプログラム
Huang et al. Split-level evolutionary neural architecture search with elite weight inheritance
Pourtaheri et al. Stability investigation of multi-objective heuristic ensemble classifiers
JP7498393B2 (ja) 情報処理装置、情報処理方法、プログラム及び情報処理システム
CN114918581B (zh) 焊接参数处理方法、装置、存储介质及处理器
WO2022230736A1 (ja) 設計支援装置、設計支援方法及び設計支援プログラム
Niskanen Statistical approach to fuzzy cognitive maps
Hannah et al. Semiconvex regression for metamodeling-based optimization
WO2023181497A1 (ja) 評価装置、評価方法、およびプログラム
CN116029193A (zh) 预测半导体设备的特性的方法和执行该方法的计算设备
Li et al. Parameters selection for support vector machine based on particle swarm optimization
US20220391737A1 (en) Evaluation device, evaluation method, and recording medium
JP5581753B2 (ja) プラント制御装置、そのモデル予測制御装置
CN113597305A (zh) 使用因果模型制造生物药物

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 23819590

Country of ref document: EP

Kind code of ref document: A1