Surrogate-based Shape Optimization of Immersion Nozzle in Continuous Casting

Abstract

In continuous casting, molten steel is fed from the tundish into the mold through the immersion nozzle. In the immersion nozzle, inclusions mainly composed of alumina present in the molten steel adhere and accumulate, it causes limitation of continuous castings. To prevent the nozzle clogging, Ar gas is blown into the immersion nozzle. However, Ar bubbles flow into the mold along with the molten steel and become trapped in the solidifying shell, causing bubbling defects of the slab. To suppress bubbling defects, it is effective to keep Ar bubbles away from the solidification interface or to use molten steel to wash away Ar bubbles that have adhered to the solidification interface. The molten steel flow in the mold is greatly affected by the shape of the immersion nozzle. In this paper, we consider the optimization of the shape of the immersion nozzle to reduce Ar bubbles trapped in the solidifying shell. A numerical model of molten steel flow and heat transfer solidification in the mold is combined with an optimization method. In the optimization process, Ar bubbles trapped in the solidifying shell are evaluated by a neural network to improve the calculation speed. The application of this method to the search for immersion nozzle shape is also reported, and the effectiveness of the obtained nozzle shape in reducing Ar bubbles is discussed.

1. 緒言

連続鋳造では，タンディッシュから浸漬ノズルを介して鋳型内に溶鋼が供給される。浸漬ノズルにおいては，溶鋼中に存在するアルミナ等を主成分とする介在物が付着堆積し，連続鋳造の連続数を制限する原因となる。このノズル閉塞を抑制するために，浸漬ノズルにArガスが吹き込まれる。Arガスを吹き込むことにより，連々数を確保できるが，Ar気泡が溶鋼とともに鋳型内に流入し，鋳片に捕捉されると気泡性欠陥の原因となる。この気泡性欠陥は，めっき不良の原因となるため，特に自動車用外装材に用いられるIF鋼で問題となる。鋳片の表層近くに捕捉されたAr気泡は，鋳片においてピンホール欠陥と呼ばれ，圧延後は引き延ばされたスジ疵となる。また，鋳片の表皮下数 mmの位置に捕捉されたAr気泡は，圧延後の鋼板表面に膨らみが生じる欠陥(フクレ疵)となる。フクレ疵の原因となる表皮下の気泡性欠陥を鋳片のブローホール欠陥と呼ぶ。これらの気泡性欠陥を防ぐためにはAr気泡を凝固界面に近づけないか，凝固界面に付着したAr気泡を溶鋼により洗い流すことが有効である。よって，鋳型内溶鋼流動の制御が重要となる。

これまでに，溶鋼流動を制御し凝固シェルへのAr気泡の捕捉を防止する技術が開発されている。電磁力を利用した電磁攪拌(Electro-Magnetic Stirring; EMS)や電磁ブレーキ(Electromagnetic Brake; EMBr)がその代表例である。鋳片表層の気泡性欠陥を抑制するには，EMSを鋳型上部に設置し，湯面近傍を攪拌することが効果的である。一方，EMBrは鋳型中央部および浸漬ノズル先端付近に設置されており，下部に向かう溶鋼流動に制動を加える装置である。特に浸漬ノズルの吐出流の下部に向かう流れを抑制することで，介在物やAr気泡を浮上分離させやすくなる。

鋳型内の溶鋼流動は浸漬ノズルの形状にも大きく影響するため，浸漬ノズルについてもこれまで検討がなされてきた^1–4)。しかしながら，ノズル詰まり，耐火物強度やコストなどの総合的なバランスから，スラブ鋳造においては，左右に1対の吐出孔を有する伝統的な2孔ノズルが用いられることが多く，近年は議論されることも少なくなった。一方，近年，鋳造から熱延までを一貫処理する薄スラブ連続鋳造機の導入が増加している。薄スラブ連鋳機では，引き抜き速度が高速となるため，浸漬ノズルの形状と鋳型内流動に関する研究が報告されている^5–9)。これらの報告を鑑みるにあたり，従来のスラブ連続鋳造機の浸漬ノズルに関しても未だ改良の余地は残されていると考えられる。そこで本稿では，従来の連続鋳造機において，凝固シェルへ捕捉されるAr気泡を低減する浸漬ノズルの形状を求めることを考える。特に，鋳片の表層~2 mmに捕捉されたAr気泡をピンホール，表皮下2 mm~15 mmに捕捉されたAr気泡をブローホールと呼び，これらを低減することを目的とする。

ピンホールとブローホールを算出するにあたり，これまでに開発，検討が重ねられてきた鋳型内の溶鋼流動および伝熱凝固の数値解析モデル^10,11)を用いる。本研究では，数値解析モデルと機械学習を使用した最適化技術の枠組みを合わせることにより，ピンホールとブローホールの最小化を目的とした浸漬ノズル形状の最適化を行う。このような数値流体力学(Computational Fluid Dynamics; CFD)と最適化技術を組み合わせた形状最適化技術は，航空機の翼設計¹²⁾や船の設計¹³⁾，高層ビルなどの構造物の設計¹⁴⁾において技術開発および適用例がみられるが，製鉄プロセスには見られない。CFDを用いた形状最適化における課題は，時間のコストが高く実用性に欠けることである。一般に最適化の過程においては数千回~数万回のシミュレーションを必要とするが，CFD計算は多くの計算時間を要するため，直接的に最適化技術を適用することは難しい。計算コストの削減と最適化の効率を図るために，近年では，代理モデルの研究が活発に行われている。代理モデルは，最適化過程において，複雑で時間のかかる計算を代替する近似モデルの一種である。クリギング¹⁴⁾や動径基底関数¹⁵⁾，(深層)ニューラルネットワーク¹⁶⁾などの機械学習モデルの適用事例が報告されている¹⁷⁾。

本稿では，ピンホールとブローホールは，溶鋼流動および伝熱凝固による複雑な現象の結果得られることを考慮し，CFD計算の代わりにニューラルネットワークモデルでピンホールとブローホールを求めることとする。すなわち，実施可能な回数のCFD計算で得られたピンホールとブローホールの評価値とその計算条件をデータセットとしてニューラルネットワークを学習し，学習済みのモデルで所与のノズル条件におけるピンホールとブローホールを求める。これにより理想とする最小のピンホールとブローホールとの誤差を求めることができるため，この誤差を目的関数として最適化アルゴリズムを適用して最適なノズル形状を求める。ニューラルネットワークによる予測は数秒程度であるため，高速な最適化が可能である。

最適化の枠組みには，反復的に最適化を行いながら最適解を求める逐次近似最適化^18,19)を適用する。本枠組みでは，反復過程で得られる計算条件とピンホールおよびブローホールからなるデータを段階的にデータセットに追加し，ニューラルネットワークモデルを更新して最適解を求めていく。ニューラルネットワークの予測精度はデータセットに強く依存する。データ数が少ない場合は予測精度が十分でない可能性があり，CFD計算の結果と乖離するため，1度の最適化の実行で最適解が得られる保証はない。一方，逐次近似最適化では，データセットが拡充されることでニューラルネットワークの予測精度が向上されるため，反復的に最適解を更新でき，ピンホールとブローホールを低減することができる。

最適化アルゴリズムには粒子群最適化²⁰⁾を使用する。これは，ランダムに位置する複数の粒子(探索点)が，決められた規則に則って探索空間を移動しながら最適解を探していく方法である。局所最適解に陥りにくく多峰性関数のような複雑な目的関数に対しても大域最適解を求める能力が高いことが特徴である。

以上に加えて，動径基底関数ネットワークによる密度関数を用いたデータ点の追加方法²¹⁾も取り入れる。逐次近似最適化の過程において，最適なノズル条件のみをデータセットに追加した場合，ニューラルネットワークの精度は最適解の周りは高まるがデータ点の無いところは不十分となる。この状態で再度最適化を行うと最初に得られた最適解周りに新しい最適解が得られるため，本来データ点の無いところに最適解があったとしても見つけることはできない。よって，データ密度の小さい場所のデータ点を求めてデータセットに追加することが有効となる。動径基底関数ネットワークは，既存のデータ点上に基底関数を置き，すべての点について和を取ったもので，これによりデータ密度を測ることができる。基底関数として，本稿ではガウス関数を適用する。動径基底関数ネットワークの最小点をデータセットに追加すれば，探索空間の全様をニューラルネットワークに学習させることができる。

以降では，鋳型内の溶鋼流動および伝熱凝固の数値解析モデルである数学モデルと浸漬ノズルの形状最適化の方法を述べ，浸漬ノズル形状最適化の実施例を報告する。

2. 数学モデル

2・1　支配方程式

鋳型内の溶鋼流動および伝熱凝固の数値解析には，Takataniらによる数学モデル^10,11)を用いる。支配方程式は，流体計算を非定常乱流モデルLES(Large Eddy Simulation)，伝熱凝固を非定常のエネルギー方程式としている。本モデルにおける仮定を下記に示す。

(1)流体は非圧縮であり，液相と固相の密度と比熱は同じく定数とする。

(2)固液共存相には流体抵抗としてDarcy則を適用する。

(3)固相は剛体として作用する(固相の速度は鋳造引き抜き速度と等しい)とし，固相の体積率が0.8より大きな領域は剛体とする。

(4)ガス気泡は非圧縮かつ分散相であり，形状は球体，直径は変わらないとする。

(5)気相の運動量方程式はBasset-Boussinesq-Oseen-Tchen方程式²²⁾とする。ただし，Basset項は無視できるとする。

(6)エネルギーバランスにおいて，気相のエンタルピーは無視できるとする。

(7)乱流モデルはLESを適用する。

(8)乱流Prandtle数は1とする。

(9)凝固温度はてこの原理で計算する²³⁾。

支配方程式は次の通りである。

  

$$f_g+f_l+f_s=1$$(1)

  

$$\frac{\partial f_g}{\partial t}+\nabla \cdot \left(f_g{u}_g\right)=0$$(2)

  

$$\frac{\partial f_l}{\partial t}+\nabla \cdot \left(f_l{u}_l\right)=-R_v$$(3)

  

$$\frac{\partial f_s}{\partial t}+\nabla \cdot \left(f_s{u}_s\right)=R_v$$(4)

  

$${\rho }_g\frac{D{u}_g}{Dt}={\rho }_l\frac{D{u}_l}{Dt}+\frac{1}{2}{\rho }_l\left[\frac{D{u}_l}{Dt}-\frac{D{u}_g}{Dt}\right]-\frac{3{\mu }_l}{4d_g^2}{\text{Re}}_g{C}_{Dg}\left(u_g-u_l\right)+\left({\rho }_l-{\rho }_g\right)g$$(5)

  

$$\begin{array}{l}\frac{\partial \left(f_l{\rho }_l{u}_l\right)}{\partial t}+\nabla \cdot \left(f_l{\rho }_l{u}_l{u}_l\right)=\\ -f_l\nabla p+\nabla \cdot \left[f_l\left({\mu }_l+{\mu }_t\right)\left(\nabla u_l+\nabla u_l^T\right)\right]+f_l\left[\beta \left(T-T_0\right){\rho }_{l}g+\gamma \left(u_s-u_l\right)+F_{lg}\right]\end{array}$$(6)

  

$$\frac{\partial T}{\partial t}+\nabla \cdot \left(f_l{u}_{l}T+f_g{u}_{g}T\right)=\nabla \cdot \left[\left(\alpha +{\alpha }_t\right)\nabla T\right]+\frac{1}{C_p}{R}_v\left(-\Delta H\right)$$(7)

  

$$T=\frac{\left(T_l-T\right)/m_l}{\left(T-T_s\right)/m_s+\left(T_l-T\right)/m_l}$$(8)

ここでtは時間，uは流速，fは体積率，pは圧力，Tは温度，ρは密度，μは粘性係数，C_pは比熱，αは熱伝導率，βは熱膨張率，gは重力加速度，R_vは溶鋼の凝固率を表す。下付き添え字のl,s,gはそれぞれ液相，固相，気相である。係数Re_gおよびC_Dgは式(9)と式(10)で与える。

  

$${\text{Re}}_g=d_g\left|u_l-u_g\right|\frac{{\rho }_l}{{\mu }_l}$$(9)

  

$$C_{Dg}=0.4+\frac{24}{{\text{Re}}_g}+\frac{6}{1+\sqrt{{\text{Re}}_s}}$$(10)

また，μ_t, α_tは乱流粘性係数および乱流熱伝導率を表わす。乱流粘性係数は式(11)で与える。

  

$${\mu }_t={\rho }_l{\left(C_s{\Delta }^{1/3}\right)}^2\sqrt{2\left|S\right|\left|S\right|}$$(11)

ここで係数C_s=0.1であり，Δはセルの体積，Sは式(12)で定義されるひずみ速度テンソルである。

  

$$S=\frac{1}{2}\left(\nabla u_l+\nabla u_l^T\right)$$(12)

式(6)中の右辺第3項は温度による浮力であり，T₀は液相の平均温度である。第4項は固液共存相の流体抵抗を表し，係数γを式(13)で与える。

  

$$\gamma =\frac{{\mu }_l}{{\rho }_l}\frac{150{\left(1-f_l\right)}^2}{D_2^2{f}_l^3}$$(13)

ここでD₂は2次デンドライトアーム間隔である。第5項は液相と気相の運動量交換であり，F_lgを式(14)で与える。

  

$$F_{lg}=f_g\frac{3{\mu }_l}{4d_g^2}{\text{Re}}_g{C}_{Dg}\left(u_g-u_l\right)$$(14)

2・2　計算条件および計算方法

計算領域はFig.1に示すように250 mm×1500 mm×3000 mm(厚さ×幅×長さ)の直方体領域とした。湯面から900 mmの範囲を鋳型として熱伝達係数を800 W/m²Kとし，それ以降は2次冷却と仮定して熱伝達係数は400 W/m²Kとした。浸漬ノズルの形状はFig.2に示す通りとし，浸漬深さd(mm)，吐出孔高さh(mm)，ノズル角度θ(degree)を連続なパラメータとして最適化により求める変数とした。計算領域の分割数は，長さ方向に516 mm以降の領域については36×80×23(厚さ×幅×長さ)で固定とした。516 mmまでの領域については，36×80(厚さ×幅)は固定であるが，長さ方向には浸漬深さと吐出孔高さの値に応じて，それぞれFig.1(b)における領域(D)および(H)の箇所の分割数を変更した。領域(D)については，格子間隔が25 mm未満となる最小の数を分割数とし，領域(H)については，25 mm以下となる最小の数を分割数とした。ただし，領域(H)の分割数が4未満の場合の分割数は4とした。

Fig. 1.

Computational Grid.

Fig. 2.

Geometry of the immersion nozzle. All dimension in mm, except angle degree.

スタッガード格子上に離散化しSOLA法²⁴⁾を基にした任意形状に対応した解法を用いて計算を行った。移流項には2次の中心差分スキーム，拡散項は2次の中心差分を用いて，移流項はEuler法により時間発展し，拡散項は陰的に計算を行った。

鋳造引き抜き速度は1.4 m/min，浸漬ノズルから鋳型へ流入するAr気泡の量は文献²⁵⁾を参考に，0.5 mm気泡径を2.668 Nℓ/min，1.0 mm気泡径を0.944 Nℓ/min，1.5 mm気泡径を0.493 Nℓ/minとした。なお，最適化の目的であるピンホールとブローホールはいずれも存在量が最も多い0.5 mm気泡径を対象とした。

2・3　気泡捕捉のメカニズム

気泡が凝固界面に捕捉されるメカニズムとして，速度勾配によるSaffman力²⁶⁾，界面活性元素の濃度勾配に起因する界面張力差^27,28)，電磁力²⁹⁾などが作用していると考えられているが，その詳細は不明である。例えば，硫黄濃度が高いほど，ピンホールが発生しやすいことが報告されており³⁰⁾，界面活性元素である硫黄の濃度勾配による界面張力差が大きく影響していると考えられる。式(15)は，界面張力差により溶鋼中の気泡に作用する力F_Iを示す²⁷⁾。

  

$$F_I=\frac{8}{3}\pi R^2{C}_0\frac{\partial \sigma }{\partial C_L}\frac{V_s}{D_L}\left(1-K_{\epsilon }\right)\exp \left(-\frac{V_s\left(x-\delta \right)}{D_L}\right)$$(15)

ここでC_Lは境界層中溶質濃度，C₀は凝固前の液相の溶質濃度，Rは気泡半径，σは界面張力，K_εは実効分配係数，V_sは凝固界面の成長速度，xは凝固界面からの距離，δは濃度境界層厚み，D_Lは溶質の拡散係数である。実効分配係数K_εは，式(16)として表される。

  

$$K_{\epsilon }=\frac{K_0}{K_0+\left(1-K_0\right)\exp \left(-V_s\delta /D_L\right)}$$(16)

ここでK₀は平衡分配係数である。

式(16)による界面張力差による気泡への作用力を連続鋳造機の数値解析において考慮するためには，硫黄濃度分布と凝固による偏析を考慮した計算を行えばよいが，0.01~1 mm程度の濃度境界層厚みを考慮することは現在の計算能力では困難である。一方，この界面張力差による作用力は，濃度境界層厚みを気泡径よりも十分に小さくすれば，その効果がなくなる。濃度境界層厚みは，Burtonらが式(17)の関係を報告している³¹⁾。

  

$$\delta =0.285D_L^{1/3}{\nu }^{1/6}{\omega }^{-1/2}$$(17)

ここでνは動粘性係数であり，ωは液相の回転による角速度である。すなわち，濃度境界層厚みは，凝固界面全面の流速のみの関数となり，流速が速いほど境界層厚みは薄くなる。

以上から，連続鋳造機の数値解析におけるピンホールの評価方法として，凝固シェルに捕捉される気泡個数を溶鋼流速に依存するとしてモデル化することが適切である。気泡の凝固シェルへ捕捉される速度η_εを式(18)に示すように仮定した。

  

$${\eta }_{\epsilon }=n_{\epsilon }{R}_s{P}_{\epsilon }\left(\left|u\right|\right)$$(18)

ここでn_εは気泡体積率，R_sは凝固速度であり，P_ε(|u|)は溶鋼流速に依存した気泡の捕捉確率である。凝固シェルに捕捉された気泡の体積率n_εは，式(19)に示す通り，凝固シェル速度u_sによって輸送する。

  

$$\frac{\partial n_{\epsilon }}{\partial t}+\nabla \cdot \left(u_s{n}_{\epsilon }\right)=-{\eta }_{\epsilon }$$(19)

気泡の捕捉流速の閾値として，0.05 m/s³²⁾や0.15 m/s³³⁾の値が報告されているが，凝固界面位置の定義や流速の評価位置が異なるため，正確な値は不明である。さらに，閾値を0.05 m/sとした場合は気泡が捕捉され，0.051 m/sの場合は全く捕捉されないというモデルは不自然と考えられる。よって，捕捉確率は速度の関数として，速度が増加するほど確率を低下させる式(20)に示す指数関数型のモデルとした。

  

$$P_{\epsilon }\left(\left|u\right|\right)=\exp \left(-C_o\left|u\right|\right)$$(20)

指数関数型モデルでは，0.15 m/sの場合の捕捉確率が，1×10^-8程度となるように，式(20)の係数をC₀=100と調整した。また，CFDにおいて，凝固が開始し固相率を有する要素の流速は，Darcy則により流速が大きく低下するため，凝固界面流速として固相率は0.1の位置の流速を用いた。

2・4　ピンホールとブローホールの算出方法

ピンホールとブローホールは，CFD計算開始から実時間換算で1000秒経過し，定常に達した後，100秒間の時間平均化を行った値を評価した。具体的に，ピンホールy_Pとブローホールy_Bは式(21)と式(22)で求めた。

  

$$y_P=\left({\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{i=1}^{M_P}\frac{C_{i,P}^n{V}_{i,P}}{100}\delta t}}\right)/\left(V_c{\displaystyle {\sum }_{i=1}^{M_P}{V}_{i,P}}\right)$$(21)

  

$$y_B=\left({\displaystyle {\sum }_{n=1}^N{\displaystyle {\sum }_{i=1}^{M_B}\frac{C_{i,B}^n{V}_{i,B}}{100}\delta t}}\right)/{\displaystyle {\sum }_{i=1}^{M_B}{V}_{i,B}}$$(22)

ここで，Nは時間の最大ステップ数，δtは時間ステップ幅，M_PとM_Bはそれぞれシェル側面の内，表層から2 mmまでのセルの総数と2 mm~15 mmのセルの総数，V_i,PとV_i,Bはそれぞれ当該領域におけるi番目のセルの体積，Cⁿ_i,PとCⁿ_i,Bはそれぞれ時間nステップ目の当該領域におけるi番目のセル内に捕捉されたAr気泡の体積率を表す。また，V_cは直径0.5 mmの球の体積，すなわちV_c=π(0.5×10^-3)³/6である。したがって，ピンホールは個数密度，ブローホールは体積率として評価している。

3. 浸漬ノズル形状最適化の方法

3・1　問題設定

本稿で考える浸漬ノズル形状最適化では，浸漬深さd，吐出孔高さh，ノズル角度θを設計変数とするため，探索空間を式(23)のように設定する。

  

$$S=\left\{{\left(d,h,\theta \right)}^T\in {ℝ}^3:d_{lb}\le d\le d_{ub},h_{lb}\le h\le h_{ub},{\theta }_{lb}\le \theta \le {\theta }_{ub}\right\}$$(23)

設計限界があるため，探索空間には下限値d_lb, h_lb, θ_lbと上限値d_ub, h_ub, θ_ubを設けた。今回用いた上下限値をTable 1に示す。

Table 1. Upper and lower bounds of the search space.

dlb(mm)	hlb(mm)	θlb(degree)
125	60	15
dub(mm)	hub(mm)	θub(degree)
225	120	40

浸漬ノズルの上記3変数が決まればCFDによりピンホールとブローホールが計算できるので，この工程を関数Fとして式(24)のように表す。

  

$${\left(y_P,y_B\right)}^T=F\left(d,h,\theta \right)$$(24)

このとき，ピンホールとブローホールを低減するための浸漬ノズル形状最適化問題は式(25)のように定式化できる。

  

$$\text{Find}{\left(d^{*},h^{*},{\theta }^{*}\right)}^T=\underset{{\left(d,h,\theta \right)}^T\in S}{\text{argmin}}L\left(F\left(d,h,\theta \right)\right)$$(25)

ここでL=L(y_P, y_B)は目的関数を表している。ピンホールとブローホールを同時に最小化することが本研究の目的であるため，問題は多目的最適化問題に分類される。これを単純化するために，本稿では，式(26)のように2つの目的を平均で定義される目的関数に落とし込むことで単一最適化問題とした。

  

$$L\left(y_P,y_B\right)=\frac{1}{2}\left|y_P\right|+\frac{1}{2}\left|y_B\right|$$(26)

一般に最適化は目的関数を評価するために複数回モデルFを呼び出すため，計算コストの観点から問題(26)を直接解くことは現実的ではない。そこで，本稿では，Fを使わず，Fを近似する計算コストの低いモデル $\tilde{F}$ としてニューラルネットワークを使用した。このとき，最適化問題(26)は，式(27)で表される近似問題となる。

  

$$\text{Find}{\left({\tilde{d}}^{*},{\tilde{h}}^{*},{\tilde{\theta }}^{*}\right)}^T=\underset{{\left(d,h,\theta \right)}^T\in S}{\text{argmin}}L\left(\tilde{F}\left(d,h,\theta \right)\right)$$(27)

さらに，最適解の精度を高めるために，データ点を段階的に追加して，ニューラルネットワークの精度を高めながら，反復的に最適解の更新を行う逐次近似最適化^18,19)の枠組みを採った。結局，式(28)で表される問題を解くこととした。

  

$$\text{Find}{\left({\tilde{d}}^k,{\tilde{h}}^k,{\tilde{\theta }}^k\right)}^T=\underset{{\left(d,h,\theta \right)}^T\in S}{\text{argmin}}L\left({\tilde{F}}^k\left(d,h,\theta \right)\right),k=1,2,\cdots$$(28)

3・2　ニューラルネットワークによるピンホールとブローホールの予測

本ニューラルネットワークモデルは，浸漬ノズルの設計変数h, d, θを入力として，ピンホールとブローホールy_P, y_Bを出力とする。教師あり学習を利用してモデルを作成した。データセットΩ={(x_j, y_j)=(h_j, d_j, θ_j, y_P,j, y_B,j)|j=1,…, m}が与えられたとき，まずデータセットΩを学習用のデータセットΩ_trainと試験用のデータセットΩ_testに7：3の割合で分割した。その後，Ω_trainの入出力変数を式(29)で標準化したデータセットΩ'_trainを作成した。

  

$$x^{\prime }=\frac{x-\mu }{\sigma }$$(29)

ここでx'はxを標準化した値，μとσは各変数の平均と分散である。また，これらの平均と分散を用いてΩ_testの入出力変数を式(29)で変換したデータセットΩ'_testを作成した。

ニューラルネットワークFは，式(30)のように複数の関数F_iが入れ子構造または層構造になった関数である。

  

$$\tilde{F}\left(x\right)={\widehat{F}}_N\left({\widehat{F}}_{N-1}\left(\cdots \left({\widehat{F}}_1\left(x\right)\right)\right)\right)$$(30)

関数F_iは，より具体的に式(31)で表される。

  

$${\widehat{f}}_i\left({\widehat{y}}_{i-1}^{\text{'}}\right)=a\left(W_i{\widehat{y}}_{i-1}^{\text{'}}\right)$$(31)

ここでW_iは重み行列， ${\widehat{y}}_{i-1}^{\text{'}}$ は式(32)で定義されるi-1層の出力ベクトル ${\widehat{y}}_{i-1}^{}$ の一行目に1を挿入したベクトルである。

  

$${\widehat{y}}_{i-1}=\{\begin{array}{cc}{\widehat{F}}_{i-1}\left({\widehat{F}}_{i-2}\left(\cdots \left({\widehat{F}}_1\left(x\right)\right)\right)\right)& ifi\ge 2\\ x,& ifi=1\end{array}$$(32)

関数aは活性化関数と呼ばれ，双曲線正接関数やランプ関数(ReLU)を始めとして様々な形状のものが提案されている³⁴⁾。

重みは事前には決定できないため，事前に決めた誤差関数を最小化することでニューラルネットワークがΩ'_trainの入出力に一致するように重みを更新すなわち学習する。学習では，確率的勾配降下法³⁵⁾を用いた。誤差関数としては式(33)で定義される平均二乗誤差に正則化項を加えたものを採用した。

  

$$E\left(y,\tilde{y};W\right)=\frac{1}{m}{\displaystyle \sum }_{j=1}^m\left\{{\left(y_{P,j}-{\tilde{y}}_{P,j}\right)}^2+{\left(y_{B,j}-{\tilde{y}}_{B,j}\right)}^2\right\}+\lambda {\left({\displaystyle \sum _{i,j}{\left|w_{i,j}\right|}^p}\right)}^{\frac{1}{p}}$$(33)

ここでyと $\tilde{y}$ はそれぞれ与えられた出力データと予測値であり，w_i,jは重み行列の成分，λは正則化の程度を表すパラメータを表している。ただし正則化項の和はニューラルネットワークすべての重みについて取られるものであり，またp≥1である。

正則化項は過学習を防止することを目的としたものであるが，加えて，ドロップアウトも使用した。ドロップアウトは各層の出力ベクトルの成分を一定の割合で間引き，次層の入力変数の次元を落とす方法である³⁶⁾。これにより，確率的勾配降下法で局所解へ落ち込むことを防ぐ効果が得られる。さらに，学習を早期に完了させるため，重みの初期値を文献³⁹⁾に倣って設定した。以上のニューラルネットワークの作成においてはPythonのライブラリであるTensorFlow³⁷⁾を用いた。

ニューラルネットワークの構造を決定するパラメータとして，層数N，各層の入力変数の次元n_unit，ドロップアウトの割合r，活性化関数の種類，正則化項におけるpおよびλ，確率的勾配降下法における学習率ηがあり，これらは学習より前に決定しておく必要がある。本稿では，ベイズ最適化を用いて各パラメータの最適化を行った。Table 2に調整対象のパラメータとその探索範囲を示す。実装では，PythonのライブラリであるOptuna³⁸⁾を用いた。各パラメータの最適候補値の決定と学習およびテストデータΩ'_testを用いた予測を繰り返し，最もテストデータの予測精度が良いものを最終的なニューラルネットワークモデルとした。なお，ベイズ最適化における試行回数は1000回とし，テストデータの予測精度は決定係数R²で測った。

Table 2. Parameters of NN structure and search range.

N	nunit	r
2 to 7	1 to 100	0.0 to 0.9
p	λ	η
1 or 2	100 to 10−3	100 to 10−4

3・3　粒子群最適化

粒子群最適化^20,40)は，ランダムに位置する複数の粒子(探索点)が，決められた規則に則って探索空間を移動しながら最適解を探していく方法である。ここではニューラルネットワークによるピンホールおよびブローホールの予測を組み込んだ方法について述べる。

一つの粒子の位置は浸漬ノズルの設計変数で構成される探索空間内のベクトルによってx_i=(d_i,h_i,θ_i)^Tと表される。ここでiは粒子の番号である。また，各粒子は位置ベクトルのほかに移動ベクトルも持っており，探索空間の次元のベクトルv_i=(v_i1, v_i2, v_i3)^Tで表される。さらに，各粒子はそれまでの探索で発見したそれぞれの最良解，すなわち目的関数が最も小さくなった時の位置lbest_i=(lbest_i1, lbest_i2, lbest_i3)^Tと，その目的関数の値L( $\tilde{F}$ (lbest_i))を記憶している。群れとしては，すべての粒子がこれまでの探索で発見した最良解gbest=(gbest₁, gbest₂, gbest₃)^Tを記憶している。

Fig.3に計算流れ図を示す。まず指定した粒子数だけ探索空間内にランダムに粒子の位置x_i⁽⁰⁾と速度v_i⁽⁰⁾を生成する。次にニューラルネットワークを用いて各粒子に対する目的関数の値L( $\tilde{F}$ (x_i⁽⁰⁾))を計算し，最良解lbest_iを記憶する。その後は，設定した最大の反復回数に到達したか収束するまで以降の処理を繰り返す。反復回数をtとするとき，現在の位置x_i^(t)から，それぞれの最良解へ向かうベクトルlbest_i^(t)-x_i^(t)，群れ全体の最良解gbest^(t)-x_i^(t)，および前回の移動ベクトルv_i^(t)の重み付き線形結合として新たな移動ベクトルv_i^(t+1)を生成する。t+1回目の移動におけるi番目の粒子の移動ベクトルは式(34)で与えられる。

  

$$v_i^{\left(t\right)}=w{v}_i^{\left(t\right)}+c_1{r}_1\left(lbes{t}_i^{\left(t\right)}-x_i^{\left(t\right)}\right)+c_2{r}_2\left(gbes{t}^{\left(t\right)}-x_i^{\left(t\right)}\right)$$(34)

Fig. 3.

Neural Network-based Particle Swarm Optimization process.

ここでr₁とr₂は0から1の間に分布する一様乱数，w, c₁, c₂は事前に設定する重みパラメータである。そして，式(35)によって各粒子の位置を次の位置x_i^(t+1)に移動する。

  

$$x_i^{\left(t+1\right)}=x_i^{\left(t\right)}+v_i^{\left(t\right)}$$(35)

なお，x_i^(t+1)のj番目の成分x_ij^(t+1)が下限値を下回った場合にはx_ij^(t+1)を下限値とし，同様に，x_ij^(t+1)が上限値を上回った場合にはx_ij^(t+1)を上限値とした。

本稿では，終了条件として最大反復回数を与えることとし，Table 3に挙げる条件で最適化を実施した。

Table 3. Particle Swarm Optimization parameters.

Number of particles	w	c1	c2	Max iteration
100	0.5	0.14	0.14	500

3・4　動径基底関数ネットワークによる密度関数を用いたサンプル点の追加

動径基底関数ネットワークは入力層，隠れ層1層，出力層の3層からなるニューラルネットワークである。設計変数からなる既存のデータ点を{x_j}_j=1,…,mとする。ただし，重複はないものとする。この時，動径基底関数ネットワークの出力は式(36)で与えられる。

  

$$z\left(x;\alpha \right)={\displaystyle \sum }_{i=1}^m{\alpha }_{i}K\left(x,x_i\right)$$(36)

ここでK(x, x_i)はi番目のユニットの基底関数，α=(α₁,…,α_m)は基底関数の重みである。本稿では基底関数として式(37)で定義されるガウス関数を用いる。

  

$$K\left(x,x_i\right)=\exp \left(-\frac{{\left(x-x_i\right)}^T\left(x-x_i\right)}{r_i^2}\right)$$(37)

ここでr_jはi番目のデータ点の半径である。半径の定義式は様々なものが提案されている^41,42,21)が，本稿では，文献⁴²⁾に倣い式(38)を用いた。

  

$$r_i=\frac{d_{i,\text{max}}}{\sqrt{n}\sqrt[n]{m-1}}$$(38)

ここでnは探索空間の次元，d_i,maxはi番目のデータ点と他のデータ点との距離の内，最大のものを表している。本稿ではd_i,maxとして式(39)に示すようにユークリッド距離を採用した。

  

$d_{i,\text{max}}=\underset{j=1,\cdots ,m}{\max }\Vert x_j-x_i\Vert =\underset{j=1,\cdots ,m}{\max }\sqrt{{\left|d_j-d_i\right|}^2+{\left|h_j-h_i\right|}^2+{\left|{\theta }_j-{\theta }_i\right|}^2}$(39)

密度関数は，動径基底関数ネットワークを用いてサンプル点の疎な領域に極小値を生成させる関数である。密度関数を用いたサンプル点の決定方法は，以下の手順で行われる。

(1)すべての成分を1とした列ベクトルを用意する。

  

$$y={\left(1,\cdots ,1\right)}_{m\times 1}^T$$(40)

(2)式(41)で係数α*=(α*₁,…,α*_m)を求める。

  

$$a^\ast=\left(K^{T}K+\mathcal{E}\right)^{-1}K^{T}y\tag{41}$$

ここでK,Eはそれぞれ式(42)と式(43)で定義される行列である。式(43)中のεは重みパラメータであり，本稿ではε=1.0×10^-2とした。

  

$$K=\left[\begin{array}{cccc}K\left(x_1,x_1\right)& K\left(x_1,x_2\right)& \cdots & K\left(x_1,x_m\right)\\ K\left(x_2,x_1\right)& K\left(x_2,x_2\right)& \cdots & K\left(x_2,x_m\right)\\ ⋮& ⋮& \ddots & ⋮\\ K\left(x_m,x_1\right)& K\left(x_m,x_2\right)& \cdots & K\left(x_m,x_m\right)\end{array}\right]$$(42)

  

$$\mathcal{E} =\begin{bmatrix} \varepsilon & 0 & \ldots & 0 \\ 0 & \varepsilon & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \ldots & \varepsilon \end{bmatrix} \tag{43}$$

(3)次の最小化問題(44)を解き，得られた最小解を追加サンプル点とする。なお，密度関数とは本問題の目的関数を指す。

  

$$\text{Find}{x}^{\text{*}}=\text{arg}\underset{x\in S}{\min }{\displaystyle \sum }_{i=1}^m{\alpha }_i^{*}K\left(x,x_i\right)$$(44)

3・5　実装方法

Fig.4に全体の計算流れ図を示す。最適化開始時点ではニューラルネットワークの学習に使用するデータがないため，CFDを用いてピンホールとブローホールを計算することで，その時のノズル条件と合わせてデータセットを作成する。ニューラルネットワークの外挿性能は高くないため，ここで用意するデータセットは探索範囲の境界に相当する条件を中心に数点選択した。具体的には，浸漬深さ(mm)が60, 100, 120，吐出孔高さ(mm)が125, 175, 155, 195, 225，ノズル角度(degree)が15, 20, 25, 30, 35, 40である時の組み合わせた条件を計算した。

Fig. 4.

Sequential Approximate Optimization process.

その後，当該データセットを用いてニューラルネットワークの最適な構造を探すとともに学習を行い，得られたモデルを用いて粒子群最適化を実行する。粒子群最適化により，その時に最適な浸漬ノズルの設計変数値が得られる。併せて動径基底関数ネットワークを用いた密度関数により，データがない領域からノズル条件を求めて追加データ点とする。粒子群最適化と密度関数により得られたノズル条件に対して，正確なピンホールとブローホールを得るためにCFDを実行する。ここまでで得られたノズル条件とピンホールおよびブローホールをまとめてデータセットに追加する。ここで，解候補をできるだけ多く確保してデータセットを拡充することでニューラルネットワークの予測精度を短工期にあげるため，粒子群最適化と密度関数によるサンプル点の追加は3回ずつ繰り返すこととした。以上の処理を反復して行い，所定の回数を終えた時点で全処理を終了とする。

4. 最適化結果と考察

Fig.5に，黒丸として最適化開始前に用意した浸漬ノズル条件のCFD計算値(base condition)，黒三角として逐次近似最適化の反復を10回行って得られた浸漬ノズル条件でのニューラルネットワークによる予測値(prediction by neural network)，白丸として同条件並びに動径基底関数ネットワークによる密度関数で求めた浸漬ノズル条件のCFD計算値(optimal conditions & search points)を示す。本稿では，理想である最小のピンホールとブローホールをいずれも0としているため，図の左下ほど良い条件である。ニューラルネットワークによる最適化結果(黒三角)は，手動での探索結果(黒丸)よりも良化する条件を予測している。この条件でCFD計算を実施した結果(白丸)とニューラルネットワークの予測値(黒三角)は一致しない場合もあるが，概ね同じ傾向を示すことが確認された。よって，ニューラルネットワークでの予測と確認のためのCFD計算を自動的に繰り返すことにより，効率的かつ短時間で手動よりも良い条件を見出すことが可能である。

Fig. 5.

Scatter plot for Pinhole and Blowhole.

Fig.5中の条件Aは手動による最良点(黒丸)，条件BはCFD確認計算(白丸)の最良点を示している。以降，これらのことを単に基本条件，最適条件と呼ぶ。Table 4に基本条件と最適条件のピンホールとブローホールを示す。いずれの数値も最適条件の場合が小さくなっており，ピンホールで約56%，ブローホールで約10%低減できた。

Table 4. Pinholes and Blowholes by condition.

	Base condition	Optimal condition
Pinhole	2322	1002
Blowhole	13.8	12.3

Fig.6にそれぞれのノズル形状を示す。基本条件の浸漬ノズルについては，吐出孔高さとノズル角度が今回考えた探索領域における下限値近くの値となっていることから，水平に近く速い吐出流を得ることがピンホールとブローホールを低減する傾向にあることが分かる。Fig.7(a)は基本条件でのCFD計算で得られた流速を示しており，鋳型短辺に衝突した吐出流が上下に分岐している。このような流れはダブルロール流動パターンとして知られている⁴³⁾。ダブルロール流動パターンで生じる上昇流は湯面に到達後，湯面を攪拌することでAr気泡が凝固界面に捕捉されることを抑制するため，ピンホールを小さくすることができる。また，吐出流とダブルロール流動パターンによる上昇/下降流により凝固界面に捕捉されたAr気泡の洗浄効果も生み，ブローホールを小さくすることもできる。

Fig. 6.

Comparison of immersion nozzle shapes.

Fig. 7.

Comparison of velocity in the mold.

最適条件は，これらの効果をより効果的に得られる結果となっている。浸漬ノズルの形状については，吐出孔高さは基本条件同等であるが，角度をより浅くしていることから，より水平方向への吐出流を生む結果となっている。Fig.7(b)に見て取れるように，この効果によって強い上昇流が得られ，湯面の流速が強くなりピンホール改善につかながったと考えられる。浸漬深さが深くなった理由としては，鋳型全面に流速を生じさせ，凝固シェルへ捕捉された気泡を流すことでブローホールを低減するためと考えられる。Fig.8に流速の大きさを示す。基本条件に対して最適条件では，全体的に流速が大きくなっていることが確認できる。また，湯面に到達した流れが浸漬ノズルのある鋳型中央部で下方に流れると，同時にAr気泡を運ぶためブローホールが悪化する。よって，浸漬ノズルを深くすることで湯面の中央部に向かう流れを強くしすぎないようにしていると考えられる。

Fig. 8.

Comparison of velocity magnitude in the mold.

5. 結言

本稿では連続鋳造機における鋳型内溶鋼流動および伝熱凝固のCFD計算とニューラルネットワークを用いた最適化法を組み合わせた浸漬ノズル形状最適化方法並びにその適用例について報告した。本法では，CFD計算に鋳型内溶鋼流動と伝熱凝固の数学モデル，最適化の枠組みに逐次近似最適化，CFD計算の代理モデルにニューラルネットワーク，最適化アルゴリズムに粒子群最適化を使用して組み合わせた。本法を，ピンホールとブローホールを低減する浸漬ノズル形状探索に適用した結果，事前にCFD計算により求めたピンホールとブローホールに対して，ピンホールで約56%，ブローホールで約10%低減することができた。また，浸漬深さを深く，吐出孔高さを低く，ノズル角度を浅くしたノズル形状が良い傾向が見られた。これにより，水平に近く速度の速い吐出流と強い上昇流が得られ，湯面がかき混ぜられることでピンホールを低減できたと考えられる。また，鋳型内全面にかけて凝固界面に付着したAr気泡を洗い流す効果も得られることでブローホールを低減できたと考えられる。

今後の課題としては，EMSやEMBrを考慮した場合のノズル形状最適化が挙げられる。また，本稿では，ピンホールとブローホールの評価値を平均化することで単一目的最適化として解いたが，多目的最適化とした場合の検討も必要であると考えられる。

謝辞

冨田克行氏には，ニューラルネットワークの作成および学習の効率化について参考文献と重要な知見をいただき，本論文の校正においてご助言いただいた。また，畠中健氏と齋藤佑丞氏にはCFD計算の実施や浸漬ノズル形状，鋳型内溶鋼流動とピンホールおよびブローホールの関係について知見をいただいた。査読者には，原稿について有益で詳細なコメントをいただいた。ここに諸氏への感謝の意を表する。

文献

• 1)  S.Yokotani, S.Haseo, Y.Asako, S.Takagi, K,Ayata, J.Szekely and S.Hara: Tetsu-to-Hagané, 82(1996), 581 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.82.7_581
• 2)  J.Yoshida, M.Iguchi and S.-I.Yokota: Tetsu-to-Hagané, 88(2002), 264 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.88.5_264
• 3)  K.Takatani: Tetsu-to-Hagané, 90(2004), 751 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.90.10_751
• 4)  Y.Tsukaguchi, H.Hayashi, H.Kurimoto, S.Yokoya, K.Marukawa and T.Tanaka: Tetsu-to-Hagané, 95(2009), 33 (in Japanese). https://doi.org/10.2355/tetsutohagane.95.33
• 5)   B.W. Li,  X.Y. Tian,  E.G. Wang and  J.C. He: Acta Metall. Sin.-Engl. Lett., 20(2007), 15. https://doi.org/10.1016/S1006-7191(07)60003-9
• 6)   H.-S. Park,  H. Nam and  J.K. Yoon: ISIJ Int., 41(2001), 974. https://doi.org/10.2355/isijinternational.41.974
• 7)   H. Liu,  C. Yang,  H. Zhang,  Q. Zhai and  Y. Gan: ISIJ Int., 51(2011), 392. https://doi.org/10.2355/isijinternational.51.392
• 8)   R.D. Morales,  Y. Tang,  G. Nitzl,  C. Eglsäeer and  G. Hckl: ISIJ Int., 52(2012), 1607. https://doi.org/10.2355/isijinternational.52.1607
• 9)   M. Xuan and  M. Chen: Metals, 11(2021), 1223. https://doi.org/10.3390/met11081223
• 10)   K. Takatani,  Y. Tanizawa,  H. Mizukami and  K. Nishimura: ISIJ Int., 41(2001), 1252.https://doi.org/10.2355/isijinternational.41.1252
• 11)   K. Takatani: ISIJ Int., 43(2003), 915. https://doi.org/10.2355/isijinternational.43.915
• 12)   L. Leifsson and  S. Koziel: J. Comput. Sci., 10(2015), 45. https://doi.org/10.1016/j.jocs.2015.01.003
• 13)   C. Yang and  F. Huang: J. Hydrodyn., 28(2016), 947. https://doi.org/10.1016/S1001-6058(16)60696-0
• 14)   E. Bernardini,  S.M.J. Spence,  D. Wei and  A. Kareem: J. Wind Eng. Ind. Aerodyn., 144(2015), 154. https://doi.org/10.1016/j.jweia.2015.03.011
• 15)   X. Song,  L. Lv,  W. Sun and  J. Zhang: Struct. Multidiscip. Optim., 60(2019), 965. https://doi.org/10.1007/s00158-019-02248-0
• 16)   H. Wu,  X. Liu,  W. An,  S. Chen and  H. Lyu: Comput. Fluids, 198(2020), 104393. https://doi.org/10.1016/j.compfluid.2019.104393
• 17)   J. Li,  X. Du and  J.R.R.A. Martins: Prog. Aerosp. Sci., 134(2022), 100849. https://doi.org/10.1016/j.paerosci.2022.100849
• 18)  H.Nakayama, Y.Yun and M.Yoon: Sequential Approximate Multiobjective Optimization Using Computational Intelligence, Springer, Berlin, Heidelberg, (2009), 113. https://doi.org/10.1007/978-3-540-88910-6
• 19)  S.Kitayama: Kougakukei No Tameno Saiteki Sekkeihou (Design Optimization for Engineering), Kyoritsu Shuppan, Tokyo, (2021), (in Japanese).
• 20)  J.Kennedy and R.Eberhart: Proc. ICNN’95 – Int. Conf. on Neural Networks, Vol.4, IEEE, Piscataway, NJ, (1995), 1942. https://doi.org/10.1109/ICNN.1995.488968
• 21)   S. Kitayama,  M. Arakawa and  K. Yamazaki: Optim. Eng., 12(2011), 535. https://doi.org/10.1007/s11081-010-9118-y
• 22)  J.O.Hinze: Turbulence, 2nd ed., McGraw Hill, New York, (1975), 460.
• 23)  Z.-G.Wang, M.Tanaka and T.Inoue: Trans. Jpn. Soc. Mech. Eng. A, 53(1987), 1735. (in Japanese). https://doi.org/10.1299/kikaia.53.1735
• 24)  C.W.Hirt, B.D.Nichls and N.C.Romero: SOLA-A Numerical Solution Algorithm for Transient Fluid Flows, Los Alamos Science Laboratory, Los Alamos, NM, (1975), LA-5852.
• 25)   T. Toh,  H. Hasegawa and  H. Harada: ISIJ Int., 41(2001), 1245. https://doi.org/10.2355/isijinternational.41.1245
• 26)   P.G. Saffman: J. Fluid Mech., 22(1965), 385. https://doi.org/10.1017/S0022112065000824
• 27)  K.Mukai and W.Lin: Tetsu-to-Hagané, 80(1994), 533 (in Japanese). https://doi.org/10.2355/tetsutohagane1955.80.7_533
• 28)   K. Mukai,  L. Zhong and  M. Zeze: ISIJ Int., 46(2006), 1810. https://doi.org/10.2355/isijinternational.46.1810
• 29)   S. Taniguchi and  J.K. Brimacombe: Tetsu-to-Hagané, 80(1994), 312. https://doi.org/10.2355/tetsutohagane1955.80.4_312
• 30)   T. Miyake,  M. Morishita,  H. Nakata and  M. Kokita: ISIJ Int., 46(2006), 1817. https://doi.org/10.2355/isijinternational.46.1817
• 31)   J.A. Burton,  R.C. Prim and  W.P. Slichter: J. Chem. Phys., 21(1953), 1987. https://doi.org/10.1063/1.1698728
• 32)  Y.Miki, H.Ohno, Y.Kishimoto and S.Tanaka: Tetsu-to-Hagané, 97(2011), 423 (in Japanese). https://doi.org/10.2355/tetsutohagane.97.423
• 33)   A.M. Mollinger and  F.T.M. Nieuwstadt: J. Fluid Mech., 316(1996), 285. https://doi.org/10.1017/S0022112096000547
• 34)   A. Apicella,  F. Donnarumma,  F. Isgrò and  R. Prevete: Neural Netw., 138(2021), 14. https://doi.org/10.1016/j.neunet.2021.01.026
• 35)  I.Goodfellow, Y.Bengio and A.Courville: Deep Learning, MIT Press, Cambridge, MA, (2016).
• 36)   N. Srivastava,  G. Hinton,  A. Krizhevsky,  I. Sutskever and  R. Salakhutdinov: J. Mach. Learn. Res., 15(2014), 1929.
• 37)  TensorFlow, Google LLC: https://www.tensorflow.org/, (accessed 2021-10-20).
• 38)  T.Akiba, S.Sano, T.Yanase, T.Ohta and M.Koyama: Proc. 25th ACM SIGKDD Int. Conf. on Knowledge Discovery and Data Mining, Association for Computing Machinery, New York, (2019), 2623. https://doi.org/10.1145/3292500.3330701
• 39)  X.Glorot and Y.Bengio: Proc. 13th Int. Conf. on Artificial Intelligence and Statistics (OMLR9), Microtome Publishing, New York, (2010), 249.
• 40)  E.Aiyoshi and K.Yasuda: MetaHeuristics To Ouyou (MetaHeuristics and Applications), IEEJ, Tokyo, (2007), 69 (in Japanese).
• 41)  S.Haykin: Neural Networks (A Comprehensive Foundation), Macmillan College Publishing, New York, USA, (1994).
• 42)   H. Nakayama,  M. Arakawa and  R. Sasaki: Optim. Eng., 3(2002), 201. https://doi.org/10.1023/A:1020971504868
• 43)   B.G. Thomas and  L. Zhang: ISIJ Int., 41(2001), 1181. https://doi.org/10.2355/isijinternational.41.1181