2019 年度
卒業論文

様々な最適化手法を用いた超伝導変圧器の機器定数推定

一木悠人(学籍番号：16232010)2020 年 2 月 17 日

九州工業大学

情報工学部

電子情報工学科

指導教員：小田部

荘司

教授

目次
第1章

序論

1

1.1

はじめに . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1

1.2

研究背景 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

1.3

本稿における解析系統 . . . . . . . . . . . . . . . . . . . . . . .

31.4

変圧器の簡易磁気飽和特性と回路方程式 . . . . . . . . . . . . . .

4

1.5

常微分方程式の計算手法 . . . . . . . . . . . . . . . . . . . . . .

6

1.6

最適化アルゴリズム . . . . . . . . . . . . . . . . . . . . . . . . .

8

1.7

先行研究 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

1.8

本研究の目的 . . . . . . . . . . . . . . . . . . . . . . . . . . . .

14

推定手法

17

2.1

計算手法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

17

2.2

差分進化法を用い機器定数の初期値を一つ固定した場合の推定 .

17

2.3

複数のアルゴリズムでの推定 . . . . . . . . . . . . . . . . . . . .

18

2.4

計算環境 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

結果及び考察20

3.1

機器定数を一つ固定し DE を用いた場合 . . . . . . . . . . . . . .

20

3.2

複数のアルゴリズムでの推定結果 . . . . . . . . . . . . . . . . .

21

3.3

焼きなましスケジュールを変化させたときの SA . . . . . . . . .

26

3.4

近傍範囲の初期値を変化させたときの SA . . . . . . . . . . . . .

26

結論

29

第2章

第3章

第4章

i

謝辞

30

研究業績

31

参考文献

32

ii

図目次
1.1

超伝導転移温度最高記録の推移. . . . . . . . . . . . . . . . . .

21.2

系統の解析モデル [2] . . . . . . . . . . . . . . . . . . . . . . . .

41.3

単相超伝導変圧器に対する等価回路 [2] . . . . . . . . . . . . . .

4

1.4

ヒステリシス曲線 [3] . . . . . . . . . . . . . . . . . . . . . . . .

5

1.5

i − ϕ 簡易磁気特性 [2] . . . . . . . . . . . . . . . . . . . . . . . .

6

1.6

ϕ − L 特性 [2] . . . . . . . . . . . . . . . . . . . . . . . . . . . .

7

1.7

DE で推定された波形 . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1

測定波形 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

3.1

Lc の初期値を固定した時の評価関数の変化 . . . . . . . . . . . .

21

3.2

Lm の初期値を固定した時の評価関数の変化 . . . . . . . . . . . .

22

3.3

Lair の初期値を固定した時の評価関数の変化 . . . . . . . . . . .

22

3.4

Bn の初期値を固定した時の評価関数の変化 . . . . . . . . . . . .

23

3.5

θ0 の初期値を固定した時の評価関数の変化 . . . . . . . . . . . .

23

3.6

ϕres の初期値を固定した時の評価関数の変化 . . . . . . . . . . .

24

3.7

各アルゴリズムで計算した励磁突入電流 . . . . . . . . . . . . . .

25

3.8

各アルゴリズムの評価関数の世代変化 . . . . . . . . . . . . . . .

26

3.9

焼きなましスケジュールを変化させたときの評価関数 . . . . . .

27

3.10

探索範囲の初期値を変化させたときの評価関数 . . . . . . . . . .

28

iii

表目次
1.1

DE での推定結果 . . . . . . . . . . . . . . . . . . . . . . . . . .

15

2.1

解候補の探索範囲 . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2.2

計算環境 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19

3.1

機器定数を 1 つ固定したときの推定結果. . . . . . . . . . . . .

20

3.2

各アルゴリズムでの評価関数と推定時間 . . . . . . . . . . . . . .

24

3.3

各アルゴリズムで推定した機器定数 . . . . . . . . . . . . . . . .

25

iv

第1章

序論

1.1 はじめに1911 年オランダの物理学者オネンスが，水銀の電気抵抗が 4.2 K 以下で消失することを発見した．これが超伝導現象の初の発見である．1933 年にはドイツのマイスナーとオクセンフィルドが超伝導体に強い反磁性があることを明らかにした．これがいわゆるマイスナー効果である．1957 年，超伝導の仕組みを解明したのはアメリカの物理学者バーディーン，クーパー，シュリーファーであった．その理論を 3 人の頭文字をとって BCS 理論と呼ぶ．BCS 理論は超伝導を量
子力学の現象として扱い量子現象のマイクロレベルでの出現という超伝導の本質を説明した．しかし，BCS 理論の枠内では臨界温度は 40 K を超えないという予
測も導き出された．臨界温度とは超伝導現象を示さなくなる温度のことである．
この後，約 30 年間臨界温度 40 K を超える超伝導体は発見されなかった．しかし 1986 年，ドイツの物理学者ベトノルツ，スイスの物理学者ミュラーにより同による酸化物が臨界温度 30 K を示し，世界中で銅酸化物の超伝導体の研究が進められた．その 3 ヶ月後に臨界温度 92 K の Y 系銅酸化物の超伝導体が発見された．さらに，YBa2 Cu3 Oy によって液体窒素 (77 K) での超伝導実験が可能になった．1988 年には臨界温度 100 K を超える Bi – Sr – Ca – Cu – O 系の超伝導体が発見されている．
図 1.1 からもわかるように超電導体の臨界温度が液体窒素の沸点より高い温度で実現されたことなどから現在，主に線材としての応用が期待されて

1

図 1.1 超伝導転移温度最高記録の推移

いる．超伝導体は電気抵抗が 0 Ω という特徴があり，エネルギーの損失な
く大電流を流すことが可能なためである．また，医療現場で使用されている

MRI(Magnetic Resonance Imaging，磁気共鳴画像) やリニアモーターカー，電力・エネルギー分野での運用が期待されている SMES(Superconducting Magnetic

Energy Storage，超伝導エネルギー貯蔵装置)，SUQID(Superconducting QUantum
Interference Device，超伝導量子干渉計) といった高感度センサーなど様々な分野での応用が期待されている．[1]1.2 研究背景
現在の電力系統において使用している変圧器を電力系統へ投入する際，励磁突
入電流が発生し，他の電力系統の電力品質に影響を与えている．これらの影響には，渦電流などの発生によって伴う保護リレーの不要動作，瞬時電圧低下による制御装置の誤作動，照明類のちらつき，消灯などがある．これらの影響について把握，対策するためには任意の投入条件の励磁突入電流発生時の挙動を解析できることが必要となる．解析に必要な機器定数のうち，定格磁束密度や空芯インダ
クタンスは励磁突入電流の大きさを左右するものでありながら，変圧器の設計に密にかかわる値であることから，系統運用や設備運用において解析を行う場合入2

手が困難となる場合が生じている．また，現在超伝導変圧器は特に，短絡電流等による巻線部分でのクエンチ時の送電停止が問題となっており実用化はまだ実現されていないが，将来的に実用化が期待されている．超伝導現象では先述したよ
うに電気抵抗が 0 Ω になることから，変圧器の巻線に超伝導体を用いることに
よって巻線部分での損失が非常に小さくなり，エネルギー効率を向上させることができる．またエネルギー損失が小さくなることで電気装荷が大きくなり，さらに鉄心の断面積を小さくすることによって小型化・軽量化が可能になるなどの利点がある．将来，超伝導変圧器が電力系統に導入される場合にも，常伝導変圧器と基本原理は同じであるため，解析に必要な機器定数が入手困難になる状況が生じると考えられる．
また一方で，電力系統の電力品質に与える影響の原因となっている変圧器の励磁突入電流のみを測定できる場合がある．変圧器高圧側には計器用変流器が設置されており，その出力回路にクランプを使用することにより非接触で測定できるからである．このような場合，逆に励磁突入電流から変圧器の機器定数推定と投入条件の値が推定できる．これまでの研究では電力系統で超伝導変圧器を運転する場合にも系統運用における解析検討が円滑に行えるよう最適ルンゲ・クッ
タ法と遺伝的アルゴリズム (Genetic Algorithm,GA) や差分進化法 (Deﬀerential

Evolution,DE) を組み合わせた推定方法を構成し，単相の Bi-2223 線材による超伝導変圧器の機器定数推定が行われている．
そこで本稿では，一つの機器定数の初期値を固定して DE を用いた推定や，他のアルゴリズムを用いた推定方法で機器定数推定を行った．その結果を比較・考察する．

1.3 本稿における解析系統
本解析方法では，図 1.2 に示す単線結線図の系統を解析モデルとする．上位系
統は 1 電源，1 回線の送電線で構成されている．解析対象の変圧器を高圧側から系統へ投入するとし，低圧側に接続した開閉器は常時開放されている．実際の設備では，変圧器の低圧側母線にコンデンサ，所内変圧器，避雷針などの機器が常時接続されている場合があるが，容易に解析ができるようにこのモデルに固定し3図 1.2 系統の解析モデル [2]

図 1.3 単相超伝導変圧器に対する等価回路 [2]た．解析で求める電流は図 1.2 の電流計での瞬時値とした．図 1.2 の系統における単相の等価回路を図 1.3 に示す．対象の変圧器は励磁突入電流に対してクエン
チしないことから，巻き線抵抗は 0 Ω とした．また磁気特性は本来，図 1.4 のようにヒステリシス曲線を描くが，鉄心の抵抗成分は無視できるほど十分に大きい
ことから，損失成分は無視し，磁気飽和特性を反映したインダクタンスの部分だ
けを扱うようにした．図 1.3 において計算に必要な系統，機器定数は以下の通りである．

• RB , LB : 上位系統のインピーダンス• Lc : 変圧器の漏れインダクタンス• Lc1 = Lc /2 : 高圧側の漏れインダクタンス• LL : 磁気飽和特性を反映した励磁インダクタンス

1.4

変圧器の簡易磁気飽和特性と回路方程式の等価回路における LL は定格運転時の日飽和時には励磁インダクタンス Lm ，飽和時には励磁インダクタンスとして考えられる巻き線だけの自己インダクタン
ス，すなわち空芯インダクタンス Lair をとる．その様子を図 1.5 に示す．その値
4図 1.4 ヒステリシス曲線 [3]は各部分で dΦ/dt となり，ほぼ一定と考えてよい．そこで最適ルンゲ・クッタ法
では，各相の磁束 ϕ の変化に対して階段状に励磁インダクタンス LL が変化する
磁気飽和特性を用いた．この特性を図 1.6 に示す．図 1.6 からわかるように Lm ，

Lair の値は非飽和時，飽和時のそれぞれの時間帯において一定の値を取る．以上のことを踏まえ，回路方程式を式 (1.1) から式 (1.3) で表す．

dim
+ RA im = V(t)dt(π )V(t) = Vm sin ωt +
θ0180
dim {V(t) − RA im }
=
dt
Lk
Lk =

(1.1)
(1.2)
(1.3)

ここで，図 1.3 の場合，回路の全抵抗 RA ，全インダクタンス Lk はそれぞれ
c1 +LLRA = RB ，Lk = LB +L
において計算に必要な系統，機器定数は以下の通りであ10000

る．LL は磁気飽和特性を表す励磁インダクタンス，Vm は電源電圧の波高値，θ0は電圧投入位相である．また，ω については ω = 2π fe で与えられ， fe は電源電圧の周波数である．磁束 ϕ は式 (1.1) を変形して式 (1.4) から得られ，LL については，式 (1.5)，(1.6) から求めることができる．
5

図 1.5 i − ϕ 簡易磁気特性 [2]∫ϕ=∫V(t)dt − RA

im dt

(1.4)





Lm , |ϕ| ≤ ϕmax
LL = 
Lair , ϕmax > |ϕ|( )BSϕmax = ϕn
Bn

(1.5)

(1.6)

ここで ϕmax は飽和磁束値，Bn は変圧器の定格磁束密度，BS は飽和磁束密度，ϕnは定格磁束密度である．以上から，励磁突入電流 im を求めることができる．

1.5 常微分方程式の計算手法
本稿で励磁突入電流の計算に用いた常微分方程式の数値計算法は，1 段階法で陽的な解法であり，ルンゲ・クッタ法に比べ誤差が小さくなるよう，A. Ralston
によって改良された最適ルンゲ・クッタ法である．この計算方法は時間の刻み幅∆t に対して係数が変化することから，Φ − L 特性の階段状の変化に対応できる．
i を時間 t の関数 i(t) とし，i の時間微分 di/dt を関数 f (t, i) とおくと最適ルンゲ・クッタ法は式 (1.7) で定式化される．

6

図 1.6 ϕ − L 特性 [2]k00 = f (tn , in )k1 = in + α10 k00
k11 = f (tn + τ1 , k1 )k2 = in + α20 k00 + α21 k11
k22 = f (tn + τ2 , k2 )k3 = in + α30 k00 + α31 k11 + α32 k22
k33 = f (tn + τ3 , k3 )in+1 = i(tn + ∆t) = in + (c0 k00 + c1 k11 + c2 k22 + c3 k33 )∆t
ここで，上記の式で用いた定数は次のようにする．

7

(1.7)

τ1 = 0.4∆tτ2 = 0.45573725∆t

τ3 = ∆tα10 = 0.4∆t

α20 = 0.29697761∆t

α21 = 0.15875964∆t

α30 = 0.21810040∆t

α31 = −3.05096561∆t

c0 = 0.17476028

c1 = −0.55148066

c2 = 1.20553560

c3 = 0.17118478

α32 = 3.8328647∆t

上記の式 (1.7) から励磁突入電流を求めることができる．

1.6 最適化アルゴリズム
1.6.1 差分進化法
差分進化法 (DE) とは，進化的アルゴリズムの一種であり，確率的な直接探索によって解集団を用いた多点探索を行うアルゴリズムである．差分進化法は非線形問題，微分不可能な問題などの現実的な実行時間では厳密な解を求めることが困難な問題に対して近似解を求めることが可能であり，様々な最適化問題に適応されている．以下にアルゴリズムの流れを示す．

1. ベクトル集団の作成指定された各次元の最小値と最大値の範囲で一様乱数によって N 個の初期個体 XRi の生成を行うことで，初期集団 {XRi ,i = 1, 2, · · · , N} を構成する．

2. 変異
各個体 XRi に対し，4 つの個体 XR0 ,XR1 ,XR2 ,XR3 を重複しないようにランダ
ムに選択する．新しいベクトル U(mutant vector) を基本ベクトル XR1 及び差分ベクトル XR2 − XR3 から式 (1.8) によって生成する．

U = XRi + F(XR2 − XR3 )(1.8)

スケーリングファクター F はユーザが値を設定する変数で [0,1] の区間で定義されるアルゴリズムの制御パラメータである．

3. 交叉
8

親ベクトル (parent vector) を XR0 とし，ベクトル U と XR0 の各次元の要素
を交叉率 CR の確率で U の要素から継承し，1 − CR の確率で XR0 から継承することで子ベクトル (child vector) を生成する．交叉率 CR はユーザが値を設定する変数で [0,1] の区間で定義されるアルゴリズムの制御パラメータ
である．

4. 選択
子ベクトルが親ベクトル XR0 よりも良い推定であった場合，親ベクトル XR0を子ベクトルに置換する．

5. 終了判定
終了条件を満たした場合，アルゴリズムは終了される．満たしていない場
合 2 から 5 を繰り返す．

1.6.2 粒子群最適化
粒子群最適化とは，差分進化法と同じように進化的アルゴリズムの一種であり，解集団を用いた多点探索を行うアルゴリズムである．以下にアルゴリズムの流れを示す．

1. ベクトル集団の作成指定された各次元の最小値と最大値の範囲で一様乱数によって N 個の初期個体 Xi の生成を行うことで，初期集団 {Xi ,i = 1, 2, · · · , N} を構成する．

2. 位置の更新
各粒子の位置 Xi を速度 Vi を用いて式 (1.9) で更新する．

Xi = Xi + Vi

(1.9)

3. 速度の更新
各粒子の Vi を式 (1.10) で更新する

Vi = w Vi + c1 r1 (Xgb − Xi ) + c2 r2 (Xlb − Xi )(1.10)

ここで，w は慣性定数，c1 ，c2 は良い解に向かう粒子の割合，r1 , r2 は [0,1]9の区間の一様乱数を取る．Xgb は解候補群の中で最も評価関数の良い個体の位置，Xlb は各個体で今までに最も評価関数の良かった位置を示す．

4. 終了判定
終了条件を満たした場合，アルゴリズムは終了される．満たしていない場
合 2 から 5 を繰り返す．

1.6.3 重力探索アルゴリズム
重力探索アルゴリズムとは，差分進化法と同じように進化的アルゴリズムの一種であり，解集団を用いた多点探索を行うアルゴリズムである．万有引力の法則と質量の相互作用からその名がつけられている．以下にアルゴリズムの流れを示す．

1. ベクトル集団の作成指定された各次元の最小値と最大値の範囲で一様乱数によって N 個の初期個体 Xi の生成を行うことで，初期集団 {Xi ,i = 1, 2, · · · , N} を構成する．

2. 個体の評価
それぞれの個体を評価関数を用いて評価する．

3. 個体間の距離の更新
個体間の距離 Ri j を式 1.11 を用いて更新する．

Ri j =√(

Xi1 − X 1j

)2

+ ··· +(

Xin − X nj

)2

(1.11)

4. 万有引力定数の更新
万有引力定数を式 (1.12) で更新する．
αt

G(t) = G(t0 ) e− T

(1.12)

ここで，α は G(t) の減衰度合い，T は総ループ回数を示す．

5. 個体に加わる力を更新
それぞれの個体に加わる力を式 (1.13) で更新する．

10

Fi j (t) = G(t)Fi (t) =∑)Mi (t) × M j (t) (
X
(t)
−
X
(t)
j
i
R2i j + ε

rand j Fi j (t)(1.13)

j,i

6. 質量の更新
それぞれの質量 Mi を以下の式 (1.14) で更新．

d − dworst
dbest − dworstmiM i = ∑ni=1 mimi =(1.14)

7. 個体の加速度を更新
それぞれの加速度を式 (1.15) で更新する．

ai (t) =Fi (t)Mi (t)(1.15)

8. 個体の速度を更新
それぞれの速度を式 (1.16) で更新する．

Vi (t + 1) = randi Vi (t) + ai (t)(1.16)

9. 個体を位置を更新
それぞれの位置を式 (1.17) で更新する．

Xi (t + 1) = Xi (t) + Vi (t + 1)

(1.17)

10. 終了判定
終了条件を満たした場合，アルゴリズムは終了される．満たしていない場
合 2 から 5 を繰り返す．

1.6.4 ネルダーミード法
ネルダーミード法とは，1965 年に J. A. Nelder と R. Mead によって発表された
最適化アルゴリズムである [5]．n 次元の最適化を行う場合，n + 1 の頂点を持つ多面体を変形させながら最適解を探索する．以下にアルゴリズムの流れを示す．
11

1. ベクトル集団の作成指定された各次元の最小値と最大値の範囲で一様乱数によって n + 1 個の頂点の初期状態を設定する．

2. ベクトル集団のソート
ベクトル集団を評価関数の大小でソートする．最良点を x1 ，第二最悪点を

xn ，最悪点を xn+1 とする．
3. 最悪点以外の点の重心を求める
最悪点を除く n 点の重心を式 (1.18) で求める．

x0 =

n
∑

xi

(1.18)

i=1

4. 反射点を求める
反射点 xreflect を式 (1.19) で求める．

xreflect = 2x0 − xn+1

(1.19)

反射点が最良点よりも悪く，第二最悪点よりも良い場合，最悪点を反射点に置換する．

5. 拡大点を求める反射点が最良点よりも良い場合，拡大点を式 (1.20) で求
める．

xexpand = x0 − 2(xreflect − x0 )(1.20)= 3x0 − 2xn+1(1.21)

拡大点が反射点よりも良い場合，最悪点を拡大点に置換する．拡大点が反射点よりも悪い場合，最悪点を反射点に置換する．

6. 収縮点を求める
反射点が第二最悪点よりも悪い場合，収縮点を式 (1.22) で求める．

xcontract = 0.5x0 − 0.5xn+1(1.22)

収縮点が最悪点よりも良い場合，最悪点を収縮点に置換する．収縮点が最悪点よりも悪い場合，最良点を除く点全てを式 (1.23) で縮小する．
12

xi = 0.5x1 + 0.5xi(1.23)

7. 終了判定
終了条件を満たした場合，アルゴリズムは終了される．満たしていない場
合 2 から 5 を繰り返す．終了条件を満たした場合，アルゴリズムは終了される．満たしていない場合 2 から 6 を繰り返す．

1.6.5 焼きなまし法
焼きなまし法とは大局的最適化問題への乱択アルゴリズムの一種である．以下にアルゴリズムの流れを示す．

1. 初期位置の設定指定された各次元の最小値と最大値の範囲で一様乱数によって初期状態を設定する．

2. 近傍探索
現在の位置の近傍を探索する．

3. 選択
現在の位置と近傍の評価を行い近傍が現状態よりも良い評価をされた場
合，現在の位置に置換する．現在の位置よりも悪くなった場合，ある確率で現在の位置に置換する．評価関数が悪くなったときの遷移確率は評価関数の差分と，現在の温度によって決まる．温度の時間変化を焼きなましス
ケジュールという．

4. 終了判定
終了条件を満たした場合，アルゴリズムは終了される．満たしていない場
合 2 と 3 を繰り返す．
今回の焼きなまし法での近傍探索範囲と焼きなましスケジュールは以下のようにした．近傍探索は現在の位置の各要素を式 (1.24) で探索する．

13

xi = xi ± n[%]( i)
n = n0 exp −I(1.24)

ここで n0 は近傍探索範囲の初期値である．
遷移確率と焼きなましスケジュールは，式 (1.25)，(1.26) で表す．

(e − e′
r = exp
T)(1.25)i

T = αI(1.26)

ここで e は現在の位置での評価関数，e′ は近傍の位置の評価関数である．また，

0 < α < 1 であり，i はこれまでのループ回数，I は終了までのループ回数である．

1.7

先行研究

先行研究では，最適ルンゲ・クッタ法と GA や DE を組み合わせて推定が行われた．DE と最適ルンゲ・クッタ法を組み合わせて行った機器定数推定の結果を以下に示す．F = 0.8, CR = 0.8 で DE を用いて 10 回推定を行い，その結果を表1.1 に示す．またその最良の推定結果から計算された励磁突入電流と測定された励磁突入電流をを求めたところ図 1.7 のようになった．図 1.7 を見ると，波形のピークの部分で特に大きな誤差が生じていることがわかる．表 1.1 を見ると，θ0以外の機器定数で分散が大きいことがわかる．そのため，十分に改善の余地があると考えられる．

1.8 本研究の目的
超伝導現象では電気抵抗が 0 Ω になるため，変圧器の巻線に超伝導体を用いる
ことで，巻線部分での損失が非常に小さくなり，エネルギー効率が向上する．
変圧器を電力系統に投入するときに励磁突入電流が生じる．その影響を把握・対策するためには波形レベルでの解析が必要であり，その解析のためには変圧器の機器定数が必要である．変圧器の設計のためには機器定数のうち，定格磁束密度と空芯インダクタンスが特に重要となるが系統運用や設備運用での解析での入手が困難となる場合がある．一方で，変圧器高圧側にある計器用変流器のクラン
14表 1.1 DE での推定結果

Lc [mH]Lm [mH]

Lair [mH]

Bn [T]

θ0 [° ]

ϕres [%]

d [A]

5.63 × 101

8.04 × 103

1.70 × 102

1.58

1.96 × 102−6.43 × 102

9.04

6.79 × 101

1.23 × 104

1.69 × 102

1.58

1.95 × 102−6.53 × 102

9.27

6.07 × 101

8.04 × 103

1.67 × 102

1.61

1.96 × 102−6.18 × 102

9.04

6.41 × 101

8.04 × 103

1.65 × 102

1.63

1.96 × 102−5.99 × 102

9.04

5.27 × 101

8.04 × 103

1.71 × 102

1.69

1.96 × 102−5.55 × 102

9.04

6.50 × 101

7.63 × 103

1.62 × 102

1.62

1.95 × 102−5.87 × 102

9.00

6.38 × 101

8.00 × 103

1.66 × 102

1.63

1.96 × 102−5.98 × 102

9.04

5.74 × 101

5.38 × 103

1.58 × 102

1.61

1.95 × 102−5.61 × 102

8.96

7.27 × 101

8.01 × 103

1.62 × 102

1.67

1.96 × 102−5.75 × 102

9.04

6.68 × 101

8.04 × 103

1.64 × 102

1.58

1.96 × 102−6.40 × 102

9.04

図 1.7 DE で推定された波形

15クを用いて励磁突入電流が得られ，それから機器定数が推定できる．
先行研究では GA や DE と最適ルンゲ・クッタ法を用いて超伝導変圧器の機器
定数推定を行っていた．しかし，まだ測定した波形と推定した波形に主にピーク部分で誤差がある．そこで，本研究では DE を用い機器定数の初期値を一つ固定した推定や，PSO・GSA・SA・NM を用いた推定を行うことで超伝導変圧器の機器定数のより精度の良い推定を行うことを目的とする．

16

第2章

推定手法

2.1 計算手法
本研究では，第 1 章で述べた各最適化アルゴリズムと最適ルンゲ・クッタ法を用いて超伝導変圧器の機器定数推定を行う．測定波形を図 2.1 に示す．測定波形
は 1.0 × 10−4 s ごとに測定しているので最適ルンゲ・クッタ法の ∆t = 1.0 × 10−4 sとし，各時間における測定値と計算値の値の距離を求め，その総和をサンプル数で割ることで誤差の平均を求める．今回はこれを評価関数とする．評価関数を式

2.1 に示す．
1∑|iMEA − iSIM |
d=
N

(2.1)

ここで，iMEA は測定値，iSIM は計算値，N はサンプル数を表す．今回推定する機器定数は，漏れインダクタンス Lc [mH]，励磁インダクタンス Lm [mH]，空芯イ
ンダクタンス Lair [mH]，定格磁束密度 Bn [T]，初期位相 θ0 [° ]，残留磁束 ϕres [%]である．これらの値を式に代入し計算し，式 (2.1) によって評価し，解を探索していく．また表 2.1 に解候補の探索範囲を示す．

2.2 差分進化法を用い機器定数の初期値を一つ固定した場合の推定
今回の研究では以下のように推定を行った．
機器定数の初期値を一つを一定間隔ごとに固定して，差分進化法 (DE) を用い
17

図 2.1 測定波形

表 2.1 解候補の探索範囲

機器定数

探索範囲

漏れインダクタンス

51.3 mH < Lc < 73.8 mH

励磁インダクタンス

Lc < Lm < 16000 mH

空芯インダクタンス

0 mH < Lair < 10 Lc

定格磁束密度1.55 T < Bn < 1.7 T

初期位相0° ≤ θ0 < 360°

残留磁束−85% ≤ ϕr ≤ 85%

機器定数の推定を行う．DE のパラメータは F = 0.8, CR = 0.8，世代数 10000 で行う．解候補の初期値によって，推定結果が変わるため，10 回推定を行う．

2.3 複数のアルゴリズムでの推定
今回の研究では以下の条件で推定を行った．

• 粒子群最適化 (PSO) では w = 0.8，c1 = 0.95，c2 = 0.95，個体数 20，世代
数 10000 で推定を行った．
18

表 2.2 計算環境

OS

Windows 10

プロセッサIntel Core i7 960

コア数

4

スレッド数

8

動作周波数

3.2 GHz• 重力探索アルゴリズム (GSA) では個体数 20，世代数 10000，α = 2 で推定を行った．

• 焼きなまし法 (SA) では，α = 0.2， step = 0.01，世代数 10000 で推定を行った．

• ネルダーミード法 (NM) は世代数 10000，推定する変数が 6 つであるため，7 点で推定を行った．
解候補の初期値によって，推定結果が変わるため，粒子群最適化，重力探索アル
ゴリズム，焼きなまし法，ネルダーミード法による推定をそれぞれ 10 回ずつ行
う．重力探索アルゴリズム，粒子群最適化では 10 回の推定の解候補群の最良の値を用い，焼きなまし法，ネルダーミード法では 10 回の推定の評価関数の最小値を用いて比較する．

2.4 計算環境
表 2.2 に計算環境を示す．

19

第3章

結果及び考察3.1 機器定数を一つ固定し DE を用いた場合
まず，超伝導変圧器の六つの機器定数の初期値を一つ固定して DE を用い推定した結果を表 3.1 に示す．表 3.1 からわかるように，すべての機器定数を推定した時は d = 8.96 A であったのに対して，機器定数の初期値を一つ固定した場合
は，最良でも Lm = 5.40 × 103 % の時に，d = 8.97 A であった．DE を用い，すべての機器定数を推定した時よりも良い推定を行うことができた機器定数はな
かった．これから，すべての機器定数の推定で十分な推定ができていることがわ
かった．

Lc を固定して推定を行った結果を図 3.1 に示す．
Lm を固定して推定を行った結果を図 3.2 に示す．
Lair を固定して推定を行った結果を図 3.3 に示す．
表 3.1 機器定数を 1 つ固定したときの推定結果

初期値を固定した機器定数

評価関数 d [A]

評価関数が最小の時の値

Lc

9.09

6.75 × 101 mH

Lm

8.97

5.40 × 103 mH

Lair

9.07

1.55 × 102 mH

Bn

9.09

1.59 Tθ0

9.75

1.95 × 102 °

ϕres

8.98

−5.10 × 101 %20

図 3.1

Lc の初期値を固定した時の評価関数の変化

Bn を固定して推定を行った結果を図 3.4 に示す．
θ0 を固定して推定を行った結果を図 3.5 に示す．
ϕres を固定して推定を行った結果を図 3.6 に示す．
Lm ，Lair ，θ0 ，ϕres では，定数の初期値によって評価関数の値が大きく変わっていることがわかる．これから，機器定数の推定では Lm ，Lair ，θ0 ，ϕres を精度よく推定しなければならないことが分かった．特に θ0 では値によって，d = 93.3 [A]となり，測定波形からの誤差が非常に大きくなってしまう．

3.2

複数のアルゴリズムでの推定結果

各アルゴリズムでの評価関数と推定時間を表 3.2 に，推定した機器定数を表3.3 に示す．表 3.2 からわかるように，DE よりも評価関数が良くなっているアル
ゴリズムはなかった．PSO・GSA では，DE よりも良い推定を行うことはできな
かったが，これは PSO は解集団の最良値を使用し，GSA は評価関数の良い値の解の加重を大きくして推定を行っているため，2 つのランダムな個体の差分を用い変異する DE よりも大域的な探索が行えなかったと考えられる．そのため，最21図 3.2

Lm の初期値を固定した時の評価関数の変化

図 3.3

Lair の初期値を固定した時の評価関数の変化

22図 3.4

Bn の初期値を固定した時の評価関数の変化

図 3.5 θ0 の初期値を固定した時の評価関数の変化

23図 3.6 ϕres の初期値を固定した時の評価関数の変化
表 3.2 各アルゴリズムでの評価関数と推定時間

アルゴリズム

評価関数 d [A]

推定時間 [ms]

DE

8.96

794

PSO

14.8

14900

GSA

11.3

16500

SA

8.99

389

NM

9.06

542

良値が局所解にあった場合，解集団全体が局所解に収束してしまうので，DE よりも評価関数の値が悪くなったと考えられる．一方，SA と NM はほぼ同じ評価関数の値で SA は DE の約 21 の時間で，NM は DE の約 23 の時間で推定できた．
これは，DE は解候補群を用いて解を探索し，毎世代，2 つの個体の差分を求め
基本ベクトルに足し合わせることをおこなうが，SA は単一の解候補で解の探索を行い，NM は一世代で最良点と第二最悪点の間ならば反射点を求めるだけで良いという違いがあるためである．今回のような 6 次元の最適化問題では DE を使
わずとも，SA や NM で十分な最適化ができ，より短い時間で推定を行うことが24

表 3.3 各アルゴリズムで推定した機器定数

アルゴリズムLc [mH]Lm [mH]

Lair [mH]

Bn [T]

θ0 [° ]

ϕres [%]

GSA

6.41 × 101

7.25 × 103

1.25 × 102

1.65

1.95 × 102−4.78 × 101

PSO

5.87 × 101

3.81 × 103

8.31 × 101

1.63

1.91 × 102−3.47 × 101

SA

5.61 × 101

6.09 × 103

1.59 × 102

1.64

1.95 × 102−5.35 × 101

NM

8.55 × 101

6.97 × 103

1.43 × 102

1.48

1.95 × 102−6.64 × 101

図 3.7 各アルゴリズムで計算した励磁突入電流できる．
表 3.8 に各アルゴリズムの評価関数の世代変化を示す．SA が他のアルゴリズ
ムと比較して解の改善速度が遅いが，これは DE，GSA，PSO，NM が複数の解候補を変化させながら探索しているのに対して，SA が一つの解候補を変化させながら解を探索しているため，解候補の初期位置によって解の改善が遅くなってしまうためであると考えられる．
表 3.7 に各アルゴリズムで推定した機器定数から計算した励磁突入電流を示す．DE と比較して PSO 以外ではピーク部分での改善が行われなかった．PSO
ではピーク部分では改善したが，立ち上がりの部分で DE よりも誤差が大きくなってしまった．
25

図 3.8 各アルゴリズムの評価関数の世代変化

3.3 焼きなましスケジュールを変化させたときの SA
SA において焼きなましスケジュールを決定する α を 0 < α < 1 と 0.1 ずつ変化させたときの評価関数の変化を図 3.9 に示す．図 3.9 より，α を変化させても評価関数の値はあまり変化しないことが分かった．これは近傍探索範囲をループを重ねるごとに小さくしているため，焼きなましスケジュールが推定に大きな影響を与えず，解候補の初期値が大きく影響したと考えられる．

3.4 近傍範囲の初期値を変化させたときの SA
SA において近傍範囲の初期値を変化させて 10 回推定を行ったときの評価関数の最小値の変化を図 3.10 に示す．現在の値からそれぞれの機器定数を 0.1 –9% 変化させた時は，0.2 – 1% は良い推定を行えたが，0.1% がとても悪く，近傍探索範囲の初期値が探索範囲の 2% を超えると，大きくなっていくにしたがって評価関数が悪くなっていった．これは，近傍探索範囲の初期値を小さくすると，

26

図 3.9 焼きなましスケジュールを変化させたときの評価関数

初期値の周囲だけで探索を行うため，良い推定を行うことができず，近傍探索範
囲の初期値が大きくすると，探索終盤での解の詳細な探索が行えないため，評価関数が悪化したと考えられる．

27

図 3.10 探索範囲の初期値を変化させたときの評価関数

28

第4章

結論
今回の研究では，DE を用い機器定数の一つを固定しての超伝導変圧器の機器定数推定や，PSO，GSA，SA，NM を用いた推定を行った．計算の結果，一つの機器定数の初期値を固定しての推定や複数のアルゴリズムを用いて推定を行ったが先行研究で行われた DE での推定よりもより良い推定を行うことができなかった．これから，現在の解析モデルではこれ以上の改善が難しいと考えられるため，評価関数やより良い推定を行うためには解析モデルの改善が必要であると考えられる．PSO，GSA では，DE よりも良い推定を行うことはできなかったがこ
れは，PSO は解集団の最良値を使用し，GSA は評価関数の良い値の解の加重を大きくして推定を行っている．そのため，最良値が局所解にあった場合，解集団全体が局所解に収束してしまうので，DE よりも評価関数の値が悪くなったと考えられる．SA，NM では，DE と同程度の精度で，より短い時間で推定を行うことできたため，今回のような 6 次元の最適化問題では SA や NM は DE よりも良
いと考えられる．

29

謝辞
本研究を進めるにあたり，ご指導いただきました私の指導教員である九州工業
大学情報工学研究院物理情報工学研究系 小田部荘司教授に深くお礼申し上げます．次に，松下照男名誉教授からは，研究の基礎となる，超伝導の基礎について教授いただきました．心から感謝申し上げます．また，岡山県立大学情報工学部の徳永義孝准教授にも深くお礼申し上げます．私の研究の助けになる実験データや資料をいただいたことを心から感謝申し上げます．最後に，小田部研究室メン
バーには公私ともに助言をいただくなど大変世話になりました．改めてここに感謝の意を表します．

本研究は JSPS 科研費 19H00771 の助成を受けたものです．

30

研究業績
国内学会

1. 一木悠人, 米中友浩, 村岡樹, 小田部荘司, 徳永義孝 「2019 年度応用物理学会九州支部学術講演会 24Cp-5 差分進化法を用いた超伝導変圧器の機器定数推定」熊本大学工学部 2019 年 11 月 23 日から 24 日

31

参考文献
[1] 超伝導技術の現状と将来展望 Business ＆ Economic Review 1996 年 09 月号[2] 徳永 義孝「励磁突入電流波形に遺伝的アルゴリズムを適用した超電導変圧器の機器定数推定」

低温工学

53 巻 1 号 2018 年 17 から 22 ページ[3] 村岡 樹「差分進化法を用いた超伝導変圧器の機器定数推定」 九州工業大学情報工学部電子情報工学科論文[4] Esmat Rashedi, Hossein Nezamabadi-pour, Saeid Saryazdi. ”GSA: A Gravitational Search Algorithm”. Information Sciences 179 (2009) 2232 ‒ 2248
[5] Nelder, John A.; R. Mead (1965).

A simplex method for function minimiza-tion . Computer Journal 7: 308 ‒ 313.
[6] E S Otabe, T Taguchi, Y Tsuruda and R Funaki. ”Determination of Pinning
Parameters in Flux Creep-Flow Model for E − J characteristics of High Temperature”. Journal of Physics: Conference Series 871. 2017, 012046

32