Activate で Modelica を使ってみる（その 12）- 最適化
Activate で Modelica を使ってみる（その 12）
最適化
前回でも書きましたが、パラメータ化されていることが多く、かつ、有限要素モデルなどに比較して規模が小さい Modelica モデルは最適化との相性が抜群です。今回は、前回で紹介したシミュレーションの反復手法をベースにした最適化の実行方法を紹介します。
Activate での最適化手法
Activate 上で作成されたモデルに対して最適化を適用する方法は以下の 4 種類用意されています。
- Activate ライブラリ内の最適化用ブロックをダイアグラムに組み込む
- OML の最適化コマンドを用いて Activate をバッチで実行する。
- Activate の最適化ユーティリティを用いる（内部的には上の 2 が実行される）。
- 最適化ツールである HyperStudy を用い、Activate をソルバーとしてバッチ実行する。
もっとも汎用的で多くの最適化問題に柔軟に対応できるのは 4) の手法です。次に汎用的なのは 2) の手法ですが、OML のある程度深い知識が必要となってきます。3) は 2) の手法を半自動化し使いやすくしたものですが、使える機能や柔軟性は 2) より限定されます。1) は Activate 内で完結でき、OML の知識もそれほど必要ではないため最も使いやすい方法であると言えます。ただ、この手法では制約条件は設定できません。
ここでは 1) に着目して詳しく紹介します。その他の手法については、Altair Community の以下の投稿を参照するか、下記の文献（英文）をご覧ください。
1) “1Dモデルのパラメータを最適化する その1（最適化ブロック）”
https://community.altair.com/community?id=community_question&sys_id=10664cf61b2bd0908017dc61ec4bcbc1&view_source=searchResult
3) “1Dモデルのパラメータを最適化する その2（最適化ツール）”
https://community.altair.com/community?id=community_question&sys_id=58664cf61b2bd0908017dc61ec4bcbc1&view_source=searchResult
4) “1Dモデルのパラメータを最適化する その3（HyperStudy）”
https://community.altair.com/community?id=community_question&sys_id=0c664cf61b2bd0908017dc61ec4bcb9b&view_source=searchResult
1) ～ 4) “Altair Activate Extended Definitions”（Activate のヘルプ内にあります）のChapter 6 “Optimization”
1) Stephen L. Campbell、Ramine Nikoukhah：”Modeling and Simulation with Compose and Activate” 8.2 章 “BobyqaOpt Block”
BobyqaOpt ブロック
最適化のために用意されているのが BobyqaOpt ブロックです。最適化では最適解を収束計算で求めるためシミュレーションを複数回実行する必要があります。BobyqaOpt は、前回（”その11”）で紹介した End ブロックと同様のシミュレーション反復実行機能を内蔵しており、最適解が見つかるまで自動的にシミュレーションを反復実行し続けます。
BobyqaOpt は以下のような入出力ポートを持っています。
入力
p0：設計変数（ベクトル）
cost：目的（費用）関数
アクティベーション入力（上面のポート）
出力
p：アップデートされた設計変数（ベクトル）
cost_opt：アップデートされた目的関数
end：アクティベーション出力（下面のポート）
BobyqaOpt はアクティベーションされた時点で、シミュレーションからの計算結果である目的関数の値を costポートから読み込み、最適化計算を１ステップ分実行して設計変数を p0 から p へアップデートします。また p を用いた時の目的関数の値も “参考に” 出力します（計算に使われることはありません）。BobyqaOpt はその実行後に収束判定を毎回行い、収束した（すなわち、最適解が得られた）と判断した場合はアクティベーション出力ポートからアクティベーション信号を出力します。
この入出力ポートの意味を理解するための非常に簡単なモデル（添付：MathExpress_opt.scm）を以下のように作成しました。
このモデルの目的は単に y = 2x2 - 4x + 5 という簡単な２次方程式の最小値（x = 1 の時に y = 3）を求める、というものです。ここでは、MathExpression ブロックで計算されたその多項式の結果（すなわち、上式の y）を最小化します。設計変数は当然ながら x です。
BobyqaOpt は、アップデートした設計変数（x）の値を ToBase ブロックを通して一旦 OML ワークスペースに保存します。次の繰り返し時に（すなわち、次に BobyqaOpt がアクティベートされた時に）、前に保存された x の値を “ダイアグラム” の “コンテキスト” のスクリプトで読み込みます（後述）。また同時にその x の値で計算された目的関数（y）の値も cost ポートから入力されるので、x の変化に応じた y の変化が分かることになり、それを最適化計算に利用します。
この BobyqaOpt ブロックをアクティベートしているのは EventGenerate ブロックです。このモデルでの EventGenerate は 1 秒おきにアクティベーション信号を出力しており、すなわち、1 秒おきに BobyqaOpt を起動することになります。通常は、シミュレーションが終了し目的関数が計算できた時点で BobyqaOpt を起動しますが、このモデルは “動的な” モデルではなくシミュレーション結果が時間的に変化しないためこのような設定にしています。極端な話、1 ステップ計算すればこの場合は十分です。
このモデルでは end ポートに End ブロック（Stop オプション）が接続されているため、BobyqaOpt が最適解を見つけた時点で最適化の反復計算は終了します。
目的関数を計算するシミュレーション部分は Constant_Input と MathExpression の二つのブロックのみで構成されています（Modelica を紹介するシリーズでありながら Modelica 部分が含まれておらず申し訳ありません）。この Constant_Input ブロックの定数 x の値は “ダイアグラム” の “コンテキスト” で以下のように設定されています。
この
x = GetFromBase('x_base',0)
という文では、ToBase ブロックによって OML ワークスペースに値が書き込まれた x_base という変数から値を取り込んで x に代入しています。要するに、GetFromBase という OML コマンドは FromBase ブロックと同じ働きをしています。このように、設計変数をシミュレーションのモデルに組み込む場合は、GetFromBase コマンドを用いて繰り返し計算ごとに変化する設計変数の値をモデルに取り込む必要があります。
この最適化計算の結果は下図のように Display_p ブロックに設計変数の（x）の結果が、Display_cot ブロックに目的関数（y）の結果が表示されており、最適化が正しく実行されていることが分かります。
なお、BobyqaOpt に組み込まれている最適化手法については詳しくは以下の文献を参照してください。
M.J.D. Powell：”The BOBYQA algorithm for bound constrained optimization without derivatives”
例１：バネマスモデルの振幅の最大化
まず1つ目の事例として。前の回（”その11”）で用いた、バネマスモデルの振動モデルを用いて最適化手順を紹介します。ここでの目的関数は定常振動の振幅の最大化であり、設計変数は入力荷重である正弦波の振動数です。前回は、設計変数である荷重の振動数を徐々に変化させながら結果のプロット上でピークが発生する振動数を観察しましたが、今回はその振動数（要するに、共振振動数）とその時のピーク値を最適化で求めます。
このモデルは以下のような構成になります（添付：SpringMass_Peak_Opt.scm）。
シミュレーション部分は前回と同じです。ここで計算された振幅を目的関数として BobyqaOpt の cost ポートに出力します。また、この問題での設計変数は荷重の振動数（wf）です。BobyqaOpt によりアップデートされていく wf を以下のように “コンテキスト” で指定してモデルに取り込んでいます。
このモデルでの最適値（すなわち、共振）は減衰固有振動数 wd（= 1.2）で起こることが予想されるため計算範囲（wf_min～wf_max）を 0.01～2.5 と設定しました。これは BobyqaOpt のブロックダイアログで以下のように設定されています。
ここでは、合わせて、繰り返し回数の上限値を 100 と設定し、かつ、最小値ではなく最大値を計算するというオプションをオンにしています。
計算結果は以下のように
振幅の最大値が 0.85、その時の振動数が 1.17 となり、前の回の結果（以下に再掲）の結果と一致することが分かります。ただ、ピークが発生する振動数が減衰固有振動数（= 1.2）とずれていますが、この原因は調べきれていません。
例２：バネマスモデルの質量と剛性の同定
ここでは最適化ではなく、最適化の機能を用いて同定を行います。同定とは、計算結果を、あらかじめ用意されたデータ（例えば、実験結果）に一致させるようにモデルのパラメータ（例えば、モデルの特性値）を調整することです。
以下のように、特性（剛性と質量。減衰係数はこの２つの特性値から計算）値が異なる２つのバネマスモデルを含むモデルを用意し（添付：SpringMass_Ident_KMB.scm）、２つの結果が一致するように片方の特性値（ここでは下側のK, M）を調整する、という同定を行います。
同定では下側のバネマスの特性値を変更して上側の結果に一致させることを目的としています。これを行うために、各時刻での２つの結果の差（Sum ブロック）の２乗（Power ブロック）を合計（Integral ブロック）し、それを最小化させることを目的関数とします。合計は必ず正の値であるため、最小化することにより合計がゼロに近づきます。各時刻での結果の差の２乗も正の値なので、合計をゼロにするためには各時刻の差もゼロに近づくことになり、２つのカーブが一致する方向に向かうことになります。
このモデルの “コンテキスト” のスクリプトは以下のように設定しています。
1 行目から 3 行目では上側のバネマスモデルのバネ定数、質量、減衰を定義しています。この値が、8、9 行目で定義されている下側のバネ定数と質量の目標値になります。7 行目では KM_base という変数（列ベクトルで初期値は [10.0; 10.0]）名の OML ワークスペースから値を読み込んで、KM というベクトル変数に代入しています。また、5, 6 行目では、 KM_base の最小値と最大値を定義しています。この値は、以下のように BobyqaOpt ブロックで設計変数の最小値／最大値として使用されます。
最後に、8 行目でベクトル KM の1 要素目を K に、2 要素目を M に代入し、B はその K と M から計算しています。
結果は以下の Display_KM が示すように、設計変数の最終値が K=3 および M=2 となって、上側のバネマスモデルに一致したことが分かります。
結果のプロットも以下のように一致しています。
例3：実験結果との同定
同定の手法がよく使われるのは、実験結果にシミュレーションの結果を一致させるためにシミュレーションモデルの特性値を調整する場合です。ここでは、例２のモデルを少し変更してその手順を紹介します。
この場合のモデルは以下のようになります（添付：SpringMass_Ident_csv.scm）。
例２のモデルの上側のバネマスモデルが csv ファイルの入力（FromCSV ブロック）に置き換えられているだけです。この csv ファイルは、例えば、実験の結果であり、これが同定のターゲットになるということになります。
このモデルの ”初期化” スクリプトは以下のように例2と同様の構成になっています。
ここでは結果を例２に合わせるために、例２の上側のバネマスを取り出した以下のモデル（添付：SpringMass_Ident_Base.scm）の結果（添付：Output.csv）を csv ファイルとして取り込むことにします。つまり、設計変数の最適値は K=3、M=2 になることを想定しています。
結果は以下の通りとなり、期待した結果が得られています。
例4：バネマスモデルの PID 制御係数の決定
PID 制御を用いた際にその最適な係数を決めるためには経験や洞察が必要となり労力を要することが多くあります。そのため、最適化手法を PID 係数の決定に用いることがよくあります。
ここでは簡単な PID 制御の例として、変位の時間的な変化パターンをターゲットとしてあらかじめ設定しておき、バネマスモデルの結果変位がそのカーブに追随するような制御を考えます。この PID の係数を最適化するために、上の例題でも用いた同定手法を使って、ターゲットカーブと結果変位の差を最小化する形で応答をターゲットカーブに近づけます。
ただ、最適化された PID の係数は設定された最大値（アクチェーターの能力で決定される値）で決められてしまうことも多く、これを避けるためには、オーバーシュートや立ち上がり時間等の制御時の応答特性を制約条件として最適化問題に加えることが望ましいと考えられます。残念ながら、ここで紹介している BobyqaOpt を用いる手法では制約条件は考慮できないため、その場合は、制約条件を考慮できる他の手法（最初に紹介した 2）～ 4) の最適化手法）を使用する必要があります。
ということで、ここでの例題は最適な PID 係数を求めるための最良の手法を紹介しているわけではなく、あくまで PID 制御のモデル化とその係数の最適化の考え方を理解するためのものであることをご理解ください。
以下がいつものバネマスモデルから減衰を除いたモデルに対して PID 制御を適用したモデルです（添付：SpringMass_PID_Base.scm）。
応答量（変位）の時間的変化のターゲットとして SquareWaveGenerator ブロックで以下のようなパルスカーブを入力しています。
PID スーパーブロックの中身は以下の通りのダイアグラムであり（ライブラリにある PID ブロックを使っても構いません）その係数は以下の “初期化” スクリプトのように “適当に” 設定しています。
このシミュレーション部分に、前の例題でも紹介した、同定のための最適化部分を加えたのが以下のモデルです（添付：SpringMass_PID_opt.scm）。
結果は上の図の Display_PID ブロックに示されているとおりであり、その時の変位結果は以下のような時間的変化となりました。
今回は、最適化のためのブロックである BobyqaOpt を Modelica モデルに適用した事例を紹介しました。ここ（”その 12”）までは Modelica 標準ライブラリの中の Spring、Mass、damper の３つのコンポーネントしか使用してきませんでしたが、次回は Translational ライブラリに含まれている接触／摩擦系のコンポーネントである、ElastoGap、SupportFriction、Brake、MassWithStopAndFriction の各コンポーネントについて、その機能や仕組みを簡単なサンプルとともに紹介します。