︿ ろ

、 ， . . . ‑ ‑ ー

なお，上式においては次の偏微分係数を用いている。

ァi﹁

lI ll 11 1l﹂

︑_I

BE Ef

k  /y

︿ 仇 ハ

u ' b

z n  

K 

AX  

/tilt¥ σ﹄口♂

一 次

﹁l

11 11 11

1﹂

一 一

一

lk 

ι n  

A  εR ^l1Xl1ヲ

B_k/k

三 [ 三 ^{( g 丸山川} ) T

^{ER川 m}

^，

φk/k

三 [ 会 ( 丸 山 UJ k / ) ] T

^{εR川 l}

^，

⁽⁴⁰⁸⁾ 

C

斗 L h ^恥

叫 九 川 ぺ 三 剥 h

^{川山恥仇(い仏料￡k川ん/九/ル叫}

k λ1)]  ]  ] )

εRPX7Y

， 

εR pxl 

次に，時点 t

= 

t k において， (405)式より時間刻み幅似の将来を求める。

︑︑B︐ノOノハUA斗/

' 't︑︑ ' k 

v 

+ 

﹁I

ll

Il

li

‑‑

﹂‑

¥I ll

八ι ノi‑‑

zK

I J J /  

4

司 、 r

このとき，状態選移行列は Pade(l))近似により以下に求められる。なお，

Pade(l， 1)近似では U の2次 (=1+ 1)の項まで正確な級数展開と一致する。

F 句 χ l I1tA~ιfルlι" ) ) Y  J 

^j 

Fklk 三

叫 吋 l d t f ^{仇 叫}

^L

ヰ ι ^叫 j L

^九/ルIk

1 =  

=1 十叫A

叫

^k 

⁺  j ^{主 が ( い}

^A

⁴ 

^k 

^{) y 2}   _+→~( ^叫

^//川k

^{r  +} 

三

I + A f A T I L + 1 1 A f A ι / f i+11AfA f I J 3 +

バ/バ 2¥  '¥1"  /  4 ¥ 円'"/ 

= ( I ‑ h J ‑ ( I + h J  

i4 

(410) 

ー- ~

ちなみに， Pade(l))近似は，連続系で安定な極 (左半平面内)を離散時間 系で必ず安定な極(単位円内)に射影する性質を有し，時間刻み毎の誤差伝播 の 観 点 で 数 値 計 算 安 定 性 向 上 に 寄 与 し て い る 。 こ の 性 質 は ，

λ [  f ( A ) 卜 f ( λ [ A

J)なる関係を利用し，

A

の j 番 目 の 固 有 値 引

A ]

に着 目すれば，以下のとおり確認できる。

イ(r一千A:/kf(r+与A~/k ⁾  1 

̲ 1 1 + 4 4 4 / k ] 1

^2 ‑(

^{1 +}   _~t ^{A ) [ A Z / k  ] ) ( 1 +} 与 ^{[ A : / k ] )}  1

‑ _~t ^{A j [ A Z / k ]}  ( ‑ 1 ミ

^λ

^{j [ A Z / k}  ] ( ) l ‑ 4 4 1 4 / k l )  

(~t)21λ J [A Z /k ] 1

1  +  ~tRe{ 十 九 l} ⁺  \~.

^J 

P ; " k / k 

JI 

I 

( 411) 

こ こ に ， 肩 文 字 * は 複 素 共 役 を 示 す。

さらに，駆動行列は， (410)式 の 適 用 ， 及 び ， 操 作 量 の O 次 ホ ー ル ド

u ( t )  

^= U 

k    t ( k  ~ t

くら

+ 1 )

仮定により，以下のとおり導かれる。

G  k l k

^{三 (}

F k l k 一山 / k ) ‑ i ^凡 z A f ( I ‑h J ‑ l B j

^( 412) 

75 

ー 、 r

このとき， (401)式 の ブ ラ ウ ン 運 動 は 次 の 白 色 雑 音 h に離散化される[7]。

¥Il

li

‑‑

︑ ノ‑

I ll ‑

J︑︑₂

11 /

4 J L 4 ' L  

/I1︑︑/

1 1i t¥

β ζ  

/I

ll l‑

‑¥  

AU K 

#付LKB ︑llt﹀l

J' 1

︑︑t 1 5 1'

'J

ず

‑

iMH 4F

︐  

/I1¥ 

LK 

#付

IK

A r ‑

‑ E 1 J︑BI'

t︑

ny

 

x 

ρ L E

 

一﹁ Il

l‑

LH

1 k l  

' b 一 一

n  

v  (413) 

こ こ に ， 駆 動 行 列 に 基 づ き 円 は 次 の 統 計 量 を 持 つ。

k 

T M  

G 

¥Illi‑‑ノ

o G  

n y  

O 

Q 

fill‑

‑

¥  

‑b MH

' K  

 

G 一 一

zK  

ιMH 

Aハ

﹀︑

︑B1l︑f 一 一

' ︐ ノ

k 

v /14l︑

ri  

QU V 

(414) 

以上の諸 式は拡張カルマンフィルタのフレームを満たし，未知パラメータ θ の常時変動の条件で5章に述べる同定・状態推定が行える。

4.  4 数 値 計 算 の 安 定 化

前節に述べたとおリ，本論では(403)，(409)式の離散時間システム拡張カル マ ン フ ィ ル タ の 手 調7]を用いて，状態 Xk の推定，パラメータ 仇 の 同 定

を行う。このとき，数値計算誤差に係わる不具合を防止するため，片山の成書 [36]を参考に U‑D分解形式による推定分散・共分散行列

P

の正定対称性の 維持をはかり さらに，観測更新 ・時間更新((408)式の遷移行列算出を含む) 過程を通じ，逆行列演算を用いない工夫の下に推定，同定を行っている。

4.4. 1観 測 更 新 過 程

本節では観測更新前 (添え字 : k/k‑l )の諸量から，時点 k の観測 値を考慮、した観測更新後(添え字 : k / k )の諸量の安定な算出法を整理す

る

76 

可r~

まず，観測方程式(403)に左より，観測雑音のコレスキー分解 Rjl/2を乗じ ると以下の関係が得られる。

v 出 ^{; l}

^( 415) 

ここに，観測量の予測値は次のとおり求められる。

ん

^1k‑1

恥/ 山

^{d/k‑l ( 416)} 

従って，対角化された分散・共分散行列に基づく観測雑音混入に際し，次の 通り逆行列を全く用いることなくカルマンゲイン

Kk

が得られ，これを用い て，求める状態、量，パラメータの観測更新を行うことができる。

¥ ︑

EEF/

LH 

/ /  

︐ K 

︿

v v 引

'bハ

m yd  

/

aE

+ 

K  ιH  E¥ 

﹁1111411li‑‑

﹂

' k ' K  

/ / / /  

ι k ' r

 

︿X︿β円﹁lili‑‑￨﹂一一

﹁l

il i‑

‑1 14

' k k  

ノ ノ ノ / ' bH

ハ

ι T AX

︿β円﹁

li li

‑‑

11_﹂

( 417) 

﹁L ノ

R  ‑LK  /

¥1 11 11 11 1/

ノいγk一_h

vk

C

一

μ

hい︑

︐ 町 一

ri

t' K

￨

C J一

μ

/ I l l i ‑

‑ ‑

¥  

k 

ノ /

lk 

p i  

‑ ‑

'k  

k  ( 418) 

ここに，次のベクトル及びスカラーを定義した。

T 4 

fl

it

‑R  

C 

三 RJl/2

白

^{klk‑l (419)} 

c y ) T  

ペ

^J)三

^{cjj)TPMA 九 1}

⁽⁴²⁰⁾ 

77 

，

・ー

また，以上の諸量を用いれば，推定分散・共分散 P の観測更新も次式に 基づいて逆行列不要となる。

LH // 

gu

n 

PA  

¥1Illi‑

‑‑

ノT一

︐FtE︑︐

A M H

‑︐ ︐

FC7tk 

lリ一

μ

︐1

17

C一 一K

刀 工 戸

/1

11 11 1l

¥ 

‑ U H  

//'' 

I HH  

D‑

‑M H  / /  

'h八

Di 一 一  

ιk 〆/LK 

PA

  (421) 

このとき，(421)式の P にU‑D分解(対角要素がすべて lの上三角行列 U と対角行列 D の積)を適用する。さらに， Pkk について ，Pkk‑1 の U 成 分を用い，次式に分解して考える。

PK/k = U K A/kuf/k=UK/kFAUf)ULK‑l  ⁽⁴²²⁾ 

このとき， P の正定性より U は正則となるから， (422)式右辺括弧内に ついて， (421)式より次の関係が得られる。

' bn  

︐K D 

︐ K  ' K 

U 

¥I ll i‑

‑‑

ノ

中J一

iL︑k

寸川 リ

C71パ

‑ w

j

一 /

Ilk

zK

一ρL一

刀 工

戸

/I

ll l1 1l

l ¥ 

gK   Ti zH  

u 

lk ︐k D 

'k  

k D 

TiιMH  一 一

U  

' bn  

D  

FK  

U

  ⁽⁴²³⁾ 

すなわち，P_k1k^‑₁の U‑D分解 Uk1k‑1，_{D k1k‑1} が既知であれば ，(423)式 の右辺の計算を実施した後，これを U‑D分解して U_k，D k を求めれば，以 下のとおり9pk/k を算出できる。

Uk1k = _Uklk-lUk~ D klk = D k   (424) 

4.4.2時間更新過程

本節では時間更新前(添え字 : k / k )の諸量から，時点 k のプラント ダイナミックスに基づく観測更新後(添え字 : k+l/k )の諸量の安定な算

78 

‑ ー

出法を整理する^c

ます¥状態、方程式(409)につき，状態遷移行列の算出に係わる(408)式を考慮、 すれば，求める状態量パラメータの時間更新は次式の連立一次方程式として 逆行列を含まない演算に帰着できる。

¥lil11ノ

I K L K   II

︐JJ 

‑ k t  

︿X︿θ d

/lill‑¥ ¥111

11

/ 1

' bn  

# 汁

+  U一2 A  ' K

Ti  

/Ill‑‑¥ 

一 一

¥111

11 1/

ノ

' b n v K   / J

︐rf 

ム一

F K ' H  

︿X

/ Aθt

Il ll 11

L¥ 

¥I

ll i‑

‑/

ノ

' K 

4汁yKA U一2

Ti  

/a

ll l‑

¥  ﹁11111111﹂

¥1

Il li

‑‑

ノFK'un 

︐ ︐

︐

︐︐

︐︐

︐  lh

ハLH

︿X︿θd

/ l l l l

¥  

k ︐ ︐ ︐ ︐  

がKA 

¥1

il il

‑‑ ノ︑︑_B

BE E‑

ノF

'b八

/ /  

︿ ハ 的

し民

u o  

ιMH 

︐ ︐ ︐ ︐ ︐ ︐ ︐  

k 

︿X〆ItZEE‑1︑σ﹄U/

fi ll

‑‑

‑¥  

﹁I

ll 11 11 11

1_﹂

LMH 

a︐

︐

︐

︐

︐  

# k  

B 

4F

ι 

A + 

(425) 

また，

P

の時間更新については， 5.2節の諸式で求めた

F

_k/k，Q k/k  に 基づき，前節と同様に U‑D分解形式 (P の正定対称性の維持のため)で次 式を演算する。

PMK=LJMFCK+OK/k  ⁽⁴²⁶⁾ 

具体的な手順として，まず¥次式を考える。

¥ー

ーー ノ

/1

J K  

ny  

‑‑︑︐K nu t 

/rllL¥ 

' K 

一D

アA

T A

hi k

.

/V K 

q q  

' K 一 一

 

/〆

+ 

TJ 'K  

U 

LH 

/ /  

十k D 

LK 

/ /  

+ 

ιMH U 

(427) 

ここに，次式の諸量を定義すれば， (427)式 は(426)式と同値になる。

¥l

il

‑‑

ノ

O I  

LK UO D 

/ell‑‑¥ 

一

D

ァJhいkn H S  

︑︑

. .  

EE

‑z

︐ ︐

︐

︐J 2川

1L 'H  

︿ ︒

V1

Lk y

K 

U 

Eり 一 一 一

(428) 

ァi

IHI︐k 

ny  

79 

~

従 っ て ， 以 下 ， (427)式 右 辺 の U‑D 分 解j去 を 示 す 。 これには，

h 川 ^{( j =} 

^{1:  のからグラム=シュ ミット法により}

^め

¹⁾

⁼  q ~l ) を起点に

直交 化 基 吋i)

べ }

^j=l

， . . n )

^{を作成する} これを用いれば当該直交化過程 の変換行列は対各要素が lの三角行列であり，以下のとおり，該手順は P の 時間更新における U‑D分解形式算出である。

TJ ︑s a ' a

''flk 

n r  

flil‑‑

¥1・Ill‑

Tj hいk

nH a 

= 

U 

k+l/k 

I : 

(429) 

TJ ︑₁

fv J1

k   ny  

t﹄

B1 E'

ht︑

iE EE EE

︐ ︐ ︐ ︐   E

TJ 

/H γι M H   nu a 

ここに

い 川 ^{( j}

^= 1， 

^{η )}

の直交性より次式が成立するから 問 式 (430)式で与えられる Uk+l/k'  D kJ‑l/k が(427)式の U‑D分解に一致する。

p(l)T  

k 、

iBKMl) pjyI))  p(fl)T l ¥ ノ

k 

= d

叫

^djl)

ょ う

^=D

(430) 

最後に得られた状態量，パラメータから操作量を求めれば制御系が完結する。 すなわち，次章で述べるプラント制御手順は次式に相当する。

Uk =

市 川 ^， e

^{k!k)  (431 )} 

80 

~

4.  5 

まとめ

石炭粉砕機，火炉の両制御サブシステムにおける同定・状態推定の枠組みと して次の特徴を持つアルゴリズムを求めた。

1)拡張カルマンフィルタに基づき，物理現象の記述が容易な連続時間の動特 性モデル，及び，石炭性状等の未知パラメータの逐次変動を仮定

2)当該拡張カルマンフィルタ演算に，

UD

分 解 ア ル ゴ リ ズ ム に よ る 推 定 分 散‑共分散行列の非負定値対称性の確保と， Pade(l，l)近似によるプラントモ デルの絶対安定な時間離散化とを用い，長期のオンライン稼動における安定 な数値演算に配慮、

81 

，・ー

5 . システム構成と性能検証

5. 1はじめに

本章では，粉砕性，燃焼速度，灰の含有量等，石炭性状の影響による当該石 炭粉砕機出口微粉炭流量と火炉出口ガス温度への外乱に対処するため，

r

新開 発の微粉炭粒度分布パラメトライズ法を用いた非線形動特性モデル，これを用 いた石炭粉砕機制御サブシステム」と「炉内送出波形と信号処理法を改良した 音響式ガス温度計これを用いた火炉制御サブシステム」の構成を具体化する。

さらに，両サブシステムにおける適応状態推定機能の動作検証結果，および，

これらを適用したボイラ制御システムの大型発電用ボイラでの稼動実績の妥

w性についても紹介する。

5.  2石炭粉砕機制御サブシステム

5.2. 1 構成と機能

石炭粉砕機制御サブシステムはFig.501の構成を用い，動特性モデルを主体 に未知パラメータ同定・状態推定によりオンライン計測困難な諸量を把握し，

これに基づく最適な制御を実現する。なお，本サブシステムにおける状態推 定・パラメータ同定演算の詳細は4章に述べてあり，本節ではフラント制御へ の応用的側面を述べる。

まず，状態推定のステップとして実機操作量u(給炭量，分級機回転数，粉 砕部加圧力，搬送空気量)の計測値および特性パラメータf)(漸次に変動する 量;石炭粉砕性，粉砕部摩耗)の仮定値を粉砕機動特性モデルに与え，観測量 y (オンライン計測可能な量;粉砕機差圧l，粉砕動力)の予測値および状態 量x (オンライン計測困難な量;粉砕機出口微粉炭流量，保有炭量，両者の粒 度分布)の算出値を求める。

1 t繍嘱託:摘控~入口と杉帰耕母国間碍証て当3務相紳嚇材茶有量と併目防漏られてし 1るこ

82 

‑ ー

pulverizer  manipulator 

burner  fuel  demand 

inlet coal  classifier  rev  grinding press  etc 

etc 

ulverizer  dynamics  model 

d x二 f(x，u，θ)  d t 

V二 h(x ^，θ )   _pulverizer 

︿UYJ︐+ L   1 d   s n e  

同

c d p r

・‑o︐d 

f t e  

‑l o r   d m p   'roller abrasion 

F i g . 5 0 1   P u l v e r i z e r  Controlling Subsystem 

このとき，特性パラメータ Oの現在値は不明であるから最善の値を同定する 必要があり， yの予測偏差(計測値とその予測値の差)を低減する方向に新た

な仮定値 Oを探索する。すなわち当該偏差最小(実フラントの挙動を最も的確 に予測した場合)を与える Oの仮定値，および，その際の動特性モデルによる xの算出値を，それぞれ求める Gの同定値， X の推定値(現時点で最も信頼で きる値)とみなす。よって適当な時間間隔で以上の手順を繰りかえし，連続的 にOの同定値およびxの推定値が得られる。

次に制御のステップとして状態推定ステッフで求めたハラメータ Gや状態 量xと，その目標値 (火炉制御サブシステムが算出する各段パーナ入口の微粉 炭流量目標値，保有炭量，粒度分布の基準値)との偏差から次の時刻の操作量

uを求めるc

すなわち石炭粉砕機制御サブシステムにおける状態推定および制御のステ ップを通じ，各パーナ入口微粉炭流量目標値への追従および保有炭量，粒度分 布の適正化(粉砕機の良好な運転状態系的)を実現できる。さらに実機特性に

83 

，・ー

追従 (状態推定によるパラメータ 8，状態量xを反映)した粉砕機動特性モデ ルを用い，し¥わゆる予測制御によるプラント負荷変化中の石炭粉砕機台数増減 時期の最適化(各粉砕機出口微粉炭流量将来値とオーバー/アンダーファイア

リンク勺こよる所要微粉炭流量との偏差を解消)を行える。

5.2.2粉砕機動特性モデル

石炭粉砕機制御サブシステムでは動特性モデルのオンライン稼働を要し解 析精度向上と並行して計算量低減を求められる。すなわち Fig.502に示す粉砕 機内諸過程の特性(粉体流量変化と粒度分布変化は相互に依存)に基づき，動 特性モデルは混合，粉砕，分級過程の粒度分布変化を忠実かつ低計算量で模擬 する必要がある。

turntable  plverizing 

F i g . 5 0 2   P r o c e s s e s  i n   P u l v e r i z e r s  

2オーバー/アンダーファイアリング :ボイラ負荷上昇/下降時において伝熱管メタル熱容量 による遅れに対処するため.過渡的に燃料をそれぞれ過大/過少に投入するこ とこ

84 

‑ ー

をFig.503に示 このため，本サブシステムでは，粒度分布密度関数

g ( c )

K4 ) 

σ，非対称 性 K3'裾の重さ μ，分散

す4次までのキュムラント(平均

で記述する高速版粉砕機動特性モデルを採用している。本モデルの導出等につ いては，基礎的事項を含めて2章に詳述しであり，以下，要点のみを述べる。

高速版粉砕機動特性モデ、ノレではFig.502の諸過程の粒度分布変化をキ まず，

4変数に集約して算出している。本 モ デ ュムラント(またはsen1i‑invariants) 

ルは石炭粒度分布の4次以下のモーメント(経験的に 5次以上のモーメントの 影響は小さく省略可能)を正確に模擬でき，従来法の粉砕機動特性モデル(粒 度分布を 30次程度の区間分割で記述)に比して 7分の lの次数であって計算

(計算量は，概ね，粒度分布の記述次数の 2乗に比例する

，> 

particle size 

D i s t r i b u t i o n  D e s c r i b e d  w i t h  s e m i ‑ i n v a r i a n t s  

85 

￨11  mean  1 

〆

~I ヰ ~kewness

¹ 

国

AH hH  

一ω c ω

℃

kF H

一 一 一

﹄ の

﹄O

﹂ 門 ]

時間は約 50分の l となる。

F i g . 5 0 3  

ため)

~

ただし，密度関数

g ( c )

概形の具体的表示を要する際は， 4次までのキユ ムラントについて次のエッジワース展開を用いる。

ドキュメント内 発電用石炭だきボイラの制御のための適応状態推定 に関する研究 (ページ 79-92)