特開 2024-131176 生成方法、生成装置、および、プログラム

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】公開特許公報(A)

(11)【公開番号】P2024131176

(43)【公開日】2024-09-30

(54)【発明の名称】生成方法、生成装置、および、プログラム

(51)【国際特許分類】

   G06F 30/20 20200101AFI20240920BHJP

【ＦＩ】

G06F30/20

【審査請求】未請求

【請求項の数】6

【出願形態】ＯＬ

(21)【出願番号】P 2023041275

(22)【出願日】2023-03-15

(71)【出願人】

【識別番号】314012076

【氏名又は名称】パナソニックＩＰマネジメント株式会社

(74)【代理人】

【識別番号】100109210

【弁理士】

【氏名又は名称】新居 広守

(74)【代理人】

【識別番号】100137235

【弁理士】

【氏名又は名称】寺谷 英作

(74)【代理人】

【識別番号】100131417

【弁理士】

【氏名又は名称】道坂 伸一

(72)【発明者】

【氏名】竹澤 昌宏

(72)【発明者】

【氏名】原 伸夫

(72)【発明者】

【氏名】中橋 昭久

【テーマコード（参考）】

5B146

【Ｆターム（参考）】

5B146AA17

5B146DJ01

5B146DJ11

(57)【要約】

【課題】複数の目的変数を用いて高い精度でパラメータを推定する推定式を生成する。
【解決手段】物理現象の条件を示す条件値から、物理現象から得られる観測値を推定する推定式を生成する生成方法は、物理現象の時系列の観測値と、時系列の観測値の複数の時間区間それぞれにおける変化率を含む複数の第一変化率とを取得し（Ｓ１０３）、物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値であって、時系列の物性値と、時系列の物性値の複数の時間区間それぞれにおける変化率を含む複数の第二変化率とを少なくとも含む１以上の第二パラメータ値を取得し（Ｓ１１２）、時系列の観測値と、複数の第一変化率と、１以上の第一パラメータ値と、１以上の第二パラメータ値とを用いて推定式を生成し（Ｓ１１３）、出力する（Ｓ１１７）。
【選択図】図２

【特許請求の範囲】

【請求項1】

物理現象の条件を示す条件値から、前記物理現象から得られる観測値を推定する推定式を生成する生成方法であって、
前記物理現象の時系列の観測値と、前記時系列の観測値の複数の時間区間それぞれにおける変化率を含む複数の第一変化率とを取得し、
前記物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値であって、時系列の物性値と、前記時系列の物性値の複数の時間区間それぞれにおける変化率を含む複数の第二変化率とを少なくとも含む１以上の第二パラメータ値を取得し、
前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて前記推定式を生成して出力する
生成方法。

【請求項2】

前記複数の第一変化率を取得する際には、
前記時系列の観測値の変化率が第一条件を満たす時刻、または、前記時系列の観測値の変化率の変化率が第二条件を満たす時刻で時間軸を区切ることで得られる複数の時間区間を、前記複数の時間区間として用いて前記複数の第一変化率を取得する
請求項１に記載の生成方法。

【請求項3】

前記複数の第一変化率を取得する際には、
前記物理現象が生じている系を観測することで、前記系の状態が変化した１以上の時刻を取得し、
取得した前記１以上の時刻で時間軸を区切ることで得られる複数の時間区間を、前記複数の時間区間として用いて前記複数の第一変化率を取得する
請求項１に記載の生成方法。

【請求項4】

前記推定式を生成する際には、
前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて応答曲面法により前記推定式を生成する
請求項１～３のいずれか１項に記載の生成方法。

【請求項5】

物理現象の条件を示す条件値から、前記物理現象から得られる観測値を推定する推定式を生成する生成装置であって、
前記物理現象の時系列の観測値と、前記時系列の観測値の複数の時間区間それぞれにおける変化率を含む複数の第一変化率とを取得するデータ取得部と、
前記物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値であって、時系列の物性値と、前記時系列の物性値の複数の時間区間それぞれにおける変化率を含む複数の第二変化率とを少なくとも含む１以上の第二パラメータ値を取得するパラメータ取得部と、
前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて前記推定式を生成して出力する生成部とを備える
生成装置。

【請求項6】

請求項１に記載の生成方法をコンピュータに実行させるためのプログラム。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、生成方法、生成装置、および、プログラムに関する。

【背景技術】

【0002】

近年の計算機技術の進展により、物理現象を模擬するための数値解析あるいは数値シミュレーション（単にシミュレーションともいう）が様々な分野で広く用いられるようになっている。

【0003】

例えば、特許文献１では、複数の材質からなる被測定体の温度を計測することにより、測定が困難な熱物性を同定する手法について記述されている。

【0004】

物性値、または、境界条件等の特性値のように、同定される対象であるパラメータが多数ある場合、それらの多数のパラメータを一意に決定するには、多数の目的変数が必要である。

【先行技術文献】

【特許文献】

【0005】

【特許文献1】特開２００５－１４０６９３号公報

【発明の概要】

【発明が解決しようとする課題】

【0006】

しかしながら、パラメータの同定のために多数の目的変数を用意することが困難であることがある。多数の目的変数を用意することが困難である場合、パラメータを推定する推定式の精度が低下するという問題がある。

【0007】

そこで、本発明は、複数の目的変数を用いて高い精度でパラメータを推定する推定式を生成する生成方法等を提供する。

【課題を解決するための手段】

【0008】

本発明の一態様に係る生成方法は、物理現象の条件を示す条件値から、前記物理現象から得られる観測値を推定する推定式を生成する生成方法であって、前記物理現象の時系列の観測値と、前記時系列の観測値の複数の時間区間それぞれにおける変化率を含む複数の第一変化率とを取得し、前記物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値であって、時系列の物性値と、前記時系列の物性値の複数の時間区間それぞれにおける変化率を含む複数の第二変化率とを少なくとも含む１以上の第二パラメータ値を取得し、前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて前記推定式を生成して出力する生成方法である。

【0009】

なお、これらの包括的または具体的な態様は、システム、装置、集積回路、コンピュータプログラムまたはコンピュータ読み取り可能なＣＤ－ＲＯＭなどの記録媒体で実現されてもよく、システム、装置、集積回路、コンピュータプログラムおよび記録媒体の任意な組み合わせで実現されてもよい。

【発明の効果】

【0010】

本発明の生成方法は、複数の目的変数を用いて高い精度でパラメータを推定することができる。

【図面の簡単な説明】

【0011】

【図1】実施の形態における生成装置の構成を示すブロック図である。

【図2】実施の形態における生成方法を示すフロー図である。

【図3A】粉体の圧縮せん断試験を示す第一の説明図である。

【図3B】粉体の圧縮せん断試験を示す第二の説明図である。

【図3C】粉体の圧縮せん断試験を示す第三の説明図である。

【図4】シミュレーションにおけるパラメータと推定式との関係を示す説明図である。

【図5】実施例１における実測した応答の一例を示す説明図である。

【図6】実施例１における系の状態の例を示す説明図である。

【図7】実施例１における区間の一例を示す説明図である。

【図8】実施例１における部分実測データの一例を示す説明図である。

【図9】実施例１における部分実測データの傾きの別の例を示す説明図である。

【図10】実施例２における実測した応答の一例を示す説明図である。

【図11】実施例２における部分実測データの第一例を示す説明図である。

【図12】実施例２における部分実測データの第二例を示す説明図である。

【発明を実施するための形態】

【0012】

（本発明の基礎となった知見）
シミュレーションは、例えば、工業的な装置または機械が、その使用目的に適った性能を発揮できるか否かを事前に評価することを目的として、装置または機械の設計段階において実施される。そして、実施されたシミュレーションの結果に基づいて、装置または機械の設計の最適化がなされ得る。

【0013】

こうした工業的な装置または機械における物理現象を対象としたシミュレーションを精度よく実施するためには、解析プログラムに与える解析パラメータが適切である必要がある。解析パラメータは、例えば、材料の特性を表す物性値、初期条件または境界条件、もしくは、解析プログラムの中で物理現象を記述するために用いられる物理モデルが含有するモデルパラメータなどである。

【0014】

例えば、電子部品の温度分布を評価するために、熱伝導方程式の推定式などに基づいた伝熱シミュレーションが行われることがある。その場合、境界条件としての流入熱量、流出熱量、材料に関する熱伝導率、比熱、または、密度などの物性値を解析パラメータとして適切に与えない限り、実用に足る解析精度を得ることはできない。

【0015】

しかしながら、一般に実際の工業的な装置または機械を対象とする場合、これらの解析パラメータを事前に全て正確に知ることは困難であることが多い。

【0016】

例えば、熱伝導率または比熱などの熱物性値は比較的容易に計測されることが可能ではあるが、一方、複雑に組み合わされた部品間の熱抵抗、または、自然対流による熱伝達の影響などのように、計測されることが難しいパラメータも多く存在する。

【0017】

また、伝熱のシミュレーションに限らず、例えば粉体のシミュレーションにおいては、一般的に個別要素法などの手法が用いられるが、そこで必要となる解析パラメータとして、粒子間の摩擦係数、転がり抵抗、または、表面エネルギーなどがある。これらは、単体としては測定が非常に難しいパラメータであり、また、複数の粉体の組み合わせ、または、粉体と構造物との組み合わせにより変化するパラメータであるので、正確な値をシミュレーションに盛り込むことは容易でない。

【0018】

シミュレーションを実施する際に、事前に正確に知ることができない解析パラメータについては、何らかの仮定値が与えられる。しかしながら、仮定値は、十分な精度を有していないこともあり得る。そのため、シミュレーション自体を実施することが可能であっても、シミュレーションにより得られた結果が目的に適う十分な精度を有しているとは限らない。

【0019】

すなわち、装置または機器の設計に足る十分な精度において物理現象を模擬するためのシミュレーションを実施するためには、未知の解析パラメータの値を的確に推定する必要があるという課題がある。

【0020】

上記の課題に対して、近年、逆問題あるいは逆解析の有用性が注目されている。逆問題または逆解析は、出力から入力を、または結果から原因を求めようとする手法である。言い換えれば、逆問題または逆解析は、入力から出力を、または、原因から結果を直接的に求める順問題または順解析とは逆の概念である。

【0021】

電子部品の温度分布を評価する上記問題についても、逆解析手法を適用することにより、出力である部品各所の温度時刻歴の測定結果から、入力である入熱および放熱の境界条件、または、電子部品の熱物性値を求めることができると期待される。

【0022】

例えば、特許文献１では、複数の材質からなる被測定体の温度を計測することにより、測定が困難な熱物性を同定する手法について記述されている。具体的には、ホットプレート上で加熱された複数材料からなる部材に対して、任意の箇所の温度履歴を計測する。一方で、有限要素法または有限体積法などの計算手法からなるシミュレーションによって、熱物性値の関数である温度推定式であって、任意の時間における温度を示す温度推定式を作成する。その上で、温度推定式から得られる信号と実測温度との差異を示す評価関数を導出し、その評価関数を最小化する解を求めることによって物性値を同定する。なお、特許文献１において、温度推定式は、実験計画法の直交表に基づいて多数のシミュレーションを行い、応答曲面法などの手法を用いて作成される、との記述がある。

【0023】

特許文献１に記述されている概念または手法は、様々な工業的な装置または機械を対象としたシミュレーションにおいて、目的に適う十分な精度が得られるように、未知の解析パラメータの値を的確に推定するために適用することが可能と考えられるが、実用上は以下の問題がある。

【0024】

特許文献１においては、熱物性値が同定される部材を含む対象物の任意の箇所における温度履歴を計測し目的変数を設定している。具体的には、対象物の任意の箇所の任意の時刻における温度Ｔ^ｅｘを測定し、実測温度信号として出力する。例えば、時刻ｔ１、ｔ２、…、ｔｎにおいて計測した温度を、それぞれ、実測温度Ｔ１^ｅｘ、Ｔ２^ｅｘ、…、Ｔｎ^ｅｘとして、得られた各実測温度を目的変数として定義している。

【0025】

しかしながら、この場合、対象物の初期の温度状態次第では、時刻ｔ１、ｔ２、…、ｔｎそれぞれにおける実測温度Ｔ１^ｅｘ、Ｔ２^ｅｘ、…、Ｔｎ^ｅｘが一意に定まらない可能性がある。

【0026】

この場合、初期温度状態を実測し、実測した初期温度状態をシミュレーションの初期状態として用いて同定を進めることが考えられるが、温度の分布を正確にシミュレーションに盛り込むことは難しい。また、ここで挙げられている温度に関しては、その初期状態のモニタリングが比較的容易であるが、流体または粉体の挙動など、他分野においては、初期状態のモニタリングが困難な場合もあり、時刻と実測値とにより特定される二次元的な定点を目的変数として用いるだけでは、同定される物性値の精度を向上させることは難しい。

【0027】

また、実測で得られた実測データは、時系列データのような連続的なデータであり、モニタリング装置の設定分解能次第では、多数の二次元的な定点として得られる。そのため、目的変数の取り方によっては、目的変数が、その系における特徴的な傾向を示していない可能性がある。その結果、推定される物性値の精度が低下することがあり得る。推定される物性値の精度の低下を抑制するには、目的変数とする定点を適切に決定することが課題となる。

【0028】

以上の例とは別に、得られた連続的なデータに対して、近似曲線を作成し、その近似曲線を表す多項式の係数を目的変数とするという方法も考えられる。この方法では、多数の実測点のデータを、例えば複数の次数を持った多項式で近似することで、実測値を平滑化することができるので、上記課題を検討する必要はない。しかしながら、場合によっては近似曲線が複雑化し、目的変数の個数が物性値の個数よりも多くなってしまう場合がある。また、そもそも近似曲線に近似することが難しい場合もある。

【0029】

したがって、高い精度でパラメータを推定するには、実測データから、その系における特徴を表す多数の目的変数を取得する手法が必要であるが、そのような手法はこれまでに知られていない。

【0030】

そこで、本発明は、物性値などの解析パラメータに計測困難な値が含まれる場合であっても、実験によって時系列データを取得し、適切な目的変数を多数取得することにより、解析パラメータを推定することが可能となる生成方法等を提供する。

【0031】

例えば、本発明によれば、物理現象を模擬するシミュレーションを所望の解析精度において実施するために、未知の解析パラメータの値を的確に推定できる。

【0032】

以下、本明細書の開示内容から得られる発明を例示し、その発明から得られる効果等を説明する。

【0033】

（１）物理現象の条件を示す条件値から、前記物理現象から得られる観測値を推定する推定式を生成する生成方法であって、前記物理現象の時系列の観測値と、前記時系列の観測値の複数の時間区間それぞれにおける変化率を含む複数の第一変化率とを取得し、前記物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値であって、時系列の物性値と、前記時系列の物性値の複数の時間区間それぞれにおける変化率を含む複数の第二変化率とを少なくとも含む１以上の第二パラメータ値を取得し、前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて前記推定式を生成して出力する、生成方法。

【0034】

上記態様によれば、生成装置は、時系列の観測値と、複数の第一変化率と、１以上の第一パラメータ値と、１以上の第二パラメータ値とを用いて推定式を生成して出力する。このとき、生成装置は、時系列の観測値から得られる複数の第一変化率を、その物理現象の特徴を表す目的変数として用いることができるので、パラメータの推定の精度を高くすることができる。このように、上記生成方法は、複数の目的変数を用いて高い精度でパラメータを推定することができる。

【0035】

（２）前記複数の第一変化率を取得する際には、前記時系列の観測値の変化率が第一条件を満たす時刻、または、前記時系列の観測値の変化率の変化率が第二条件を満たす時刻で時間軸を区切ることで得られる複数の時間区間を、前記複数の時間区間として用いて前記複数の第一変化率を取得する、（１）に記載の生成方法。

【0036】

上記態様によれば、生成装置は、時系列の観測値の変化率、または、時系列の観測値の変化率の変化率に特徴がある時刻によって区切られた時間区間を用いて時系列の観測値から複数の第一変化率を算出する。時系列の観測値の変化率、または、時系列の観測値の変化率の変化率に特徴がある時刻の前後において、その物理現象が変化していることがある。言い換えれば、物理現象の変化によって、時系列の観測値の変化率、または、時系列の観測値の変化率の変化率の特徴が変化している可能性がある。そこで、時間区間ごとに算出される、その時間区間に生じている物理現象による影響を受けて生じた第一変化率を目的変数として用いることで、より高い精度でパラメータを推定することができる。このように、上記生成方法は、時間区間を用いて算出した第一変化率を目的変数として用いることで、より高い精度でパラメータを推定することができる。

【0037】

（３）前記複数の第一変化率を取得する際には、前記物理現象が生じている系を観測することで、前記系の状態が変化した１以上の時刻を取得し、取得した前記１以上の時刻で時間軸を区切ることで得られる複数の時間区間を、前記複数の時間区間として用いて前記複数の第一変化率を取得する、（１）に記載の生成方法。

【0038】

上記態様によれば、生成装置は、物理現象が生じている系の状態が変化したことが観測された時刻によって区切られた時間区間を用いて時系列の観測値から複数の第一変化率を算出する。物理現象が生じている系の状態が変化したことが観測された時刻の前後において、その物理現象が変化していることがある。そこで、時間区間ごとに算出される、その時間区間に生じている物理現象による影響を受けて生じた第一変化率を目的変数として用いることで、より高い精度でパラメータを推定することができる。よって、上記生成方法は、時間区間を用いて算出した複数の第一変化率を目的変数として用いることで、より高い精度でパラメータを推定することができる。

【0039】

（４）前記推定式を生成する際には、前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて応答曲面法により前記推定式を生成する、（１）～（３）のいずれかに記載の生成方法。

【0040】

上記態様によれば、推定式を生成するのに応答曲面法を用いるので、より容易に、推定式を生成することができる。よって、上記生成方法は、より容易に、複数の目的変数を用いて高い精度でパラメータを推定することができる。

【0041】

（５）物理現象の条件を示す条件値から、前記物理現象から得られる観測値を推定する推定式を生成する生成装置であって、前記物理現象の時系列の観測値と、前記時系列の観測値の複数の時間区間それぞれにおける変化率を含む複数の第一変化率とを取得するデータ取得部と、前記物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値であって、時系列の物性値と、前記時系列の物性値の複数の時間区間それぞれにおける変化率を含む複数の第二変化率とを少なくとも含む１以上の第二パラメータ値を取得するパラメータ取得部と、前記時系列の観測値と、前記複数の第一変化率と、前記１以上の第一パラメータ値と、前記１以上の第二パラメータ値とを用いて前記推定式を生成して出力する生成部とを備える、生成装置。

【0042】

上記態様によれば、上記生成方法と同様の効果を奏する。

【0043】

（６）（１）～（５）のいずれかに記載の生成方法をコンピュータに実行させるためのプログラム。

【0044】

上記態様によれば、上記生成方法と同様の効果を奏する。

【0045】

なお、これらの包括的または具体的な態様は、システム、装置、集積回路、コンピュータプログラムまたはコンピュータ読み取り可能なＣＤ－ＲＯＭなどの記録媒体で実現されてもよく、システム、装置、集積回路、コンピュータプログラムまたは記録媒体の任意な組み合わせで実現されてもよい。

【0046】

以下、実施の形態について、図面を参照しながら具体的に説明する。

【0047】

なお、以下で説明する実施の形態は、いずれも包括的または具体的な例を示すものである。以下の実施の形態で示される数値、形状、材料、構成要素、構成要素の配置位置及び接続形態、ステップ、ステップの順序などは、一例であり、本発明を限定する主旨ではない。また、以下の実施の形態における構成要素のうち、最上位概念を示す独立請求項に記載されていない構成要素については、任意の構成要素として説明される。

【0048】

（実施の形態）
本実施の形態において、複数の目的変数を用いて高い精度でパラメータを推定する推定式を生成する生成装置および生成方法などについて説明する。

【0049】

図１は、本実施の形態における生成装置１０の構成を示すブロック図である。

【0050】

生成装置１０は、実測で得られた時系列データから、物性値などのパラメータを同定するために用いる目的変数を生成し、物理現象の条件を示す条件値から、その物理現象から得られる観測値を推定する推定式を生成する。

【0051】

図１に示されるように、生成装置１０は、機能部として、データ取得部１１と、パラメータ取得部１４と、推定式生成部１５を備える。生成装置１０が備える機能部は、生成装置１０が備えるプロセッサ（例えばＣＰＵ（Ｃｅｎｔｒａｌ Ｐｒｏｃｅｓｓｉｎｇ Ｕｎｉｔ））（不図示）が、メモリ（不図示）を用いてプログラムを実行することで実現され得る。

【0052】

データ取得部１１は、実験において、温度センサまたはロードセルなどによる測定器によって出力された時系列の観測値を含む、時系列データを取得する。

【0053】

データ取得部１１は、区間決定部１２と、目的変数生成部１３とを有する。

【0054】

区間決定部１２は、データ取得部１１が取得した時系列データにおいて複数の時間区間を決定する。複数の時間区間は、複数の目的変数を生成するためのデータ範囲として用いられる。

【0055】

目的変数生成部１３は、データ取得部１１が取得した時系列データのうち、区間決定部１２が決定した複数の時間区間それぞれに含まれる観測値の時間に対する変化率（第一変化率に相当）を算出する。なお、「時間に対する変化率」を単に「変化率」ともいう。以降においても同様である。

【0056】

観測値の変化率の算出は、当該観測値と、当該観測値が観測された時刻と異なる時刻に観測された観測値との差分を、その時刻間の差分で除算することでなされ得る。当該観測値が観測された時刻と異なる時刻に観測された観測値は、例えば、当該観測値に時間的に隣接する観測値、言い換えれば、当該観測値より１つ前、または、１つ後に観測された観測値であり得るが、これに限られない。上記変化率は、目的変数の１つとして用いられる。上記変化率を算出することは、目的変数を生成することに相当する。目的変数生成部１３は、算出した変化率を目的変数としてメモリに格納する。

【0057】

パラメータ取得部１４は、物理現象を模擬したシミュレーションを、１以上の第一パラメータ値それぞれを条件値として実施した結果を示す１以上の第二パラメータ値を取得する。１以上の第二パラメータ値は、時系列の物性値と、時系列の物性値の複数の時間区間それぞれにおける変化率とを含む。

【0058】

なお、シミュレーションを実施する主体は、生成装置１０であってもよいし、生成装置１０と異なる装置（シミュレータ装置ともいう）（不図示）であってもよい。

【0059】

シミュレーションを実施する主体がシミュレータ装置である場合、パラメータ取得部１４は、複数の第一パラメータ値を取得する。また、パラメータ取得部１４は、複数の第一パラメータ値をシミュレータ装置に送信し、送信した複数の第一パラメータ値を用いてシミュレータ装置が実施したシミュレーションの結果である第二パラメータ値を受信することで、第二パラメータ値を取得する。

【0060】

推定式生成部１５は、推定式を生成して出力する。推定式生成部１５は、推定式を生成する際には、データ取得部１１が取得した時系列データ（つまり時系列の観測値）と、目的変数生成部１３が算出した観測値の変化率と、パラメータ取得部１４が取得した１以上の第一パラメータ値と、１以上の第二パラメータ値とを用いて、推定式を生成する。推定式生成部１５は、推定式を生成する際には、時系列の観測値と、複数の第一変化率と、１以上の第一パラメータ値と記１以上の第二パラメータ値とを用いて応答曲面法により推定式を生成することができる。

【0061】

目的変数生成部１３は、複数の第一変化率を取得する際には、時系列の観測値の変化率が第一条件を満たす時刻、または、時系列の観測値の変化率の変化率が第二条件を満たす時刻で時間軸を区切ることで得られる複数の時間区間を用いて複数の第一変化率を取得することができる。ここで、時系列の観測値の変化率が第一条件を満たす時刻は、例えば、時系列の観測値の変化率がゼロである時刻、または、極大値もしくは極小値である時刻である。また、時系列の観測値の変化率の変化率が第二条件を満たす時刻は、例えば、時系列の観測値の変化率の変化率がゼロである時刻、または、極大値もしくは極小値である時刻である。

【0062】

また、目的変数生成部１３は、複数の第一変化率を取得する際には、物理現象が生じている系を観測することで、系の状態が変化した１以上の時刻を取得し、取得した１以上の時刻で時間軸を区切ることで得られる複数の時間区間を用いて複数の第一変化率を取得することもできる。系の状態は、例えば、系としてのワークの表面の溶融状態等であり得る。ワークの表面の溶融状態は、ワークの表面の温度の計測により定められる。

【0063】

以上のように構成された生成装置１０が実行する処理（生成方法ともいう）を以下で説明する。

【0064】

図２は、本実施の形態における生成方法を示すフロー図である。

【0065】

ステップＳ１０１において、データ取得部１１は、実験にて得られる時系列の実測データである観測値を取得する。実測データは、例えば、熱伝導に関する実験であれば温度であり、流体に関する実験であれば圧力または流速であり、粉体に関する実験であれば粉体の嵩密度または流動性に関する力などである。

【0066】

ステップＳ１０２において、区間決定部１２は、ステップＳ１０１にて取得した実測データを分析し、複数の時間区間を決定する。区間決定部１２が決定する複数の時間区間は、その複数の時間区間を用いて目的変数生成部１３が目的変数を生成するのに適した時間区間であり得る。

【0067】

ステップＳ１０３において、目的変数生成部１３は、観測値の変化率を算出する。具体的には、目的変数生成部１３は、ステップＳ１０２で区間決定部１２が決定した時間区間を用いて、ステップＳ１０１でデータ取得部１１が取得した実測データを分割する。分割により、データ取得部１１が取得した実測データから、各時間区間に含まれる実測データ（部分実測データともいう）が生成される。そして、目的変数生成部１３は、部分実測データのそれぞれについて、観測値の変化率を算出する。目的変数生成部１３は、算出した変化率を、後述の実験計画における目的変数としてメモリに格納する。

【0068】

ステップＳ１１１において、パラメータ取得部１４は、シミュレーションの条件値として用意された実験計画を取得する。実験計画は、例えば、不明な解析パラメータである物性値または特性値（条件値に相当）を含む初期実験計画モデルであり、実験計画法によって作成されたものであり得る。シミュレーションの実施に必要な１以上の条件値の組を条件とも言う。実験計画には、複数の条件が含まれ得る。

【0069】

実験計画法（Ｄｅｓｉｇｎ ｏｆ Ｅｘｐｅｒｉｍｅｎｔｓ：ＤＯＥ）には、古典的な直交計画法、中心複合計画法およびＳｐａｃｅＦｉｌｌｉｎｇ計画法というように、目的に応じて多種の方法がある。直交計画法は、交互作用が強く働く場合に精度が低い傾向があり、また、ＳｐａｃｅＦｉｌｌｉｎｇ計画法は、実験回数が多い傾向がある。複雑なプロセスにおいては、一般に、交互作用が強く働き、実験またはシミュレーションに多くの時間を要することが想定される。そのような場合には、初期実験計画モデルの作成は、例えば中心複合計画法を用いて行われる。以降では、中心複合計画法を用いて初期実験計画モデルを作成する場合を例として説明する。

【0070】

ステップＳ１１２において、パラメータ取得部１４は、ステップＳ１１１で取得した初期実験計画モデルに基づいて、物性値または特性値からなる解析パラメータを条件値として設定し、実験条件の全ての組合せが完了するまでシミュレーション（初期シミュレーションともいう）を実施する。そして、パラメータ取得部１４は、初期シミュレーションの結果（つまり応答）を取得する。例えば、伝熱に関するシミュレーションの場合、応答は、材料における各部の温度などである。

【0071】

ステップＳ１１３において、推定式生成部１５は、パラメータ取得部１４にてステップＳ１１２で得られた応答に対して応答曲面法を適用し、応答曲面法に基づいて、解析パラメータから得られる応答を予測する推定式を生成する。また、推定式生成部１５は、生成した推定式と解析パラメータとを用いて、応答を推定する。

【0072】

また、パラメータ取得部１４は、生成した推定式を用いて、実際の物理現象から得られた観測値が、シミュレーションの結果として得られるような解析パラメータを導出する。解析パラメータの解として、複数の解が存在する場合には、例えば、全ての解析パラメータに対して数万から数百万条件のランダムサンプリングを行い、観測値に比較的近い推定値を得ることができる条件を複数抽出する。

【0073】

ステップＳ１１４において、パラメータ取得部１４は、ステップＳ１１３で解として得られた複数の解析パラメータの中から推定式の精度を確認するためのシミュレーション（確認シミュレーションともいう）に用いる１以上の条件を抽出する。条件の抽出の方法に関しては、同様の解析パラメータをできるだけ選ばないよう、クラスタリングしてから抽出することが望ましい。クラスタリングの手法は多々あるが、例えばＫｍｅａｎｓ法などを用いる。

【0074】

ステップＳ１１５において、パラメータ取得部１４は、ステップＳ１１４で抽出した複数の条件を用いて、確認シミュレーションを実施する。

【0075】

ステップＳ１１６において、パラメータ取得部１４は、ステップＳ１１５で得た確認シミュレーションの結果と、ステップＳ１１３で推定式から得た推定値とを比較することで、推定式の精度を算定する。また、パラメータ取得部１４は、算定した推定式の精度が所定より高いか否かを判定する。ここで、推定式の精度が所定より高い、とは、推定式が物理現象の条件から観測値を推定する推定式として十分な精度を有していること、言い換えれば、推定式が物理現象の条件から観測値を推定する推定式として実用に足ることを意味し得る。推定式の精度が所定より高いと判定した場合（ステップＳ１１６でＹｅｓ）、ステップＳ１１７に進み、そうでない場合（ステップＳ１１６でＮｏ）、ステップＳ１２１へ進む。なお、推定式の精度は、推定式により推定された推定値とシミュレーションの結果を示す値との誤差、または、ＲＭＳＥ（Ｒｏｏｔ Ｍｅａｎ Ｓｑｕａｒｅｄ Ｅｒｒｏｒ）などの評価指標を用いて評価され得る。

【0076】

ステップＳ１２１において、パラメータ取得部１４は、推定式の精度を上げるために、実験計画の修正を行う。例えば、パラメータ取得部１４は、確認シミュレーションの結果を推定式の生成に用いる応答に含めることで、実験計画の修正を行う。また、パラメータ取得部１４は、説明変数を増やすように実験計画を拡張することにより、新たな実験計画を生成することで、実験計画の修正を行ってもよい。

【0077】

ステップＳ１２２において、パラメータ取得部１４は、追加シミュレーションが必要であるか否かを判定する。パラメータ取得部１４は、例えば、ステップＳ１２１での実験計画の修正において、説明因子を拡張した場合には、追加シミュレーションが必要であると判定することができる。一方、ステップＳ１２１での実験計画の修正において、確認シミュレーションの応答を追加するのみの場合には、追加シミュレーションが必要でないと判定することができる。追加シミュレーションが必要であると判定した場合（ステップＳ１２２でＹｅｓ）には、ステップＳ１１２に進み、そうでない場合（ステップＳ１２２でＮｏ）には、ステップＳ１１３に進む。ステップＳ１１２に進む場合、ステップＳ１２１で修正された新たな実験計画における応答をシミュレーションで取得する。

【0078】

ステップＳ１１７において、推定式生成部１５は、推定式を出力する。目的変数生成部１３が出力する推定式は、ステップＳ１１３で推定式生成部１５が生成した推定式であって、ステップＳ１１６でパラメータ取得部１４により精度が所定より高いと判定された推定式である。

【0079】

図２に示される一連の処理により、生成装置１０は、複数の目的変数を用いて高い精度でパラメータを推定する推定式を生成することができる。

【0080】

（実施例１）
以下に上記実施の形態の具体的な実施例を示す。本実施例では、粉体の圧縮工程を事例として記述するが、開示技術がこれにより限定されるものではない。

【0081】

圧縮装置の１つであるロールプレス装置は、工業用途で幅広く利用され、フィルム、紙、不織布、金属箔または鋼板などの材料を圧延成形する目的で用いられている。ロールプレス装置は、その他にも、シート表面の艶出し、張り合わせまたは繊維材料の絞り（脱水）など、成形以外の用途に適用されることも多い。また、単純な板材の圧延以外にも、ステンレス、鉄、ニッケル、アルミニウムなどの金属の粉体をロールプレス装置に連続的に供給することで成形する方法もよく知られている。

【0082】

このロールプレス装置を用いた粉体圧縮プロセスに関して、個別要素法による粉体のシミュレーションを適用するためには、粉体にまつわる物性値（具体的には、粒子間の摩擦係数、もしくは、粒子と構造物材料との間の摩擦係数、転がり抵抗、または、粒子のヤング率もしくは真密度など多数のパラメータ）が必要である。

【0083】

本発明を、このロールプレス装置を用いた粉体圧縮プロセスに対して適用するのも可能ではあるが、プロセスが大規模であり、多数のシミュレーションを実施するには効率的ではないので、問題を簡素化し、より簡易的な実験において物性値を同定する方法をとる。

【0084】

例えば、少量の粉体をサンプルとして用いた圧縮せん断試験で物性値を同定する方法がある。上記方法について、以下で図３Ａ、図３Ｂおよび図３Ｃを参照しながら説明する。

【0085】

図３Ａ、図３Ｂおよび図３Ｃは、粉体の圧縮せん断試験を示す説明図である。

【0086】

図３Ａに示されるように、圧縮せん断試験機は、杵２１と臼２２と底板２３とを有する円筒状のセルを備える。圧縮せん断試験機は、杵２１が垂直に押し出されることでセル内に充填させた粉体２４が圧縮される過程（圧縮過程ともいう）を模擬することができる（図３Ｂ参照）。また、圧縮せん断試験機は、底板２３が水平に押し出されることでセル内に充填させた粉体２４がせん断される過程（せん断過程ともいう）を模擬することができる（図３Ｃ参照）。

【0087】

杵２１、および、底板２３には、ロードセルが設置されており、それぞれにかかる荷重が測定できるようになっている。

【0088】

まず、杵２１が一定速度で粉体２４を圧縮していく（図３Ｂ参照）。このとき、杵２１にかかる反力がロードセルによって計測され、結果として反力の時系列のデータが取得される。反力が、指定された圧縮荷重に達した時点で、その荷重が維持された状態で底板２３が一定速度で水平に移動を始める（図３Ｃ参照）。このとき、底板２３に設置したロードセルによって、せん断力の時系列データが取得される。

【0089】

図４は、シミュレーションにおけるパラメータと推定式との関係を示す説明図である。

【0090】

図４に示されるように、粉体の圧縮せん断試験の入力は、粉体材料の粉体物性と、壁面物性と、粉体－壁面間特性とである。また、粉体の圧縮せん断試験の出力は、粉体材料の圧縮特性およびせん断特性である。

【0091】

ここで、入力である粉体物性、壁面物性および粉体－壁面間特性と、出力である圧縮特性およびせん断特性との関係を示す関係式を以下のように設定する。

【0092】

第１関係式：圧縮特性＝ｆ１（粉体物性，壁面物性，粉体－壁面間特性）
第２関係式：せん断特性＝ｆ２（粉体物性，壁面物性，粉体－壁面間特性）

【0093】

ここで、粉体物性と壁面物性とは、それぞれの材料のヤング率、ポアソン比、粒子同士の摩擦係数または転がり抵抗などを示す。粉体－壁面間特性は、粒子と壁面との間の摩擦係数または転がり抵抗などを示す。

【0094】

圧縮特性とせん断特性とは、それぞれロードセルによって測定された力のデータである。圧縮特性とせん断特性とは、粉体物性と、壁面物性と、粉体－壁面間特性との関数として表現される。

【0095】

圧縮特性とせん断特性とを示すパラメータを推定すべく、実験計画に基づくシミュレーション結果から、この関数と同等の意義を有する推定式を精度よく求めることが目的である。

【0096】

本来であれば、求めたい物性値または特性値の数に合わせて、プロセス条件を変えた場合の実測値または別の観測点での実測値を収集し、それぞれの実測値に対する応答をシミュレーションにより取得し、得られた複数の推定式から物性値を一意に決めることが望ましい。ここでは簡単のため、粉体の、ある圧縮量に対する圧縮力を応答の一例として、図２のフロー図に則って推定式を生成する処理を説明する。

【0097】

図５は、本実施例における実測した応答の一例を示す説明図である。図５には、圧縮せん断試験により得られた、圧縮特性の一例である圧縮荷重の実測データがグラフとして示されている。図５に示される実測データは、ステップＳ１０１にてデータ取得部１１が取得した実測データに相当する。

【0098】

図５のグラフにおいて、横軸が圧縮量［ｍｍ］であり、縦軸が圧縮力［Ｎ］である。なお、杵２１が一定速度で粉体２４を圧縮するので、圧縮量［ｍｍ］は、時間［ｓｅｃ］に対して線形であり、言い換えれば、図５のグラフの横軸を時間［ｓｅｃ］で表現することもできる。

【0099】

図５のグラフは、前述の粉体２４の圧縮荷重を表している。グラフに示される曲線は、理論的な式、または、近似式に近似することが困難な曲線であり、粉体挙動の複雑性が顕著に表されている。

【0100】

図６は、本実施例における系の状態の例を示す説明図である。図６には、図５に示した圧縮荷重において、粉体２４の３つの状態Ａ、ＢおよびＣを示した模式図が示されている。図５の曲線が表す、粉体２４の複雑な挙動の概要は次のように推測される。

【0101】

図６の（ａ１）、（ｂ１）および（ｃ１）はそれぞれ、状態Ａ、ＢおよびＣの各状態における粉体２４に含まれる粒子の様子を示しており、図６の（ａ２）、（ｂ２）および（ｃ２）はそれぞれ、上記各状態における粉体２４の体積を模式的に示している。

【0102】

状態Ａ（つまり（ａ１）および（ａ２））は、圧縮を開始した当初における粉体２４の状態である。圧縮を開始した当初において、粉体２４は、初期状態であり、空隙を多く含んだ状態である。

【0103】

状態Ｂ（つまり（ｂ１）および（ｂ２））は、状態Ａから圧縮を進めている途中における粉体２４の状態が示されている。状態Ａから、杵２１による圧縮を進めると、粉体２４に含まれる粒子間の空隙が徐々に埋まっていくことで、粒子の再配置が生じ、次第に最密充填状態に近づいていく。

【0104】

状態Ｃ（つまり（ｃ１）および（ｃ２））は、状態Ｂからさらに圧縮を進めたときの粉体２４が示されている。状態Ｂからさらに圧縮を進めていくと、粉体２４は、粉体２４に含まれる粒子間の空隙がほぼ埋まりきり、杵２１により圧縮されたときに粒子自体が圧縮される状態となる。

【0105】

言い換えれば、粉体２４が圧縮されたときに生ずる現象は、比較的初期段階では、粒子の再配置が支配的であり（例えば状態Ａ～状態Ｂ）、圧縮が進むにつれて粒子自体の圧縮が支配的となる（例えば状態Ｂ～状態Ｃ）。粉体２４が圧縮されたときに生ずる現象に差異があれば、圧縮により変化する粉体２４の状態または特性にも差異がある。

【0106】

このような粉体２４の状態変化が、図５に示した圧縮荷重の曲線の曲率に反映されていると推測される。すなわち、粉体２４が圧縮されるにつれて、粉体２４の状態が、状態Ａ、ＢおよびＣと変化し、その状態変化を反映して圧縮荷重の曲線の曲率が徐々に変化していくと推測される。言い換えれば、圧縮荷重の曲線は、粉体２４の特徴的な状態を全て含めて表しているともいえる。

【0107】

前述の通り、圧縮荷重の曲線を理論的な式、または、近似式に近似することが困難である。また、粉体２４の初期状態はシミュレーション上で再現しづらく、２次元的な定点を目的変数として定義することが困難である。そのため、実測したデータにおいて、何を目的変数として定義するかを検討する必要がある。

【0108】

本発明においては、実測データを分割し、分割した範囲に含まれる部分実測データである観測値の変化率を目的変数として定義する。ここでは図６における３つの状態Ａ、ＢおよびＣを想定して、実測データを３つに分割する。実測データの分割範囲は、区間決定部１２が決定する、圧縮量の区間により定められる（ステップＳ１０２）。なお、上記の通り、圧縮量が時間に対して線形であるので、区間決定部１２は、時間区間を決定することで、圧縮量の区間を決定する。

【0109】

図７は、本実施例における区間の一例を示す説明図である。図７には、図５に示される圧縮荷重の変化率の変化率がグラフとして示されている。図７のグラフは、横軸が圧縮量［ｍｍ］であり、縦軸が、圧縮荷重の変化率の変化率である。圧縮荷重の変化率の変化率は、図５に示される圧縮荷重のグラフにおける傾きの増減として表されている。

【0110】

図５に示される圧縮荷重の曲率の増加幅が大きい圧縮量ほど、当該圧縮量における、圧縮荷重の変化率の変化率が大きい。このような圧縮量は、粉体２４に含まれる粒子の様子に特徴的な変化が起こった点であると推定される。そこで、区間決定部１２は、圧縮荷重の変化率の変化率が極大値をとる圧縮量として、例えば図７に示される圧縮量Ｖ１、Ｖ２およびＶ３を、圧縮量の区間の境界として実測データを分割することができる。

【0111】

図８は、本実施例における部分実測データの一例を示す説明図である。図８には、ステップＳ１０２で区間決定部１２が決定した区間が、圧縮荷重の曲線（図５に示されるものと同じ）上に示されている。

【0112】

具体的には、区間決定部１２は、時系列データのうち、圧縮量Ｖ１から圧縮量Ｖ２までの区間に含まれる時系列データを、部分実測データＤ１とする。区間決定部１２は、時系列データのうち、圧縮量Ｖ２から圧縮量Ｖ３までの区間に含まれる時系列データを、部分実測データＤ２とする。区間決定部１２は、時系列データのうち、圧縮量Ｖ３以上の区間に含まれる時系列データを、部分実測データＤ３とする。

【0113】

目的変数生成部１３は、部分実測データＤ１、Ｄ２およびＤ３のそれぞれについて、圧縮量に対する圧縮荷重の傾きを算出し、取得する。なお、各部分実測データについての上記傾きは、各部分実測データに含まれる両端の点だけから求められてもよいし、各部分実測データに含まれる任意の一部の点、または、全部の点から求められてもよい。

【0114】

目的変数生成部１３は、例えば、部分実測データＤ２についての圧縮荷重の傾きとして、部分実測データＤ２に含まれる両端の点Ｐ２およびＰ３を結ぶ直線Ｌの傾きを算出することができる。

【0115】

なお、目的変数生成部１３が各部分実測データに含まれる一部または全部の点から、各部分実測データについての傾きを算出する例を、図９を参照しながら説明する。

【0116】

図９は、本実施例における部分実測データの傾きの別の例を示す説明図である。図９には、図８における部分実測データＤ２に含まれる観測値を示す点Ｐ２、Ｐ２１、Ｐ２２、Ｐ２３、Ｐ２４、Ｐ２５、・・・、Ｐ３が示されている。なお、点Ｐ２５と点Ｐ３との間の点の図示は省略されている。

【0117】

目的変数生成部１３は、部分実測データＤ２についての圧縮荷重の傾きとして、部分実測データＤ２に含まれる一部の点としての例えば点Ｐ２、Ｐ２２、Ｐ２４、・・・の傾き（つまり、直線Ｍ１、Ｍ３、Ｍ５、・・・の傾き）の平均値を算出することができる。なお、一部の点は、上記に限られない。

【0118】

また、目的変数生成部１３は、部分実測データＤ２についての圧縮荷重の傾きとして、部分実測データＤ２に含まれる全部の点Ｐ２、Ｐ２１、Ｐ２２、・・・の傾き（つまり、全部の直線Ｍ１、Ｍ２、Ｍ３、・・・の傾き）の平均値を算出することができる。

【0119】

このように算出された傾きが、目的変数として実験計画に含められ（ステップＳ１１１）、実験計画に基づいて推定式の生成および出力がなされる（Ｓ１１２～Ｓ１２２）。

【0120】

（実施例２）
次に、別の実施例について説明する。本実施例では、レーザを用いた溶接工程を事例として記述するが、開示技術がこれにより限定されるものではない。

【0121】

図１０は、溶接装置によるレーザを用いた溶接工程における実測データである。横軸はレーザ照射時間［ｓｅｃ］、縦軸はワークの表面温度［Ｋ］である。ワークの表面温度は、レーザが照射されているワークの温度状態を示しており、ひいては、ワークの表面の溶融状態を示している。

【0122】

図１０には、ワークの相状態の変化点、つまりワークの沸点および融点が破線により示されている。図１０に示されるように、融点より低い温度における相状態（つまり固体状態）でのワークの温度の変化率と、融点と沸点との間の温度における相状態（つまり液体状態）でのワークの温度の変化率とは、大きく異なっている。また、それぞれの相状態の中でも温度の変化率（つまり、曲線の傾きまたは曲率）が複雑に変化している。そのため、レーザ照射によるワークの状態変化について、伝熱および相変化を含むシミュレーションを精度よく実行するには、物性値を合わせこむことが重要である。

【0123】

本発明によれば、このような実測データに関しても、目的変数を取得し、シミュレーションにおける入力となる物性値を決定することが可能である。

【0124】

図１１は、本実施例における部分実測データの第一例を示す説明図である。図１１には、部分実測データＤ５およびＤ６が示されている。

【0125】

図１１に示される部分実測データＤ５およびＤ６は、ワークの温度が融点に等しい時間Ｔ１を境界として、ワークの温度が融点より低い時間区間と、ワークの温度が融点より高い時間区間とに分割されている。これは、融点が、ワークの状態変化の分岐点として化学的に定義されていることを根拠としている。

【0126】

なお、複合材などで構成されるワークにおいて、状態変化の分岐点が不明である場合、状態変化の分岐点は、ワークの溶融状態をリアルタイムに観測することで定義されてもよい。

【0127】

図１２は、本実施例における部分実測データの第二例を示す説明図である。

【0128】

前述の通り、図１０に示される表面温度の曲線は、融点前後における温度の変化率が複雑であり、状態変化に依存した伝熱パラメータ（具体的には熱伝導率等）の変化をシミュレーションに盛り込むためには、電熱パラメータの変化を正確に模擬することも必要となる。その場合、図１１に示される部分実測データＤ５およびＤ６では、物性値の同定のための目的変数が不足する可能性がある。

【0129】

その場合、図１１に示される時間Ｔ１に加えて、実施例１で示したように、実測データの変化率の変化率を用いて時間区間を決定することができる。図１２には、このように決定される時間区間の境界の例として、時刻Ｔ２、Ｔ３およびＴ４が示されている。

【0130】

区間決定部１２は、このように時間区間を増やすことで、部分実測データＤ５１、Ｄ５２、Ｄ６１、Ｄ６２およびＤ６３を生成し、これにより、目的変数を増やすことができる。これにより、シミュレーションにおける状態変化に依存した伝熱パラメータの変化を、より正確にすることに寄与し得る。

【0131】

以上、一つまたは複数の態様に係る生成装置などについて、実施の形態に基づいて説明したが、本発明は、この実施の形態に限定されるものではない。本発明の趣旨を逸脱しない限り、当業者が思いつく各種変形を本実施の形態に施したものや、異なる実施の形態における構成要素を組み合わせて構築される形態も、一つまたは複数の態様の範囲内に含まれてもよい。

【産業上の利用可能性】

【0132】

本発明に係る生成方法は、物性値などの解析パラメータが不明な場合において、解析パラメータと実測値の関係式を生成する生成装置に利用可能である。

【符号の説明】

【0133】

１０ 生成装置
１１ データ取得部
１２ 区間決定部
１３ 目的変数生成部
１４ パラメータ取得部
１５ 推定式生成部
２１ 杵
２２ 臼
２３ 底板
２４ 粉体
Ｄ１、Ｄ２、Ｄ３、Ｄ５、Ｄ５１、Ｄ５２、Ｄ６、Ｄ６１、Ｄ６２、Ｄ６３ 部分実測データ
Ｐ２、Ｐ２１、Ｐ２２、Ｐ２３、Ｐ２４、Ｐ２５、Ｐ３ 点
Ｌ、Ｍ１、Ｍ２、Ｍ３、Ｍ４、Ｍ５ 直線

【図1】

【図2】

【図3A】

【図3B】

【図3C】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11】

【図12】