JP7467376B2 - 制御装置、制御方法およびプログラム - Google Patents

Description

本発明の実施形態は、制御装置、制御方法およびプログラムに関する。

例えば、風力発電、潮流発電、および、海流発電などのエネルギシステム、並びに、インフラシステムの制御装置、または、精密医療向けマイクロプラントにおける流路・ソーティング・反応制御システム、などの制御装置に関して、センシングデータおよび運用データなどのモニタリングデータを用いてシステム制御を行う技術が知られている。また、これらのシステム制御に関する物理量を、物理現象シミュレーション、および、非線形数理モデルに基づくシミュレーションにより推定する技術が知られている。

特許第４４３９５３３号公報

Sam Greydanus et al., "Hamiltonian Neural Networks", in arXiv: 1906.01563 5 Sep 2019． Miles Cranmer et al., "LAGRANGIAN NEURAL NETWORKS", in arXiv: 2003.04630 30 Jul 2020．

しかしながら、従来技術の流体シミュレーションおよび流体・構造連成シミュレーションなどの計算負荷を要する物理現象シミュレーション、および、非線形数理モデルに基づくシミュレーションでは、物理量の推定に多大な時間（例えば数時間以上）を要する場合がある。このため、モニタリングデータを用いて、リアルタイムにシステム制御を行うことが困難である。

実施形態の実施形態の制御装置は、取得部と、算出部と、制御部と、を備える。取得部は、制御対象を示す解析領域を離散化した複数の要素ごとの第１時刻の第１物理量をそれぞれ表す複数の入力データを取得する。算出部は、複数の入力データを推定モデルに入力して得られる、複数の要素のエネルギを表すエネルギ汎関数の値に基づいて、第１時刻より後の第２時刻の複数の要素ごとの第２物理量を算出する。制御部は、第２物理量に基づく制御量が目標値となるように制御対象を制御する。

実施形態にかかる制御装置のブロック図。 解析領域に含まれる要素の例を示す図。 実施形態における制御処理のフローチャート。 ＦＥＭＬＮを用いた解析の概念を説明するための図。 ＦＥＭＬＮの構成例を示す図。 ＦＥＭＬＮの構成例を示す図。 損失関数で用いられる各勾配の差を説明するための図。 実施形態における学習処理のフローチャート。 実施形態の適用例を示す図。 実施形態の適用例を示す図。 制御部による制御処理の例を説明するための図。 自励振動する振動物体の例を示す図。 自励振動する振動物体の例を示す図。 振動物体表面の変位の計測方法の例を説明するための図。 翼構造を有する振動物体の制御方法の例を示す図。 翼構造を有する振動物体の制御方法の例を示す図。 翼構造を有する振動物体の制御方法の例を示す図。 複数の振動物体の相互作用の例を示す図。 複数の振動物体の相互作用の例を示す図。 複数の振動物体の相互作用の例を示す図。 引き込み現象の誘起例を示す図。 励起力係数の変化の例を示す図。 振動物体の変位と変位速度との関係の例を示す図。 可変速度フィードバック機能の構成例を示す図。 可変速度フィードバック機能の構成例を示す図。 可変速度フィードバック機能の構成例を示す図。 電磁誘導型の振動発電の構成例を示す図。 発電機構の構成例を示す図。 実施形態にかかる制御装置のハードウェア構成図。

以下に添付図面を参照して、この発明にかかる制御装置の好適な実施形態を詳細に説明する。

本実施形態では、サイバーフィジカルシステム（Cyber Physical Systems）におけるデジタルツインを行う際に、センシングデータ（計測、センシングまたはモニタリングしたデータ）から、対象とする系（対象システム）の制御に関する物理量（変位場および変位速度場など）を推定し、制御に活用する手法の例を説明する。以下では変位場および変位速度場を物理量とする例を主に説明する。物理量はこれらに限られず、他の物理場および物理速度場（例えば、圧力場、および、圧力変化速度場、または、ひずみ場、および、ひずみ速度場）であってもよい。

具体的には、偏微分方程式の離散化数値計算手法、および、機械学習手法であるラグランジアン・ニューラルネットワーク（Lagrangian Neural Networks）を導入することにより、システム制御等で必要となるセンシングデータから時間的および空間的な物理量分布への変換を高速化および高精度化する。

離散化数値計算手法は、例えば有限要素法、有限体積法および差分法である。以下では、有限要素法（ＦＥＭ：Finite Element Method）を活用したラグラジアン・ニューラルネットワークモデル（以下、ＦＥＭＬＮ：Finite Element Model－based Lagrangian Networks）を用いる例を主に説明する。

図１は、本実施形態にかかる制御装置１００の構成の一例を示すブロック図である。図１に示すように、制御装置１００は、取得部１０１と、演算処理部１１０と、制御部１０２と、記憶部１２１と、を備えている。

取得部１０１は、制御装置１００で用いられる各種データを取得する。データの取得方法はどのような方法であってもよいが、例えば、ネットワークを介して接続された外部装置から受信することにより取得する方法、および、記憶媒体に記憶されたデータを読み出して取得する方法などを適用できる。

例えば取得部１０１は、ＦＥＭＬＮに入力する複数の入力データ、および、ＦＥＭＬＮの学習に用いる学習データを取得する。複数の入力データそれぞれは、例えば、システム制御の対象（制御対象）となる解析領域を離散化した複数の要素ごとのある時刻ｔａ（第１時刻）の物理量Ｐｈａ（第１物理量）を表す。

図２は、解析領域に含まれる要素の例を示す図である。各要素は、曲面で囲まれた解析領域２０１を空間的に離散化した情報に相当する。節点は、各要素を構成する多次元空間上（図２では３次元空間）の点である。

以下では、解析領域２０１の各節点の空間的な変位をｕと表す。ｕの微分が変位速度（速度ベクトル）である。以下では、ある変数の上にドットを付すことにより、その変数の微分を表す。例えば変位速度は、ｕの上にドットを付した記号で表される。

本実施形態では、要素ごとにエネルギ汎関数が算出される。例えばエネルギ汎関数は、以下の（１）式のように、蓄積エネルギ、損失エネルギ、および、仕事量により表される。

σは相当応力、εは相当ひずみ、Ｆは境界に作用する外力ベクトル、Ｖは物体の体積、Ｓは表面積を示す。第２項は、相当応力と、相当ひずみ速度の増分と、の積分に相当する。第３項は、物体が外力に対してなす仕事（外力ベクトルと速度との積）に相当する。

図１に戻り、取得部１０１は、システム制御に用いるセンシングデータを入力データとしてさらに取得してもよい。

センシングデータは、例えば、物体周辺媒体の流速、流れの向き、圧力、荷重、温度、加速度、変位、電流、電圧、振動、および、ひずみなどを示すデータである。これらのセンシングデータを検出するセンサや上記センシングデータの物理量に変換可能な画像取得装置は、例えば、解析対象とする構造または周辺構造の予め定められた個数のサンプル点に配置される。

センシングデータは、解析対象とする構造のパフォーマンス特性であってもよい。例えば解析対象が風力発電プラントの場合、パフォーマンス特性は、以下のようなデータである。
・発電機性能
・インバータ性能
・回転数
・バッテリ性能
・発電機負荷率
・インバータ負荷率
・電力伝送特性
・信号伝送特性
・ノイズ特性
・冷却性能
・運転履歴

パフォーマンス特性は、例えば、Bios（Basic input/output system）等とのコミュニケーションを行うプロファイリングツールまたはモニタリングツール（監視ツール）により取得することができる。

センシングデータは、負荷・境界条件に相当する。取得部１０１は、センシングデータとともに、解析対象とする構造の設計変数を取得してもよい。設計変数は、例えば、境界条件、材料特性、および、構造変数（構造の形状および寸法など）などを示す情報である。以下では、センシングデータおよび設計変数を条件データという。

演算処理部１１０は、取得部１０１により取得された入力データに対して、ＦＥＭＬＮを用いた演算を行い、システム制御に関する物理量を推定する。演算処理部１１０は、算出部１１１と、学習部１１２と、を備える。

算出部１１１は、取得された複数の入力データから、システム制御に用いることができる物理量を算出する。例えば算出部１１１は、複数の入力データを推定モデルに入力し、推定モデルが出力する出力データを算出する。出力データは、解析領域を離散化した複数の要素のエネルギを表すエネルギ汎関数の推定値である。推定モデルは、統計モデル、確率モデル、および、機械学習モデルとして構成することができる。推定モデルは、ＦＥＭＬＮなどのニューラルネットワークモデルに限られず、階層ベイズモデルなどであってもよい。

算出部１１１は、さらに、出力データを用いて、時刻ｔａより後の時刻ｔｂ（第２時刻）の複数の要素ごとの物理量Ｐｈｂ（第２物理量）を算出する。時刻ｔｂは、例えば時刻ｔａの次の時間ステップに相当する時刻である。時間ステップをΔｔとすると、ｔｂ＝ｔａ＋Δｔと表すことができる。

学習部１１２は、演算処理部１１０が用いる推定モデルを学習する。学習部１１２は、取得部１０１により取得された学習データを用いて、ＦＥＭＬＮを学習する。学習部１１２は、例えば、出力データの勾配と、正解データである勾配との差を最小化するように推定モデルを学習する。

制御部１０２は、制御対象を制御する。例えば制御部１０２は、算出部１１１により算出された物理量Ｐｈｂに基づく制御量が目標値となるように制御対象を制御する。制御部１０２による制御処理の詳細は後述する。

上記各部（取得部１０１、演算処理部１１０、および、制御部１０２）は、例えば、１または複数のプロセッサにより実現される。例えば上記各部は、ＣＰＵ（Central Processing Unit）などのプロセッサにプログラムを実行させること、すなわちソフトウェアにより実現してもよい。上記各部は、専用のＩＣ（Integrated Circuit）などのプロセッサ、すなわちハードウェアにより実現してもよい。上記各部は、ソフトウェアおよびハードウェアを併用して実現してもよい。複数のプロセッサを用いる場合、各プロセッサは、各部のうち１つを実現してもよいし、各部のうち２以上を実現してもよい。

記憶部１２１は、制御装置１００による各種処理で用いられる各種データを記憶する。例えば記憶部１２１は、取得部１０１により取得された入力データ、演算処理部１１０による演算結果などを記憶する。

記憶部１２１は、フラッシュメモリ、メモリカード、ＲＡＭ（Random Access Memory）、ＨＤＤ（Hard Disk Drive）、および、光ディスクなどの一般的に利用されているあらゆる記憶媒体により構成することができる。

次に、このように構成された本実施形態にかかる制御装置１００による制御処理について説明する。図３は、本実施形態における制御処理の一例を示すフローチャートである。

取得部１０１は、解析領域を離散化した複数の要素ごとの入力データ（時刻ｔａの物理量Ｐｈａを取得する（ステップＳ１０１）。算出部１１１は、入力データを推定モデルに入力し、時刻ｔｂの物理量Ｐｈｂを算出する（ステップＳ１０２）。制御部１０２は、物理量Ｐｈｂから、制御対象を制御するための制御量を算出する（ステップＳ１０３）。制御部１０２は、算出した制御量を用いて制御対象を制御し（ステップＳ１０４）。制御処理を終了する。

次に、ＦＥＭＬＮの詳細について説明する。図４は、ＦＥＭＬＮを用いた解析の概念を説明するための図である。ＦＥＭＬＮの出力データであるエネルギ汎関数は、変分原理とのアナロジーから、例えば連続体力学シミュレーション（非弾性応力シミュレーションなど）の場合には、離散化した各要素における蓄積エネルギ、損失エネルギ、および、仕事量により表される。非弾性応力シミュレーションの変位場および変位速度場（ひずみ場およびひずみ速度場）は、変分原理により、これらのエネルギ汎関数の総和が極小となる経路（動的可容速度場）として決まる。図４の矢印は、初期状態から、エネルギ汎関数の総和が極小となるように決定される経路の例を示す。

ＦＥＭＬＮは、以下のように構築できる。
（Ｓ１）ＦＥＭＬＮの定義：
・時間的および空間的に離散化された変位場および変位速度場をモデルへの入力データとする。このとき、条件データを入力データとして加えてもよい。モデルの出力データは、対象システムにおける各要素のエネルギ汎関数とする。このとき、次の時間ステップ（次の時刻）の変位場および変位速度場を出力として加えてもよい。
・入力データから出力データへ変換するモデル（変換モデル）、または、入力データから出力データを推定するモデル（推定モデル）を定義する。ニューラルネットワークの場合は、モデルを定義するパラメータは、層数、各層の素子数、および、各素子の活性化関数の構成などが挙げられる。なお、階層ベイズモデルの場合は、モデルを定義するパラメータは、中間層の潜在変数、データ分布、事前分布、および、ハイパーパラメータの構成などが挙げられる。
（Ｓ２）ＦＥＭＬＮの学習データの準備：
・対象システムを記述する偏微分方程式を離散化し、条件データに関する解析条件を設定する。学習用の入力データに対して、変位場、変位速度場、エネルギ汎関数、および、エネルギ汎関数の勾配を、離散化した各要素および各節点について時間ステップごとに数値計算する。数値計算の結果を、ＦＥＭＬＮにおける損失関数を算出するための正解データとして準備する。ここで、条件データに関してパラメータサーベイを行った数値解析結果を正解データとしてもよい。入力データ、および、正解データを含む学習データが学習処理で使用される。
（Ｓ３）ＦＥＭＬＮの学習：
・ＦＥＭＬＮによりエネルギ汎関数近似モデル（パラメータθ）を作成する。
・ＦＥＭＬＮにより、学習データに含まれる入力データに対して、エネルギ汎関数のスカラー値φ_θ、並びに、次の時間ステップの変位場および変位速度場を出力する。φ_θは、パラメータθにより定められるＦＥＭＬＮにより出力されるエネルギ汎関数の値を表す。
・ＦＥＭＬＮを用いて、以下の（２）式および（３）式でそれぞれ表される勾配を算出する。

・ＦＥＭＬＮを用いて次の時間ステップの変位場および変位速度場の時間微分を算出する。
・エネルギ汎関数、変位場の勾配、および、変位速度場の勾配を用いて、損失関数を最小化するようにＦＥＭＬＮのパラメータを学習する。

次に、ＦＥＭＬＮの構成例について説明する。図５は、ＦＥＭＬＮの構成例Ｎ１を示す図である。構成例Ｎ１は、エネルギ汎関数（蓄積エネルギ、損失エネルギ、仕事量）と、次の時間ステップにおける変位場および変位速度場との両方が出力データに含まれる構成の例である。図６は、ＦＥＭＬＮの構成例Ｎ２を示す図である。構成例Ｎ２は、エネルギ汎関数のみが出力データに含まれる構成の例である。

出力データがエネルギ汎関数のみの場合は、次の時間ステップの変位場および変位速度場は、ＦＥＭＬＮの出力のエネルギ汎関数から、変分原理により推定することで算出される。

λ_Ｇは、材料特性および構造変数についての条件データを示す。λ_Ｆは、負荷条件および境界条件についての条件データを示す。λ_ｔは、時間に関する条件データを表す。これらの条件データのうち一部が入力されるように構成してもよい。

図５に示すように、構成例Ｎ１のＦＥＭＬＮは、入力層と、２つの中間層と、出力層と、を含む。なお中間数の個数は２個に限られず、１個または３個以上であってもよい。

入力層は、要素ごとおよび節点ごとの変位場および変位速度場、並びに、条件データλが、入力データとして入力される。出力層は、要素ごとのエネルギ汎関数と、要素ごとおよび節点ごとの次の時間ステップにおける変位場および変位速度場と、を出力データとして出力する。

次に、ＦＥＭＬＮの学習に用いられる損失関数について説明する。損失関数は、例えば以下のように定義される。
・ＦＥＭＬＮから算出した勾配と、事前のＦＥＭ解析結果（離散化要素または離散化節点ごとに算出）である正解データの勾配と、の差を最小にできる関数

図７は、損失関数で用いられる各勾配の差を説明するための図である。図７に示すように、勾配は、エネルギ汎関数のスカラー値φ_θの変位場についての偏微分、エネルギ汎関数のスカラー値φ_θの変位速度場についての偏微分、変位速度場の時間についての偏微分、および、変位場の時間についての偏微分を含む。これらの勾配それぞれについて、ＦＥＭＬＮから算出された勾配の値（勾配値）と、正解データの勾配値との差が算出される。

損失関数は、例えば、解析領域における各要素または各節点の勾配値の差の二乗和の総和などである。正解データの勾配値は、例えば、離散化要素または離散化節点ごとに算出した勾配値の集合データである。正解データの勾配値は、以下のような関係モデルを用いて算出されてもよい。
・変位場を変数としたエネルギ汎関数に関する近似モデルなどの関係モデル
・時間を変数とした変位場に関する近似モデルなどの関係モデル

損失関数には、エネルギ汎関数に関する制約を加えてもよい。例えば、対象システムがなす仕事量が蓄積エネルギと損失エネルギとの和に等しくなるという制約条件に基づく損失関数である。

また、損失関数は、エネルギ汎関数についてもＦＥＭＬＮから算出したエネルギ汎関数値と、事前のＦＥＭ解析結果から取得したエネルギ汎関数値との差を最小にできる関数を含んでもよい。

また、ＦＥＭＬＮの出力データが変位および変位速度などの物理量を含む構成（構成例Ｎ１）については、変位および変位速度などの物理量についてもＦＥＭＬＮから算出した値と事前のＦＥＭ解析結果から取得した値との差を最小にできる関数を損失関数として追加してもよい。これらの各損失関数の和における各項の重み係数を変えてもよい。

また、損失関数におけるエネルギ汎関数の変位場または変位速度場に関する勾配について、以下の（４）式および（５）式で表される、条件データλに関する偏微分の連鎖律を適用してもよい。

エネルギ汎関数の条件データλに関する偏微分データと、条件データλの変位場および変位速度場に関する偏微分データは、事前に学習データ（正解データ）として準備される。ＦＥＭＬＮの学習の際には、エネルギ汎関数の条件データλに関する偏微分データ（以下の（６）式）、条件データλの変位場に関する偏微分データ（以下の（７）式）、および、条件データλの変位速度場に関する偏微分データ（以下の（８）式）も算出される。

また、損失関数は、エネルギ汎関数の条件データλに関する偏微分データと、条件データλの変位場および変位速度場に関する偏微分データについて、ＦＥＭＬＮから推定される勾配と、事前の学習データの整合性に関する損失関数を含んでもよい。

次に、本実施形態にかかる制御装置１００によるＦＥＭＬＮの学習処理の流れについて説明する。図８は、本実施形態における学習処理の一例を示すフローチャートである。

取得部１０１は、学習に用いる学習データを取得する（ステップＳ２０１）。例えば取得部１０１は、上記（Ｓ２）で説明したような手順で準備された学習データを取得する。

算出部１１１は、取得された学習データに含まれる入力データをＦＥＭＬＮに入力し、ＦＥＭＬＮが出力する出力データを得る（ステップＳ２０２）。例えば出力データは、エネルギ汎関数の値、次の時間ステップの変位場および変位速度場を含む。算出部１１１は、エネルギ汎関数の値、変位場および変位速度場の勾配をそれぞれ算出する（ステップＳ２０３）。

学習部１１２は、算出された勾配と、学習データに含まれる正解データの勾配との差に基づく損失関数を最小化するように、ＦＥＭＬＮを学習する（ステップＳ２０４）。

学習部１１２は、学習が終了したか否かを判定する（ステップＳ２０５）。例えば学習部１１２は、勾配の差が閾値より小さくなったか、および、学習の回数が上限値に達したか否か、などにより、学習の終了を判定する。

学習が終了していない場合（ステップＳ２０５：Ｎｏ）、ステップＳ２０２に戻り、新たな学習データに対して処理が繰り返される。学習が終了したと判定された場合（ステップＳ２０５：Ｙｅｓ）、学習処理を終了する。

以上のように、本実施形態では、離散化の各要素から構成されたエネルギ汎関数のラグラジアンを活用した変分原理に基づく考え方で、ニューラルネットワークの出力層にエネルギ汎関数が組み込まれる。そして、ラグラジアンの変位場および変位速度場の勾配が、変位場および変位速度場の時間的変化と整合（事前に物理現象シミュレーションした数値実験の結果集合との整合など）するように、ＦＥＭＬＮが学習される。このようにして構築されたＦＥＭＬＮにより、時間依存性の物理現象にも活用できる超高速シミュレーション手法を実現できる。このため、後述するように、ＦＥＭＬＮにより得られる物理量を用いてシステム制御を行うことにより、モニタリングデータなどを用いたリアルタイムのシステム制御を行うことが可能となる。

ＦＥＭＬＮに基づく超高速シミュレーションを行う際には、ＦＥＭＬＮにおける出力の蓄積エネルギと損失エネルギの和が仕事量に等しくなるという制約を設けて、入力データの組み合わせを事前に選定してもよい。または、超高速シミュレーションの入出力データの組み合わせは事前に設定されてもよい。

図９および図１０は、本実施形態の適用例を示す図である。図９は、解析の対象とする領域の全体について、本実施形態を適用する例を示す。対象とする系（対象システム）の解析モデルは、物理現象シミュレーション、または、最適化シミュレーションのためのモデルである。現象解析モデルとして、本実施形態のＦＥＭＬＮなどの推定モデルが適用される。

図１０は、解析の対象とする領域（全体領域）のうち、ＦＥＭなどのシミュレーションにおいて多大な時間を要する一部の領域のみに本実施形態を適用する例を示す。一部の領域に適用するモデル（部分モデル）として、本実施形態のＦＥＭＬＮなどの推定モデルが適用される。この場合、制御装置１００は、全体領域のシミュレーションと、条件データ（境界変位、負荷条件など）をやり取りするインターフェース（データ取得部）を備えればよい。

（変形例）
これまでは連続体力学問題の解析に適用した例を説明したが、適用可能な解析処理はこれに限られない。例えば、偏微分方程式の数理モデルで記述される、以下のような物理現象の解析に本実施形態の手法を適用できる。
・電磁場解析
・構造および磁場の連成解析
・Ｇｉｎｚｂｕｒｇ－Ｌａｎｄａｕ現象問題（超伝導現象などの相転移現象解析）
・デバイスシミュレーション（半導体の電子密度およびホール密度挙動解析など）

以下、各解析で用いられるエネルギ汎関数について説明する。
（Ａ１）電磁場解析で使用するエネルギ汎関数：以下の（９）式

変位誘起誘導電流による外部仕事：Ｗ_ｗ＝Ｊ×（ｕ_ｔ×Ｂ） ［Ｗ］
発熱エネルギ：Ｗ_ｊ＝Ｊ^２／σ ［Ｗ］
ポインティングベクトル：Ｓ＝Ｅ×Ｈ［Ｗ／ｍ^２］
磁束密度：Ｂ［Ｔ］
磁場：Ｈ［Ａ／ｍ］
電束密度：Ｄ［Ｃｍ^－２］
電場：Ｅ［Ｎ／Ｃ］
時間:ｔ［ｓ］
体積:Ω［ｍ^３］
表面積:ｄΩ［ｍ^２］
外向き単位法線ベクトル：ｎ
（Ａ２）構造および磁場の連成解析で使用するエネルギ汎関数：以下の（１０）式

弾性ひずみエネルギ：Ｅ_ｋ＝σ_ｉｊ×ε_ｉｊ／２ ［Ｊ］
運動エネルギ：Ｅ_ｅ＝ρｕ_ｉ，ｔ×ｕ_ｉ，ｔ／２ ［Ｊ］
構造減衰による散逸エネルギ：
Ｗ_ｃ＝α×ρｕ_ｉ，ｔ×ｕ_ｉ，ｔ＋β×σ_ｉｊ，ｔ×ε_ｉｊ ［Ｗ］
コイルの渦電流などによるによる外部仕事：Ｗ_Ｆ＝Ｆ_ｉ×ｕ_ｉ，ｔ ［Ｗ］
変位：ｕ［ｍ］
応力：σ［Ｐａ］
密度：ρ
ひずみ：ε
レイリー減衰の係数：α［ｓ^－１］
レイリー減衰の係数：β［ｓ］
ｉ方向の力：Ｆ_ｉ
（Ａ３）Ｇｉｎｚｂｕｒｇ－Ｌａｎｄａｕ方程式：以下の（１１）式（（１１）式内の「ｉ」は虚数を表す）、（１２）式

Ｇｉｎｚｂｕｒｇ－Ｌａｎｄａｕ方程式で使用するエネルギ汎関数：以下の（１３）式

超伝導エネルギ：以下の（１４）式

磁場によるエネルギ：以下の（１５）式

相互作用エネルギ：以下の（１６）式

外部磁場：Ｂａ
磁場のベクトルポテンシャル：Ａ
Ｇｉｎｚｂｕｒｇ－Ｌａｎｄａｕ パラメータ：κ
超伝導状態を表す秩序変数；ψ
常伝導状態の電気伝導率：σ
（Ａ４）デバイスシミュレーション：
半導体のデバイスシミュレーションにおける結晶欠陥挙動解析で用いられるエネルギ汎関数のＨｅｌｍｈｏｌｔｚ自由エネルギＦは以下の（１７）式のように表現できる。

ケミカルポテンシャル：ｆ_ｃｈｅｍ
欠陥エネルギ：ｆ_ｓｓｆ
弾性ひずみエネルギ：ｆ_{ｅｌａｓｔ}
勾配エネルギ：ｆ_ｇｒａｄ
結晶学的エネルギ：ｆ_{ｃｒｙｓｔ}
外部応力からの仕事量：Ｗ

位置ベクトルｘに関して対象とする領域Ωにおいて、それぞれの離散化要素のエネルギ汎関数ｆの総和（積分）をとることにより、全体のエネルギ汎関数を求めることができる。

半導体のデバイスシミュレーションにおける電子密度およびホール密度分布の挙動解析は時間を要するため、結晶欠陥挙動解析におけるケミカルポテンシャル計算のための電子密度およびホール密度分布の挙動解析のみにＦＥＭＬＮを適用することも可能である（上記図１０の例）。

ケミカポテンシャルｆ_ｃｈｅｍは、以下の（１８）式、（１９）式により算出される。

ケミカルポテンシャル関数：μ
結晶欠陥有無の状態量：φ

電子密度とホール密度のケミカルポテンシャル関数は以下の（２０）式および（２１）式で算出できる。

電場Ｅにおける電子密度：Ｄ_ｅ（Ｅ）
電場Ｅにおけるホール密度：Ｄ_ｈ（Ｅ）
電子のＦｅｒｍｉ－Ｄｉｒａｃ分布関数：Ｆ_ｎ（Ｅ）
ホールのＦｅｒｍｉ－Ｄｉｒａｃ分布関数：Ｆ_ｐ（Ｅ）
電子密度：ｎ
ホール密度：ｐ
ボルツマン定数：ｋ_Ｂ
温度：Ｔ
伝導帯：Ｅ_ｃ
価電子帯：Ｅ_ｖ
真性半導体バンド：Ｅ_ｉ

電子密度分布およびホール密度分布は、ボルツマン方程式、ポアソン方程式、および、電流連続方程式をセルフコンシステントに解くことにより求めることができる。

次に、制御部１０２による制御処理の詳細について説明する。図１１は、制御部１０２による制御処理の例を説明するための図である。なお図１１は、フィードフォワード制御を行う制御部１０２ａと、フィードバック制御を行う制御部１０２ｂと、の２つの制御部、および、対応する各部（記憶部１２１ａ、１２１ｂ、演算処理部１１０ａ、１１０ｂ）を含む構成の例を示す。制御部１０２ａ、１０２ｂのいずれか一方および対応する各部のみを含むように制御装置１００を構成してもよい。また、３つ以上の制御部１０２および対応する各部を含むように制御装置１００を構成してもよい。

演算処理部１１０ａ、１１０ｂは、それぞれ制御部１０２ａ、１０２ｂによる制御処理で用いられる物理量などのデータの演算を行う。

記憶部１２１ａ、１２１ｂは、それぞれ演算処理部１１０ａ、１１０ｂによる処理で用いられるデータを記憶する。例えば記憶部１２１ａ、１２１ｂは、演算処理部１１０ａ、１１０ｂがそれぞれ用いる推定モデル（ＦＥＭＬＮなど）を定める情報、および、対応する推定モデルによる処理で用いられるデータ（入力データ、演算結果、学習データなど）を記憶する。

制御対象２００ａは、例えば、外乱要因からセンサまでの応答特性に相当する。制御対象２００ｂは、例えば、対象となる物体を動作させる制御アクション部からセンサまでの応答特性に相当する。

演算処理部１１０ａの推定モデルには、外乱の影響を受けた変位場および変位速度場を表す入力データが入力される。制御部１０２ａは、演算処理部１１０ａにより算出される時刻ｔｂの物理量Ｐｈｂから、フィードフォワード制御のための制御量を求め、求めた制御量によりフィードフォワード制御を行う。

制御部１０２ｂは、演算処理部１１０ｂにより算出された変位場および変位速度場と、計測された変位場および変位速度場と、の差が小さくなるようにフィードバック制御を行う。言い換えると、制御部１０２ｂは、制御量の一例である差が目標値の一例であるゼロとなるようにフィードバック制御を行う。

制御対象はどのような対象であってもよいが、以下では、例えば発電システムに用いられる物体であって、自励振動または自励回転を誘起する物体（振動物体または回転物体）を制御対象とする例を主に説明する。自励回転は、自励振動を回転に変換するクランク機構（カウンターウェイトを設置してもよい）などの機構により実現することができる。以下では自励振動を行う振動物体への適用例を主に説明するが、「振動」を「回転」に置き換えれば、自励回転を行う回転物体に対しても同様の手順を適用できる。自励振動（自励回転）を行う振動物体（回転物体）は、自励振動子（自励回転子）とも呼ばれる。

例えば制御部１０２ｂは、安定かつロバストな自励振動（リミットサイクル）になるように振動物体を制御する。より具体的には、制御部１０２ｂは、振動物体を振動させるための励起力と、振動物体自体の変位速度（または迎角α）との理想的な関係を目標値として、振動物体を制御する。励起力は、無次元化、正規化または標準化した値である励起力係数で表されてもよい。以下では、励起力を励起力係数で表す例を主に説明する。制御部１０２ｂは、振動物体の周辺（例えば風上など）の媒体（以下、周辺媒体。例えば風および水など）の流速が、安定的なリミットサイクルを励起する臨界流速Ｕ_ｃを超える条件などの制約条件を設定してもよい。

振動物体に対して用いる推定モデル（例えばＦＥＭＬＮ）の入力データは、例えば以下のようなデータを含む。
・周辺媒体の流速（１つでも複数でもよい）などの流れ場（渦度など他の流れ場を含めてもよい）
・温度場および湿度場（周辺媒体の密度などの流れ場に関する物性値の推定に活用可能）
・振動物体の表面の圧力分布（各節点または各要素）などの変位場
・振動物体の表面の圧力変化速度分布（各節点または各要素）などの変位速度場
・振動物体自体の振動変位（例えば振動物体の重心位置の変位）
・振動物体自体の振動変位速度（例えば振動物体の重心位置の変位速度）
・時間に関する条件データλ_ｔ

さらに入力データは、センシングデータを活用したデータ同化またはキャリブレーションのために、振動物体の材料モデルに関する変数（例えば部材の弾性率など）、振動物体の構造モデルに関する変数（例えば板厚など）を含んでもよい。

例えば上記の周辺媒体の流速などの流れ場は、図５および図６の負荷条件および境界条件についての条件データλ_Ｆに相当する。例えば材料モデルに関する変数および構造モデルに関する変数は、図５および図６の材料特性および構造変数についての条件データλ_Ｇに相当する。

流速、変位、および、圧力などの物理量は、ベクトルの各成分でもよいし、縮約した値またはノルムでもよい。また、上記物理量は、ある時刻ｔの数値でもよいし、ある時刻ｔより前の時刻ｔ－ｎΔｔまでの平均または移動平均をとった数値であってもよい。

推定モデル（例えばＦＥＭＬＮ）により、時刻ｔ（時刻ｔａ）での変位場、変位速度場、負荷・境界条件、および、材料・構造変数の入力のもと、出力であるエネルギ汎関数を通じて、時刻ｔから次の時間ステップ後の時刻ｔ＋Δｔ（時刻ｔｂ）での変位場および変位速度場が推定される。また、この処理を繰り返すことにより、変位場および変位速度場の時間的変化（時刻歴）が推定される。これにより、励起力係数と振動物体の変位速度との関係、または、励起力係数と周辺媒体の流れ場（例えば風速）との関係が推定できる。

制御部１０２は、これらの関係の目標値に整合するように、振動物体の表面形状、振動物体の向き、および、気流制御機構などの制御量（制御変数、制御パラメータ）を制御する。目標値は、例えばデータテーブル、関係式、および、数理モデルなどにより事前に用意される。

計測された周辺媒体の流速、温度、および、湿度データなどにより、周辺媒体の流れ場に関する物理量が変化する。このため、制御部１０２は、事前に用意している目標値を、センシングデータに合わせて、逐次、修正してもよい。

図１２は、自励振動する振動物体３０２の例を示す図である。振動物体３０２は、翼形状の構造（翼構造）を有し、例えば流速Ｕの気体（気流）に応じて振動する。この場合、振動物体３０２は例えば風力発電に用いることができる。振動物体３０２の形状は翼型に限らないが、気流の抵抗を低減できる形状が望ましい。なお気体（気流）の代わりに液体（水流）、または、気体と液体の混合物（混相流）などの流体が用いられてもよい。すなわち、振動物体３０２は、潮流発電および海流発電などに用いられてもよい。以下では、主に気体を流体として用いる場合を例に説明する。

振動物体３０２は、粘性パラメータｒ、および、弾性パラメータｋなどにより表される。振動物体３０２は、振動により図１２の上下方向に変位する。以下では、振動の方向をｙ方向とする座標を用いるものとする。図１２では、変位をｙとしたときの速度（ｙドット）が表されている。制御部１０２は、振動物体３０２を制御対象とし、自励振動を誘起するように振動物体３０２の翼の形状、翼の向き、粘性パラメータｒ、および、弾性パラメータｋなどを制御する。これらの制御量は、制御入力データとして制御部１０２に入力される。

また図１２では、振動物体３０２の前方に、気流制御機構３０１が備えられる例が示されている。気流制御機構３０１は、振動物体３０２の周辺の気流を制御する機構（流れ場の増速および整流など）である。図１２の例では、気流制御機構３０１は、流速ｕの気流を、流速Ｕ（ｕ＜Ｕ）の気流に増速する。気流制御機構３０１は、振動物体３０２の前方に備えられる必要はなく、振動物体３０２の後方を含む、振動物体３０２の周辺のいずれの方向に備えられてもよい。

気流制御機構３０１は、例えば、複数の翼構造を含み、翼の形状および翼の向きを制御可能とされてもよい。すなわち制御部１０２は、振動物体３０２が振動を誘起するように、気流制御機構３０１の翼の形状等を制御してもよい。

図１３は、自励振動する振動物体の他の例を示す図である。図１３の振動物体３０４は、回転プロペラ翼形状の振動物体の例である。部分翼３０５は、回転プロペラ翼形状の１枚の翼の一部を切り出した部分に相当する。振動物体３０４の周辺の気流を制御する気流制御機構３０３が備えられてもよい。

振動物体３０２、３０４の変位（変位分布）または変形（変形分布）は、どのような方法で計測されてもよい。例えば、撮像装置（カメラなど）により撮影される画像を用いる方法、および、振動物体３０２の内部に備えられた複数のばねの変位を計測することにより振動物体３０２の表面の変位・変形分布を計測する方法を適用できる。図１４は、画像を用いて振動物体３０２の表面の変位を計測する計測方法の例を説明するための図である。なお部分翼３０５の表面の変位についても同様の計測方法を適用できる。

センサの１つとして機能するカメラ１４０１、１４０２は、振動物体３０２の表面を撮影し、時系列の画像を出力する。拡大部１４１１に示すように、振動物体３０２の表面には、撮影された画像で識別可能な模様が付与される。画像は例えば取得部１０１により取得され、演算処理部１１０に出力される。演算処理部１１０は、取得された画像に対して、例えばデジタル画像相関法などを適用することにより、振動物体３０２の変位、変形、または、ひずみの時間的分布および空間的分布を算出する。

次に、翼構造を有する振動物体（振動物体３０２、部分翼３０５など）の制御方法の例について説明する。図１５は、ＮＡＣＡ（National Advisory Committee for Aeronautics）翼型の振動物体３１１－１、３１１－２、３１１－３を制御対象とする制御の例を示す図である。振動物体３１１－１、３１１－２、３１１－３は、制御部１０２により、相互に異なる形状となるように制御された振動物体に相当する。

制御部１０２は、翼の形状および翼の傾きの少なくとも一方を変更するように、例えば形状を変更させるアクチュエータ、および、傾きを変更させるアクチュエータの動作を制御する。図１５の５つの矢印は、翼の形状が５自由度で制御可能であることを示す。自由度は５に限られず、２～４または６以上であってもよい。制御部１０２は、振動物体３１１－１～３１１－３の表面がなめらかになるなどの制約を満たすように翼の形状等を変更するように制御してもよい。

図１６は、翼構造を有する振動物体の制御方法の他の例を示す図である。図１６は、振動物体の表面に備えられる複数の羽の傾きを変更することにより、振動物体の翼の形状を制御する例を示す。図１６では、振動物体が５枚の羽３２１－１～３２５－１を備える例が示されているが、羽の枚数は５枚（５自由度）に限られない。羽３２１－２～３２５－２、および、羽３２１－３～３２５－３は、それぞれ羽３２１－１～３２５－１とは異なる傾きに変更された羽の組の例である。

制御部１０２は、複数の羽の傾き、および、翼全体の傾きの少なくとも一方を変更するように、例えば羽の傾きを変更させるアクチュエータ、および、翼全体の傾きを変更させるアクチュエータの動作を制御する。

図１７は、翼構造を有する振動物体の制御方法の他の例を示す図である。図１７は、図１６の複数の羽（３２１－１～３２５－１など）の上部を覆う弾性材料３３１－１（３３１－２、３３１－３）をさらに備える振動物体の制御の例である。

この例の場合も、制御部１０２は、図１６の同様の方法により、複数の羽の傾き、および、翼全体の傾きの少なくとも一方を変更することができる。

図１５～図１７の制御方法は一例であり、これらに限られるものではない。例えば制御部１０２は、アクチュエータにより翼形状の振動物体の表面に誘起流を発生させることで、振動物体に自励振動を誘起させるように制御してもよい。誘起流は、例えば、プラズマ誘起流、電気浸透流、温度・濃度勾配流、および、振動物体の内部から表面の小さな孔を介して出力される噴流などである。

これまでは１つの振動物体の制御について説明したが、制御装置１００は、複数の振動物体を制御するように構成されてもよい。このような構成の場合、制御部１０２は、さらに複数の振動物体の相互作用により、位相引き込み現象を誘起するように制御を行ってもよい。相互作用は、複数の振動物体を共通に支持する支持台からの振動伝播による相互作用、および、気流による各自励振動の相互作用などを含む。

複数の振動物体から発生する励起力を、例えばクランク機構により回転運動に変換する構成の場合、クランクシャフトが、カウンターウェイトによるバランサーを備えてもよい。また、回転を滑らかにするために、一般的なクロスプレーン形式のように、角度位相をずらして各連結部が設置されてもよい。さらに、複数の振動物体から安定的で滑らかな回転力を得るために、クラッチ機構が設けられてもよい。

図１８～図２０は、複数の振動物体の相互作用の例を示す図である。図１８は、支持台１８０１の振動伝播により、複数の振動物体３０２－１～３０２－４に相互作用が生じる例を示す。図１９は、流れ場（周辺の気流）により複数の振動物体３０２－１～３０２－３に相互作用が生じる例を示す。

図２０は、周辺の場からの作用を介して複数の振動物体３０２－１～３０２－４に相互作用が生じる例を示す。周辺の場からの作用は、例えば、振動物体３０２－１～３０２－４の周辺の流れ場の増速および整流を行う制御機構２０１０により制御される流れ場の作用である。制御機構２０１０は、例えば、増速および整流のための翼型の表面形状の制御機構２０１１を備える。

複数の振動物体が備えられる場合は、位相引き込み現象を誘起するように制御することにより、より安定的な励起力を発生し、外乱に対してロバストな制御を行うことが可能となる。この場合、目標値は、励起力係数と振動物体の動作に影響する物理量との関係が、安定かつロバストな自励振動を誘起し、かつ、複数の振動物体の相互作用により引き込み現象を誘起する関係となるように定められる。図２１は、引き込み現象の誘起例を示す図である。図２１の横軸は時間を表し、縦軸は励起力係数を表す。

グラフ２１０１は、複数の振動物体間に相互作用が生じていない場合の励起力の変化を表す。グラフ２１０２は、複数の振動物体間に相互作用が生じている場合の励起力の変化を表す。複数の振動物体間の相互作用の大きさを定める変数（相互作用係数など）を適切な値に制御することにより、複数の振動物体の振動の位相を同期させ、大きな励起力を安定的に発生させることが可能となる。すなわち、構造安定な自励振動（リミットサイクル）の相互引き込み現象を発生させることができる。

次に、励起力により自励振動する振動物体（以下、自励振動物体という場合がある）の力学モデルについて説明する。

（２２）式は、励起力を表す式である。Ｃは励起力係数、ρは周辺媒体（空気など）の密度、ａは振動物体正面の断面の等価面積、Ｕは流速、αは迎角を表す。

図２２は、振動物体の振動の速度（ｙドット）または迎角αに対する励起力係数Ｃの変化を表すグラフの例を示す図である。なお振動の速度と迎角αとは、図２２の上部に示すような関係を有する。Ｕ_αは迎角α方向の流速に相当する。なお、以下では主に振動の速度と励起力係数との関係を用いて説明するが、振動の速度を迎角αに置き換えても同様の手順を適用できる。図２２の下部は、翼の形状および翼の向きなどが相互に異なる２つの振動物体に対するグラフを表す。これらのグラフが示すように、例えば振動物体の翼の形状および翼の向きなどの制御に応じて、励起力係数の変化を表すグラフは変動しうる。

励起力係数と迎角との関係は、例えば７次関数で近似すると以下の（２３）式のように表される。ｍは振動物体の質量を表す。

（２３）式を３次関数までで打ち切ると、変数変換により、励起力係数と迎角との関係はＶａｎ ｄｅｒ Ｐｏｌ方程式に帰着できる。（２４）式は、励起力係数と迎角との関係を３次関数で表した式の例である。

（２３）式は、励起力に相当する励起項（右辺第１項）が、リミットサイクルを生じさせるような変位速度に関する１次関数と３次関数とを含む式となっている。励起項は、発散性を伴う振動（１次関数）を非線形項（３次関数）が抑制する項であると解釈することができる。なお励起項は、５次関数、および、７次関数などの３次以上の関数をさらに含んでもよい。

周辺媒体の流速Ｕ（例えば風速）が臨界流速Ｕ_ｃを超えると、負の空力減衰が正の構造体減衰より大きくなり、安定的なＨｏｐｆ分岐（動的分岐）にリミットサイクルが生成される。流速Ｕがさらに大きくなると２つの安定リミットサイクルの間に１つの不安定リミットサイクルがある状態が現れる。

図２３は、振動物体の変位（横軸）と、変位速度（縦軸）との関係の例を示す図である。図２３の線２３０１は、初期状態に依存しない構造安定なリミットサイクルを表す。制御部１０２は、安定なリミットサイクルを目標として振動物体の振動を制御する。なお、リミットサイクルが生じる条件は以下の（２５）式で表される。

（２６）式および（２７）式は、複数（ｎ個、ｎは２以上の整数）の振動物体の引き込み現象を含む力学モデルの例を示す。ｙ_ｉ（１≦ｉ≦ｎ）は、ｉ番目の振動物体の変位を表す。

（２６）式は、１つの相互作用係数εにより相互作用項（右辺第２項）が表される例である。相互作用係数εは、例えば振動物体を支持する支持台の振動伝播係数である。（２６）式の右辺第３項は、外乱に相当する揺動項を表す。（２６）式は、例えば図１８に示す相互作用のモデルに相当する。

（２７）式は、振動物体ごとに定められる相互作用係数ε_１～ε_ｎが用いられる例を示す。（２７）式は、例えば図１９、図２０に示す相互作用のモデルに相当する。

なお（２６）式および（２７）式は、複数の振動物体における全体の挙動をＧｉｎｚｂｕｒｇ－Ｌａｎｄａｕ方程式に帰着できるようにモデル化した式である、と解釈することができる。

演算処理部１１０（算出部１１１）は、周辺の気体の流速、流量、および、流れの方向（温度、湿度を含めてよい）などを条件データとして入力して、ＦＥＭＬＮを用いて、振動物体の制御に用いる物理量を算出する。物理量は、例えば、振動物体の変位場および変位速度場である。

以下、物理量の算出、および、算出した物理量による制御の例について説明する。

（制御例１）
制御例１は、振動物体の表面の変位分布をセンサなどによりセンシングできない場合、例えば、上記の撮像装置により撮影される画像を用いた変位分布の計測方法を適用できない場合の例である。

算出部１１１は、ＦＥＭＬＮとは異なる推定モデルなどを用いることにより、変位場および変位速度場から振動物体の表面の圧力分布を算出する。ＦＥＭＬＮが、物理量として振動物体の表面の圧力分布を算出するように構成してもよい。

算出部１１１は、さらに、算出された圧力分布を振動物体の表面に渡って積分することで、振動物体に負荷される励起力係数の時間的変化および振動物体の変位速度に対する依存性（変位速度依存性）を算出する。

制御部１０２は、例えば以下の領域の圧力が向上するように、振動物体の表面形状および向きの少なくとも一方を変化させるように制御する。
・振動物体の表面の圧力が低下した領域
・圧力分布が乱れる領域
・圧力分布が細動する領域
・上記の各領域の周辺領域

例えば、圧力が向上する表面形状および向きと、圧力分布との関係を記述したデータベースまたはモデルベースを記憶部１２１に記憶しておく。制御部１０２は、このデータベースまたはモデルベースを用いて、圧力が向上するような振動物体の形状の変化量を求める。

制御部１０２は、目標とする励起力係数と振動物体の挙動（制御対象の動作に影響する物理量）との間の関係を示す以下のような制御量を事前に設定しておき、目標との誤差が小さくなるように、振動物体の形状を制御してもよい。制御量は、例えば、データテーブル、関係式、または、数理モデル（例えば、３次関数の係数Ａ_１、Ａ_２など）で表される。
・目標とする励起力係数（または振動物体の表面の圧力分布）と風速との関係
・目標とする励起力係数（または振動物体の表面の圧力分布）と振動物体の振動速度との関係
・目標とする励起力係数（または振動物体の表面の圧力分布）と迎角との関係

（制御例２）
制御例２は、振動物体の表面の変位分布を撮像装置（センサの一例）によりセンシングできる場合の例である。

演算処理部１１０（算出部１１１）は、撮像装置などにより取得された振動物体の表面の画像履歴を、デジタル画像相関法により変位分布の履歴に変換する。算出部１１１は、得られた変位分布を、例えばＦＥＭＬＮなどの推定モデルを用いて振動物体の表面全体の圧力分布（圧力分布の時間的変化）を算出する。

算出部１１１は、さらに、算出された圧力分布を振動物体の表面に渡って積分することで、振動物体に負荷される励起力係数の時間的変化および変位速度依存性を算出する。その他の処理は、制御例１と同様である。

（制御例３）
制御例３は、振動物体の表面全体の圧力分布を圧力センサによりセンシングできる場合の例である。

演算処理部１１０（算出部１１１）は、圧力センサによりセンシングされる振動物体の表面の圧力分布の時間的変化を取得する。算出部１１１は、取得した圧力分布を振動物体の表面に渡って積分することで、振動物体に負荷される励起力係数の時間的変化を算出する。その他の処理は、制御例１と同様である。

なお制御例３は、ＦＥＭＬＮを用いずに制御する例と解釈することができる。このように、より高速に物理量を得られる他の方法を適用できる場合は、ＦＥＭＬＮなどの推定モデルを用いなくても、モニタリングデータなどを用いたリアルタイムのシステム制御を実現可能である。

（制御例４）
制御例４は、振動物体の表面一部の変位分布を、変位計測計または撮像装置（画像）によりセンシングできる場合の例である。

演算処理部１１０（算出部１１１）は、例えばＦＥＭＬＮにより振動物体の表面の圧力分布を推定する。算出部１１１は、センシングできる領域のノイズを含む圧力分布と、推定した圧力分布とから、データ同化により振動物体の表面の圧力分布を推定する。算出部１１１は、センシングした変位分布と整合（誤差最小化など）するようにＦＥＭＬＮモデルを修正してもよい。算出部１１１は、得られた圧力分布を振動物体の表面に渡って積分することで、振動物体に負荷される励起力係数の時間的変化を算出する。その他の処理は、制御例１と同様である。

（制御例５）
制御例５は、複数の振動物体がある場合の例である。制御部１０２は、上記制御例１～４のような制御に加え、複数の振動物体の振動位相の引き込み現象が誘起されるように、複数の振動物体間の相互作用係数εを制御する。

制御部１０２は、以下のようなパラメータを制御する制御量に含めてもよい。
・弾性パラメータｋ
・粘性パラメータｒ
・弾性パラメータｋと振動物体の質量ｍとの比ｋ／ｍ
・粘性パラメータｒと振動物体の質量ｍとの比ｒ／ｍ

また、各振動物体に関するｍ、ｋ、ｒは、一部または全部の振動物体の間で異なっていてもよい。ここで、粘性パラメータは、発電機構、ギアおよび軸受などの、摩擦および非弾性変形による機械的なエネルギ損失、並びに、電気的なダンピングなどのエネルギ損失も含んでいる。

（制御例６）
上記の各制御例は、振動物体の自励振動の励起力を制御する方法である。これらの制御方法に加えて、自励振動物体の現在位置と速度とをフィードバックする可変速度フィードバックを行ってもよい。図２４は、可変速度フィードバックを行う機能の構成例を示す図である。

可変速度フィードバックを行う制御部１０２は、例えば、外乱（励起力）を加振源とする系において、自励振動物体の速度に比例した制御力を与える速度フィードバックループ、振幅に応じて速度フィードバックゲインＫを変化させる振幅コントローラ２４０１、および、振幅検出器２４０２を含む。

可変速度フィードバックの具体的な実現方法は、どのような方法であってもよいが、例えば以下の方法を適用できる。

（２８）式は、強制外力がない場合の自励振動物体の運動方程式を表す。

以下のように、粘性パラメータｒと速度フィードバックゲインＫとの差の正負で自励振動物体の安定性すなわち、（２８）式の解の安定性が決まる。
ｒ－Ｋ＜０：自励発振（不安定）
ｒ－Ｋ＝０：持続振動（安定限界）
ｒ－Ｋ＞０：アクティブ制振（安定）

速度フィードバックゲインＫが粘性パラメータｒより大きければ（ｒ－Ｋ＜０）、振動物体は自励振動を発生する。振幅が大きくなると、制御部１０２は、速度フィードバックゲインＫの大きさを振動の大きさに応じて変化させ、振幅を一定に保たせる。制御部１０２は、振動が小さい場合には、速度フィードバックゲインＫを大きくし、自励振動を発生させてもよい。制御部１０２は、振動が大きいときには、速度フィードバックゲインＫを小さくしてアクティブ制振を行ってもよい。

速度フィードバックゲインＫの値を調整する振幅コントローラ２４０１としては、振幅偏差の比例項、振幅偏差の３乗項、および、振幅偏差を時間に関して積分した項、の３つの項から構成する方法などが挙げられる。

図２５および図２６は、複数の自励振動物体が相互引き込みを起こす構成に対する可変速度フィードバック機能の構成例を示す図である。図２５は、複数の自励振動物体が相互作用する構成（図１８、図１９など）に対する可変速度フィードバック機能の構成例である。図２６は、複数の自励振動物体が周辺の場から相互作用を受ける構成（図２０など）に対する可変速度フィードバック機能の構成例である。

図２５の制御部１０２－１および１０２－２は、例えば複数の振動物体のうち、それぞれ振動物体３０２－１および３０２－２に対する制御を行う。ｙ_１およびＫ_１は、それぞれ振動物体３０２－１の変位、および、振動物体３０２－１に対する速度フィードバックゲインを表す。ｙ_２およびＫ_２は、それぞれ振動物体３０２－２の変位、および、振動物体３０２－１に対する速度フィードバックゲインを表す。Ｋ’は、相互引き込み現象に関する速度フィードバックゲインを表す。上記のように、ｍ、ｋ、ｒは、振動物体ごとに異なる値であってもよい。

制御部１０２－１に含まれる振幅コントローラ２４０１ａおよび振幅検出器２４０２ａは、図２４の振幅コントローラ２４０１および振幅検出器２４０２と同様の機能を備える。制御部１０２－２に含まれる振幅コントローラ２４０１ｂおよび振幅検出器２４０２ｂは、図２４の振幅コントローラ２４０１および振幅検出器２４０２と同様の機能を備える。

なお図２５は２つの振動物体に対する制御部１０２－１および１０２－２の例を示しているが、３つ以上の振動物体に対する制御部１０２も同様に構成することができる。

また図２６は、説明を簡単にするため２つ（ｎ＝２）の振動物体に対する制御部１０２－１および１０２－２の例を示しているが、３つ以上（ｎ≧３）振動物体に対する制御部１０２も同様に構成することができる。

上記の各制御例は、振動物体の形状制御だけではなく、気流制御機構内に備えられる、流れ増速などのための翼型形状の物体の表面形状および向きの制御にも適用できる。

次に、振動物体の振動に基づき発電を行う機構の構成例を説明する。図２７は、電磁誘導型の振動発電（回転発電）の構成例を示す図である。振動発電に用いられる振動物体３０６は、連結部２７０１を介して発電機構２７１０（後述の可動子２７２２）に接続される。

図２８は、発電機構２７１０の構成例を示す図である。発電機構２７１０は、ばね２７１１と、鉄心コア２７１２と、磁石２７１３と、コイル２７１４と、ばね２７１５と、を含む。また発電機構２７１０は、固定子２７２１と、可動子２７２２とに分けることができる。固定子２７２１は、コイル２７１４を含む。可動子２７２２は、ばね２７１１、鉄心コア２７１２、磁石２７１３、および、ばね２７１５を含む。可動子２７２２は、連結部２７０１を介して振動物体３０６に接続され、振動物体３０６の振動に応じて振動する。このとき固定子２７２１は固定状態であり、振動しない。

発電機構２７１０を用いる場合、制御部１０２は、発電機構２７１０の電気的な物理量を制御量として制御してもよい。

また制御部１０２は、振動物体および周辺支持部の変形量、変位、および、加速度を制御量として制御してもよい。

このように、本実施形態にかかる制御装置では、発電に用いる振動物体の制御などに用いることができる物理量をより高速に推定することができる。

例えば、本実施形態では、物体周辺の流れ場は複雑になることが多く、時間的および空間的な流速分布を推定することは困難であることに着目し、このような推定を可能とする推定モデルを実現する。具体的には、本実施形態の制御装置は、周辺の流れ場により物体の表面に変位分布を生じさせる、または、変形させる機構（構造）と、変位に関する物理量（変位場、変位速度場）を計測する機構と、変位場および変位速度場から物体が受ける圧力（時間的分布、空間的分布）を推定する推定モデルと、を備える。これにより、物体が受ける圧力分布から物体周辺の流れ場を推定することができ、推定された情報を元に最適な気流制御を行うことができる。

次に、本実施形態にかかる制御装置のハードウェア構成について図２９を用いて説明する。図２９は、本実施形態にかかる制御装置のハードウェア構成例を示す説明図である。

本実施形態にかかる制御装置は、ＣＰＵ５１などの制御装置と、ＲＯＭ（Read Only Memory）５２やＲＡＭ５３などの記憶装置と、ネットワークに接続して通信を行う通信Ｉ／Ｆ５４と、各部を接続するバス６１を備えている。

本実施形態にかかる制御装置で実行されるプログラムは、ＲＯＭ５２等に予め組み込まれて提供される。

本実施形態にかかる制御装置で実行されるプログラムは、インストール可能な形式又は実行可能な形式のファイルでＣＤ－ＲＯＭ（Compact Disk Read Only Memory）、フレキシブルディスク（ＦＤ）、ＣＤ－Ｒ（Compact Disk Recordable）、ＤＶＤ（Digital Versatile Disk）等のコンピュータで読み取り可能な記録媒体に記録してコンピュータプログラムプロダクトとして提供されるように構成してもよい。

さらに、本実施形態にかかる制御装置で実行されるプログラムを、インターネット等のネットワークに接続されたコンピュータ上に格納し、ネットワーク経由でダウンロードさせることにより提供するように構成してもよい。また、本実施形態にかかる制御装置で実行されるプログラムをインターネット等のネットワーク経由で提供または配布するように構成してもよい。

本実施形態にかかる制御装置で実行されるプログラムは、コンピュータを上述した制御装置の各部として機能させうる。このコンピュータは、ＣＰＵ５１がコンピュータ読取可能な記憶媒体からプログラムを主記憶装置上に読み出して実行することができる。

本発明のいくつかの実施形態を説明したが、これらの実施形態は、例として提示したものであり、発明の範囲を限定することは意図していない。これら新規な実施形態は、その他の様々な形態で実施されることが可能であり、発明の要旨を逸脱しない範囲で、種々の省略、置き換え、変更を行うことができる。これら実施形態やその変形は、発明の範囲や要旨に含まれるとともに、特許請求の範囲に記載された発明とその均等の範囲に含まれる。

１００ 制御装置
１０１ 取得部
１０２ 制御部
１１０ 演算処理部
１１１ 算出部
１１２ 学習部
１２１ 記憶部

Claims (20)

1. 制御対象を示す解析領域を離散化した複数の要素ごとの第１時刻の第１物理量と、前記制御対象の負荷条件、境界条件、材料の条件、および、構造条件の少なくとも一部を表す条件データと、を含む複数の入力データを取得する取得部と、
   複数の前記入力データを推定モデルに入力して得られる、複数の前記要素のエネルギを表すエネルギ汎関数の値を用いた変分原理により、前記第１時刻より後の第２時刻の複数の前記要素ごとの第２物理量を算出する算出部と、
   前記第２物理量と計測された物理量との差に基づく制御量が目標値となるように前記制御対象を制御する制御部と、
   を備える制御装置。
2. 前記算出部は、前記エネルギ汎関数の値の総和が極小になる条件により、前記第２時刻の複数の前記要素ごとの第２物理量を算出する、
   請求項１に記載の制御装置。
3. 前記エネルギ汎関数の値は、
   複数の学習用の入力データと複数の前記学習用の入力データそれぞれに対して数値計算により算出した複数の前記要素のエネルギを表すエネルギ汎関数の値とを含む正解データに基づき学習された前記推定モデルに、複数の前記入力データを入力し、前記推定モデルが出力する値であって、複数の前記要素のエネルギを表す、
   請求項１に記載の制御装置。
4. 前記制御対象は、励起力を発生させる自励振動または自励回転を行う物体であり、
   前記制御部は、前記励起力と前記制御対象の動作に影響する物理量との関係を示す前記制御量が、前記目標値となるように前記制御対象を制御する、
   請求項１に記載の制御装置。
5. 前記目標値は、前記自励振動または前記自励回転を誘起する前記励起力と前記物理量との関係に基づいて定められる、
   請求項４に記載の制御装置。
6. 前記目標値は、前記励起力と前記物理量との関係を表すデータテーブル、関係式、または、数理モデルにより表される、
   請求項５に記載の制御装置。
7. 前記制御対象は、それぞれ前記自励振動または前記自励回転を行う複数の物体であり、
   前記目標値は、前記自励振動または前記自励回転を誘起し、かつ、複数の前記物体の相互作用により引き込み現象を誘起する前記励起力と前記物理量との関係に基づいて定められる、
   請求項４に記載の制御装置。
8. 前記エネルギ汎関数の値の勾配と、勾配の正解データとの差を最小化するように前記推定モデルを学習する学習部をさらに備える、
   請求項１に記載の制御装置。
9. 前記学習部は、さらに、前記エネルギ汎関数の値と、前記エネルギ汎関数の値の正解データとの差を最小化するように前記推定モデルを学習する、
   請求項８に記載の制御装置。
10. 前記推定モデルは、ニューラルネットワークモデル、または、階層ベイズモデルである、
    請求項１に記載の制御装置。
11. 前記推定モデルは、前記第１時刻の複数の前記入力データを入力し、前記第２時刻の前記エネルギ汎関数の値と、前記第２物理量を出力する、
    請求項１に記載の制御装置。
12. 前記エネルギ汎関数は、連続体力学の解析に用いられるエネルギ汎関数であって、複数の前記要素ごとの、蓄積エネルギ、損失エネルギ、および、与えられる仕事量により算出されるエネルギを表す、
    請求項１に記載の制御装置。
13. 前記エネルギ汎関数は、電磁場解析に用いられるエネルギ汎関数であって、複数の前記要素ごとの、発熱エネルギ、および、誘導電流による仕事量により算出されるエネルギを表す、
    請求項１に記載の制御装置。
14. 前記エネルギ汎関数は、構造および磁場の連成解析に用いられるエネルギ汎関数であって、複数の前記要素ごとの、弾性ひずみエネルギ、運動エネルギ、散逸エネルギ、および、渦電流による仕事量により算出されるエネルギを表す、
    請求項１に記載の制御装置。
15. 前記エネルギ汎関数は、相転移現象の解析に用いられるエネルギ汎関数であって、複数の前記要素ごとの、超電導エネルギ、磁場によるエネルギ、および、相互作用エネルギにより算出されるエネルギを表す、
    請求項１に記載の制御装置。
16. 前記エネルギ汎関数は、電子密度およびホール密度の解析に用いられるエネルギ汎関数であって、複数の前記要素ごとの、ケミカルポテンシャル、欠陥エネルギ、弾性ひずみエネルギ、勾配エネルギ、結晶学的エネルギ、および、外部応力からの仕事量により算出されるエネルギを表す、
    請求項１に記載の制御装置。
17. 制御対象を示す解析領域を離散化した複数の要素ごとの第１時刻の第１物理量と、前記制御対象の負荷条件、境界条件、材料の条件、および、構造条件の少なくとも一部を表す条件データと、を含む複数の入力データを取得する取得部と、
    複数の前記入力データを推定モデルに入力して得られる、複数の前記要素のエネルギを表すエネルギ汎関数の値と、前記第１時刻より後の第２時刻の複数の前記要素ごとの第２物理量と、を含む出力データを求める算出部と、
    前記第２物理量と計測された物理量との差に基づく制御量が目標値となるように前記制御対象を制御する制御部と、
    を備える制御装置。
18. 制御装置で実行される制御方法であって、
    制御対象を示す解析領域を離散化した複数の要素ごとの第１時刻の第１物理量と、前記制御対象の負荷条件、境界条件、材料の条件、および、構造条件の少なくとも一部を表す条件データと、を含む複数の入力データを取得する取得ステップと、
    複数の前記入力データを推定モデルに入力して得られる、複数の前記要素のエネルギを表すエネルギ汎関数の値を用いた変分原理により、前記第１時刻より後の第２時刻の複数の前記要素ごとの第２物理量を算出する算出ステップと、
    前記第２物理量と計測された物理量との差に基づく制御量が目標値となるように前記制御対象を制御する制御ステップと、
    を含む制御方法。
19. コンピュータに、
    制御対象を示す解析領域を離散化した複数の要素ごとの第１時刻の第１物理量と、前記制御対象の負荷条件、境界条件、材料の条件、および、構造条件の少なくとも一部を表す条件データと、を含む複数の入力データを取得する取得ステップと、
    複数の前記入力データを推定モデルに入力して得られる、複数の前記要素のエネルギを表すエネルギ汎関数の値を用いた変分原理により、前記第１時刻より後の第２時刻の複数の前記要素ごとの第２物理量を算出する算出ステップと、
    前記第２物理量と計測された物理量との差に基づく制御量が目標値となるように前記制御対象を制御する制御ステップと、
    を実行させるためのプログラム。
20. 励起力を発生させる自励振動または自励回転を行う物体である制御対象の第１時刻の第１物理量と、前記制御対象の負荷条件、境界条件、材料の条件、および、構造条件の少なくとも一部を表す条件データと、を含む入力データを取得する取得部と、
    前記入力データを推定モデルに入力して得られる、前記制御対象のエネルギを表すエネルギ汎関数の値を用いた変分原理により、前記第１時刻より後の第２時刻の前記制御対象の第２物理量を算出する算出部と、
    前記第２物理量と計測された物理量との差に基づく制御量であって、前記励起力と前記制御対象の動作に影響する物理量との関係を示す前記制御量が、前記自励振動または前記自励回転を誘起する前記励起力と前記物理量との関係に基づいて定められる目標値となるように前記制御対象を制御する制御部と、
    を備える制御装置。

Images