階乗！から三角関数が顔をだす

　最近であった数学の話題です。

ご存じ階乗 2! =2,　3!=6, 4!= 24 　・・・となりますが、 n ! = (n-1)! × n の関係があります。　ではこの n を有理数・実数などに拡張し、　1.5 ! の計算はできないものでしょうか？

　この計算を行うのが、Γ（ガンマ）関数。平たくいうと、自然数のみだった「階乗」を拡張したものです。　Γ(1) ,Γ(2),Γ(3) ・・・が階乗に相当するのですのですが、正しくは Γ(n+1) = n! となります。　

　階乗を拡張したΓ関数ですが、　別の計算をすると、三角関数が表れます。その式は

　Γ(x)Γ(1－x)=π/sin(πx)　というシンプルなものです。

　ζ（ゼータ）関数というのも、興味をそそり、

 1+2+3+4+5 ・・・ =∞ではなく、 = - 1/12 となるのです。下記の書籍での出会い。

（参考）　へんな数式美術館 ~世界を表すミョーな数式の数々~ (知りたい!サイエンス)

　　　　　https://ja.wikipedia.org/wiki/1%2B2%2B3%2B4%2B…

---

 今取り組もうとしているのは、「多重三角関数」。私の研究成果である、Synchro PRIMO法を複数成分に拡張しようとすると三角関数の展開は超難解。また、自身の提案した方法の「一般化」もめざして「パルス」の位相差にも取り組んでいるところ。

　数学のネタを探していたら、Γだのζなど、読めない関数にであい。引き込まれてしまいました。　大学２年で「複素解析」をやっていれば、入門としてはついていけると思います。

 

調波を重ねると・・・奇妙な波形

　一風変わった信号を紹介します。

X(t) = sin ωt + (1/4) * sin ( 4ωt ) + (1/9) * sin (9ωt ) + ... +  (1 / m^2 ) sin ( 1/2^m * ωt)

これは  「リーマン関数」と呼ばれています。　https://ja.wikipedia.org/wiki/リーマン関数　

基本周波数に、4, 9, 16,25 ..... 高調波　が、振幅 1/4, 1/9, 1/16 , 1/25 となって加わっていくだけです。

波形は奇妙な形になり、無限に重ねると「いたるところで微分不能」な「連続関数」になるそうです。このような波形を「病的関数」というらしく、「ワイエルシュトラス関数」　https://ja.wikipedia.org/wiki/ワイエルシュトラス関数　も奇妙な波形になってしまします。

・・・　ちょっと考えると、フーリエ変換の定義の範囲で、特別なことをせず成分を重ねていくと、波形が微分不能になる、つまり、ある特定点での微係数を求めることができないということです。　当然２階微分も求められないわけなので、周波数も定義できない？　ことになるのでしょうか。

　非線形の世界なら、こういうことは起こりえるのかもしれません。音で考えるなら、リーマン関数のルールで調波を重ねていくと、どんな音なのでしょうか。 興味あるかたは実験してみてください。

（追記）　2020.1.23

　前述のリーマン関数を使い、逆フーリエ変換すると、「あちこちで微分できない」信号ができてしまいます。従来からの伝統的な考えでは、既知の「周波数成分」をもつ、ごく普通の波形でしょうが、微分できないとなると、その波形はどう解釈すべきでしょうか。　関数の連続性と微分可能性については、フーリエ変換は、あまり明快な答えを示していません。　一度、整理をしてみたいと思います。

　

「基本周波数」定義のジレンマ

　２つの、いわゆる「周波数成分」からなる波形の「基本周波数」は、どう定義すればいいでしょうか。例として、100 Hz の成分１、300 Hzの成分２　を考えます。　

それぞれの振幅、 100Hz成分振幅＝A1,　300Hz成分振幅＝ A2　としましょう。

（１）　A1 >> A2 の場合

　常識的に、300 Hzが "ちょっとだけ"混入しているが、基本周波数 = 100 Hz　と考えてよさそうです。300 Hz の成分２は第３高調波として処理すればなんてことないです。

（２）　A1 << A2 の場合

　前と動揺に、 基本周波数＝300 Hz と考えてよさそうです。1/3 の分数波が混じります

（疑問）

・　振幅A1, A2 がどのような条件の時に、上記（１）（２）の境界となるのでしょうか？

・ A1 = A2 のときの基本周波数はどう考えるべきでしょうか。三角関数の和積変換から、周波数が 200 Hz であると見ることもできます。

・　瞬時位相角から導出される瞬時周波数は元の成分以外の値をとり得ます。基本周波数とは、これらの平均周波数と一致するものなのでしょうか。

・・・　さて、どう解釈すべきでしょうか。基本周波数はどう定義すべきでしょうか。

 

電力系統周波数の「急降下」（続報）

　ちょっと前に投稿した記事、電力系統周波数の「急降下？」 の続報です。

　Rapid Mains Frequency Decsent in Tokyo   https://youtu.be/ejQsBP8EW80

 

　データを解析してみたところ、周波数変化率・変化の幅、周波数低下してからの挙動がどうも発電機の「系統からの解列」に思えてなりません。そこで思い切って東京電力に電話。（注：最近の用語では電凸て言うらしい）

　最初に電話した部署はカスタマーセンターみたいなところ。事情を丁寧に話したら内容は伝わったみたいで、あとで専門の部署からご返事しますとのこと。なにせ発電設備（または送電）に関する事項なので簡単には教えてくれないだろうと思っていましたが、なんと、大船制御所の方から直接ご連絡いただきました。観測した事象を詳細に説明したところ、詳細を調べていただき、後日わざわざ電話でご連絡をいただきました。丁寧なご対応にこの場を借りてお礼申し上げます。

　結果は、予想どおりで、私の観測した時間帯、"近隣の"火力発電所で「ちょっとしたこと」があったらしく、私の想定した内容どおりでした。もっと細かい話を伺いたかったのですが、運用のプロに対して素人が立ち入る話ではないと思い、丁寧にお礼を述べ電話を切りました。　

---

 　話変わって九州電力での話。

こんな「フリッカー」が九州電力管内で広域に発生したみたいです。

http://www.sankei.com/region/news/170503/rgn1705030054-n1.html

　原因は太陽光発電についている「パワーコンディショナー」（以下パワコン）（=交流をつくって系統に流し込む装置）が「無効電力」を大量に流し込んだためと判明。電力会社は、実際に消費される「有効電力」、電圧の維持や系統の安定のために「無効電力」の調整を行っています。東京電力の例では新富士変電所（富士サーキット近く）に系統安定化のための「同期調相機」が置かれています。ここで無効電力のコントロールが行われています。

　九州電力は系統容量が小さく、隣の中国電力と連系している連系線が細いため、太陽光発電が増えると何かしらの問題があるかも、と以前から感じておりました。　2015年3月に別府で研究会に参加した際もAC100Vモニターと近隣の変電所（西大分変電所 ＜地図＞）を視察しています。北陸電力とくら朝7時〜9時頃の変動のパターンが特徴的だった印象があります。

　この問題の本質は、パワコンの「単独運転検出機能」が未熟だったこと。少数のパワコンだけでは問題なくても大量のパワコンが同じ検出動作(アクティブ方式）を行い系統にダメージを与えてしまったこと。九州電力ではこんな「お願い」をしていますが、単価の高い電力を買い取らされ、今回のような系統が不安定になりかけて事業者に「設定変更をお願い」はおかしな話に感じます。系統に影響を及ぼす装置に対しては、「系統側の判断で遮断する」くらいの気概で望まないとパワコンも進化しないでしょう。

　電気工学的にみた「理想的なパワコン」は・・・　三相同期発電機とまったく同じふるまいをする発電機だと思います。近隣の発電機と周波数・位相同期を行うことができ、アンバランス時は「横流」で負荷分散する・・・　周囲の電圧が上昇してきたなら（無効電力の増大）自身は誘導性負荷に変身すればいいのです。　安定供給・系統安定化する機能をもっている系統側の発電所は遠慮する必要はなく、分散電力側は、電力を供給するからには、最低限の有効電力・無効電力をきっちりと制御して欲しいものです。ちなみに「遅れ無効電力の消費＝進み無効電力の供給」といったはて？な関係があります。

　この「有効電力(real power)」「無効電力(reactive power)」は電気工学を勉強していないと馴染みのない用語です。スマートグリッドExpo で説明員に意地悪な質問をしてみましたが、無効電力が系統安定のために重要な働きをしていることを理解できている方はほとんどいらっしゃいませんでした。。。交流は「送っているうちに電圧が上昇」することがあるという摩訶不思議なことだらけなのです（フェランチ効果といいます）　情報工学全盛の時代、コイルだの電線だのというレトロな？技術は興味が失われているのでしょうか・・・　

　三相同期発電機と同等の特性をもつパワコンが開発されれば世界的な表彰ものだと思うのです。今の時代、「基準発電所」からのフェーザー情報をパワコンが受信し、所定の制御指示にしたがって、所定の電圧・力率（＝つまり所定の有効／無効電力）で発電するなど簡単だと思うのです。そういう装置にこそ、高い単価で買電してあげれば良いと思います。

瞬時周波数の定義でおこるジレンマ

　 まったく同じ振幅で異なる周波数の２つの信号を混合した際の「瞬時周波数」はどう計算すべきか。

　三角関数の和積の公式から sinθ1 + sinθ2 =sin (θ1+θ2)/2  * cos(θ1-θ2)/2 　となるのですが、これは、「２成分の信号」と「１成分に振幅変調された信号」は、本質的に区別できないこと意味しているのではないでしょうか。

　つまり、瞬時周波数と瞬時振幅を正確に測定しようとしても、このジレンマの答えを整理することが先決だと考えています。実際、振幅が同一で周波数が異なる２成分の瞬時周波数はややこしい形をした公式になっています（解析信号で考えると導出できます。面倒です・・・）　しかしこのアプローチでは、「位相角が微分できない箇所」が発生してしまいます。　ここでもう一歩踏み込むと・・・「１サイクルでの平均値」は計算でき、これは周波数１、周波数２の平均値になり、前述の「和積変換」と同じ答えがでてきます。

　つまり、Fourier変換で複数成分の和と考えているかぎり、上記のジレンマには気づきません。最近はやった解析信号で得られる「瞬時位相角θ」の微分から ω=dθ/dt　とやってもやっかいな問題に出くわしてしますのです。基本周波数F0は、瞬時周波数を「１周期」で平均して計算できるなら、既知の常識と一致します。しかし、1/2 周期の信号が交じると、解析信号の軌跡（閉曲線）が２周で閉じる形となり、「どこまでが１周期？」という問題が発生します。

　このような瞬時周波数に対する疑問・考察は 1990年後半から2010年あたりまでIEEEでちらほら見かけましたが最近はめっきり。

　もしやと思って「三角関数」の発展の歴史を調べてみたのですが、チェビシェフ級数でさえ １９世紀末。1900年をすぎて電子回路技術と理論が急成長。その間に、その次代で知られていたさまざまな手法が導入されたのでしょう。想像ですがその１つがリサージュ図。図形の外形から位相差がわかる便利なものでししたが、その閉曲線に着目した数理は登場していません。　その理由も考えてみたのですが、電磁気でならう「グリーンの定理・ストークスの定理」で面積分・閉曲線の積分がひとつの傍証。ストークスはマックスウェルと同時代の人で、マックスウェルの良き相談相手だったらしい。なので体系がまとまってくるのは19C末だったでしょう。

（追記）　上記の考えは、D論を書くときから気付いており、"Die Singularitaten Der Lissajous'schen Stimmgabelcurven" (私訳：音叉のような形状をとるリサージュの特異点）という古書を購入、読んでみました。x(t) = sin (at+ b) のように「周波数項」と「初期位相項」にわける考えはありますが、Lissajous軌跡の形状の分類にページを割いており。当時の数学では「閉曲線」や「閉曲線面積」に着目するには無理があったような印象です。

　https://www.bookdepository.com/Die-Singularitaten-Der-Lissajousschen-Stimmgabelcurven-1875-Wilhelm-Braun/9781168311795

　こう考えると、（私が試みた）ある軌跡の「面積」「軌跡からの面積算出」というのは、19C後半から20Cにかけて可能になってきたのであって、オシログラフの登場した1920年にはほとんど知られていなかったことになります。　その後・・・もし・・・位相差測定を時間差測定で代用できる技術が一般化していたなら、位相差計算の数理をまともにあつかった研究が無かったとしても不思議ではありません。

　現代の研究者・エンジニアは先人たちの考えた手法を身に着けています。しかし、時代背景や技術開発のながれに取り残された分野はどこかにあるはずです。もしあるとしたら、「みんな知っていること」ではなく常に疑問をもつ・・・心がけで真理をみつける心構えが必要だと思うのです。

　もうひとつ苦言を呈するなら、FFTの乱用。任意の波形を正弦波の和として考えると実用上の威力は抜群です。しかし、もともとFourier級数は微分方程式を解くための方法であって、その後の数学者達によって現在のFourier積分となり、Freuency-Timeドメインの解析が可能になったのです。Wiener, Gaborのように数式だけでやっていればよかったのかもしれませんが、1950年以降、計算機科学で有用性が見直され、1964年? のFFTの登場で一気にブレークします。　しかし元々周波数を測るための道具でなかったFourier Transform がF0推定に使われることには、未だに違和感を覚えます。「基本周期とは何か」という Fundamental な議論が不完全のように私には見えるのです。。。

 　上記の「２成分の瞬時周波数」、三角関数の和積変換からは、「F0は２つの周波数の平均」となりますが、一カ所で位相角φが微分できない点をどう解釈するのか・・・　今後成分を増やしていくと「特異点」がどんどんでてくるでしょう。

 

５点でもとめる周波数計算方法 ・３点で求める位相差計算の動画

　本ブログで発表した Excel による計算例を、動画(HQ)にしたものをYoutubeに掲載しました。後日Anotation で解説を加えます。

５点のデータで周波数を求める手順を動画で制作し、Youtubeに掲載しました。周波数のみです。

https://youtu.be/MMvmpj9JbpU

 

２つの時系列波形の、１）それぞれの周波数、３）それぞれの振幅、３）位相差　の計算手順の動画を公開しました。

https://youtu.be/ZG7wL3z1rHg

　シナリオ無しで、Excelに手打ち操作しながら録画しました。計算式は覚えているのですが時たま勘違いしてNG。。。。一発撮りの録画は大変つかれました。多少もたついているところがあります。

電力系統周波数の「急降下？」

　東京電力の電力波形(50Hz)をコンセントからモニタする実験を続けていますが、溜まりに溜まったデータを解析していて、珍しいシーンを発見しました。

研究会で発表しそこなったので、YouTubeに測定の様子を公開しました。

Rapid Mains Frequency Decsent in Tokyo   https://youtu.be/ejQsBP8EW80

　周波数変化率(ROCOF)が瞬間的に -0.05 Hz/s になっており、Synchro PRIMO の目盛りからすっ飛んでしまいました。数秒で低下は収まり -49.8 Hz で安定したのですが、もし低下が止まらなければ、あっという間に49Hzを割るところだったでしょう。 この回復はどうやったのか？

　電力系統工学では「周波数制御」として解説されています。今回の現象は、「ガバナーフリー運転」が働き、複数発電機の間の「速度調定率」できまる周波数に落ち着いたのだと思います。急激な負荷の変動を、発電機間で「うまく分担」して「ちょっと低めの周波数」で（自立的に！）安定させたのです。

　データを細かくみると、現象の発生する数十秒前に、わずかな「（進み方向への）位相飛び」がありました。推測ですが、どこかの系統で工事等で計画的に切り離したのでしょう。周波数低下のカーブは綺麗なexpカーブだったです・・・

　動画では公開しませんでしたが、この後、数分くらいでの50Hzに何事もなく戻っています。ご安心くださいませ。

　こんな技は三相同期発電機のみなせることで、パワーエレクトロニクス技術満載？のパワーコンディショナーでこんな制御はできるのでしょうか？　分散電源は大ブーム、でもやばくなるとさっさと「解列」して逃げてしまう・・・　「３相同期発電機」は歴史も古くレガシー技術ではありますが、こんな秘めた能力について、皆さんにも知っていただきたいと思います。

　下記はイギリスでのワンシーン。一瞬周波数が -0.3 Hz まで低下しています。フランスからの緊急送電でしのいだみたいです。　Synchro PRIMOが検出したのは、これと同じ現象です。

https://www.youtube.com/watch?v=UTM2Ck6XWHg

 

参考図書

　荒井純一（監修）ほか,  ”基本からわかる電力システム講義ノート" , オーム社.

　

微分と差分

　デジタル信号処理では、微分の代わりに差分が用いられ、アナログでは s をつかったラプラス変換、デジタルでは z 変換がいわゆる周波数領域での設計に用いられています。

　z 変換に登場する z^(-1) はデジタル処理では、「一段の遅延」と等価で、この原理が Synchro PRIMO では大きな意味をもっています。

　２階微分を考えてみます。　x(t) = A * sin(ωt)　という連続関数では ω 自乗が飛び出してきて、d2/dt2 x(t) = - (ω^2) sin(ωt) となりますが、「離散系」（＝デジタル）ではどうなるか？

　このあたり、詳しく解説されている専門書・教科書は少ない気がしますが・・・差分演算は連続系の問題を解く際の「近似解法」と考えられているような気配も・・・

計算は難しくありません。一応、「１階差分」をとり、差分の差分として２階差分を計算すると、

Δ(2) x[n] = 1/4 * { x[n+1] + x[n-1] - 2* x[n] } * A

となります。Aは振幅項です。連続系では  - ω^2 となりますが、デジタルではどうなるか。　z 変換 を習ったあとの演習問題としてはよくできていますね。半角の公式が見事にヒットし綺麗にまとまり、

 A = - sin^2(Ω/2) 　　

となります。　

上の式をみたとき・・・積分ではどうなる？　・・・　との疑問から、積分でやってみましたが、「積分」はデジタルでは、どうなるのでしょうか？どう呼べばいいのでしょうか。差分がΔなら、いっそのこと「累積」=Σ と呼んでみることにします。暫定的にやってみると、

Σ(2) x[n] = 1/4 * { x[n+1] + x[n-1] + 2* x[n] } * A

形としては初学の時に習う「平均化」です。 A を計算してみると、今度は、

 A = cos^2(Ω/2)

となります。あれれ？　Δ(2)とΣ(2) を演算子としてみると、振幅項だけでは、　Δ(2) + Σ(2) = 1 ではないですか・・・微分と積分（平均）に着目した演算子では、３項の和差の１箇所しか違わないのに。。。なにかしらの「秘密の恒等式」が隠れている？？？

　これらの演算、振幅項を厳密に求めていることと、一応線形演算なので、私の 「Lissajous 外積」よりははるかにシンプル。積分とはいっても、非常に短い区間での、（一応）「瞬時演算」なのでエンジニアリング向けに何かの用途がありそうです。

　平均化・スムージングといった「実用演算」ではなく、「微分の代用」といった近似演算でもなく、時系列波形をもっと厳密に考えてみたいと思っていました。　下記文献の影響です。もっと理解したいですが、超難解 です ... orz 　

　参考書籍：　

• 飯島泰蔵,  "自然観測法 - 瞬時性に着目した新しい波形解析法,"  森北出版.
• 飯島泰蔵,  "デジタル自然観測法 - 時系列解析のための新しい理論,"  森北出版. 

 PS: 計算式に誤りがあり修正しました。 (2/11)

欧州の電力系統、各地で波形収集

　欧州は巨大な単一系統(50Hz)となっています。チェコあたりからスペインの先まで、１つの系統になっています。地図でみると線路のように１つに繋がっていますが、すべてが単一系統として運用されているかどうかは不明です。

　系統図は ENTSO-E という機関で公開されています。　

https://www.entsoe.eu/map/Pages/default.aspx

　欧州では、電力の売買が自由化されているのですが、安定供給のためにはロードバランスや潮流を適切に管理する運用が大事なのです。国によっては発電過剰だったり消費過剰だったり・・・その埋め合わせのために「国際連系線」とよばれる国をまたぐ送電線で電力のやりとり（売買）を行っています。とは言うものの遠距離で売買をしようとしても実際の電力を送るのではなく、「最終的な供給・消費の収支」が合うように発電量・送電量が調整されるのです。上記の地図をみると国内はメッシュ状に多数の送電線がありますが、国をまたぐ送電線はそれほど多くないことが分かります。

　2015年、国際会議での発表がてら、欧州各地をまわり、電力波形・周波数変動データを収集してきました。しかし、その場所の選定をどうするか。。。そこで、上記のENTSO-Eの巨大地図のお世話になっています。

　この系統図から選んだ地点（訪問日程順）

1. スイス：Laufenburg ・・・スイスのオペレータ(Swiss Grid)の本拠地であり、太い回線が多数集結、ドイツとのやりとりをしている。近場に原発もあり周囲の様子をみたかった。
2. ドイツ：Simbach/Braunau am Inn ・・・ドイツ-オーストリアを結ぶ国際連携線。Inn川をはさんでドイツ・オーストリアの変電所が接続
3. オーストリア：Villach ・・・かつて仕事で住んだところ。周囲はほとんどが水力発電で、大消費地からも遠い。Pureな波形ではないかと予想
4. オーストリア：Wien郊外　・・・東欧諸国と接続される巨大変電所。一帯にはLNG発電所が多数。消費地に近い
5. ドイツ：Aachen　・・・　ドイツ・オランダ国境近く。波形はあれているかもと予想
6. ベルギー：ベルギー・フランス間にある連携線　・・・　帰国直前の静養もかねてベルギーのDinant という町を起点に現地を視察。

　このツアーですが・・・欧州で発生した難民の大移動に見事に遭遇してしまい、当初予定を変更、上記のように落ち着いています。当初はチェコ方面も計画していたのですが、難民のハンガリーからオーストリアへの流入、クロアチア経由オーストリアへの流入に遭遇しキャンセル。

 

　上記Simbachという町はInn川を渡ってオーストリアからドイツへ徒歩で移動できます。なんと国境の橋(Simbach-Braunau間)では検問があり救護所の設営の真っ最中でした。またWien中央駅でも1000人を超えているであろう集団に遭遇しています。そんな事情もあり移動は電車・バスのみとし、治安に不安を生じそうなエリアへの移動は控えました。

 

　この時のモニター波形については、次回以降、記事にしたいと思います。

 

 

 

 

 

 

ゆらゆら？ 系統電圧波形

　添付の画像(GIF)をご覧ください。複数画像をあつめてGIFアニメーションにしたものです。

　この波形は、フランス・ニース市内で取得した電圧波形(欧州の系統は50Hz)です。なんと、奇数／偶数サイクルごとに図のようにゆらゆらと波形がゆれています。ミクロ的には、ゼロクロス点が毎サイクルで前後しており、Synchro PRIMOで周波数を計測すると、２値が突然現れるように観測されます。

 

 

　詳細な原因は不明ですが、おそらく１サイクルおきに電流が大きく変動する負荷がつながったのでしょう。しかし抵抗負荷だけだと仮に非線形であったとしても、ゼロクロス点が変動するのは理解できません。おそらく・・・ですが、かなり大きい誘導性の負荷が、かつ１サイクルおきに変動していた・・・と想像しています。

　見かけの周波数のゆれは、別の言い方をすれば、位相角が波打っているともとれます。時変の誘導性負荷が接続されていたと考えれば、「インピーダンスの変化率」できまる「位相角の変化率」が系統周波数の変動として観測されているのでしょう。

　現在、欧米で研究開発がすすめられているＰＭＵ(=Phasor Measurment Unit)ですが、こういった電圧フリッカーの検出はできるのでしょうか？　PMU仕様や観測標準として検討されているのでしょうか・・・　

　　Synchro PRIMOは半サイクルのデータでも瞬時周波数・位相角、瞬時振幅の測定が可能です。

　余談ですが、2015年9月にニースで開催された国際会議(EUSIPCO-2015)で発表した際に計測したものです。ニースの配電系統はわりとシンプルで、ニース郊外のCarrosという町にそこそこの変電所があり、ここから供給されているようです。https://www.google.co.jp/maps/@43.80259,7.1916265,16z　近隣には工場が多数あり、変電所に近いところでの誘導負荷の影響があったのかと想像しています。

　CarrosまではNice市内からレンタカーを使用しましたが・・・現地の運転は短気なドライバー？が多く、市内への戻りは大変だったです・・・Carros から上流にも足を伸ばしましたが・・・ドライブを堪能できました。ただし欧州ではMT車が標準・・・私はMT乗りなので問題なかったですが、シフトミスでエンスト多数。少々後ろをイライラさせたかも。