スパースモデルの基礎

モデルオブジェクト

sparss オブジェクトと mechss オブジェクトを使用して、スパースな 1 次システムと 2 次システムをそれぞれ表します。このようなスパースモデルは有限要素解析 (FEA) から発生し、構造解析、流体の流れ、熱伝達、電磁などの分野で役立ちます。FEA には有限要素法 (FEM) を使用した問題の解析が含まれます。この方法では、大規模なシステムが多数の小さいコンポーネントまたは有限要素 (FE) に分割され、その後、個別に解析されます。これらの多数のコンポーネントを組み合わせると大規模なスパースモデルになり、ss などの従来の密なモデルオブジェクトで表現するには計算量が多く、非効率的です。

スパース行列は大部分がゼロの double または logical データの効率的なストレージを提供します。非スパース行列 (密行列) は、値に関係なく個々の要素をメモリに保存しますが、スパース行列は 0 以外の要素とその行インデックスだけを保存します。このため、スパース行列を使用すると、データストレージに必要なメモリ量を大幅に削減できます。詳細については、スパース行列の計算上の利点を参照してください。

次の表は、表現できるスパースモデルのタイプを示しています。

モデルタイプ	数学的表現	モデルオブジェクト
連続時間スパース 1 次モデル	$\begin{array}{r} E \frac{d x}{d t} = A x (t) + B u (t) \\ y (t) = C x (t) + D u (t) \end{array}$	`sparss`
離散時間スパース 1 次モデル	$\begin{array}{r} E x [k + 1] = A x [k] + B u [k] \\ y [k] = C x [k] + D u [k] \end{array}$	`sparss`
連続時間スパース 2 次モデル	$\begin{array}{r} M \ddot{q} (t) + C \dot{q} (t) + K q (t) = B u (t) \\ y (t) = F q (t) + G \dot{q} (t) + D u (t) \end{array}$	`mechss`
離散時間スパース 2 次モデル	$\begin{array}{r} M q [k + 2] + C q [k + 1] + K q [k] = B u [k] \\ y [k] = F q [k] + G q [k + 1] + D u [k] \end{array}$	`mechss`

sparssdata および mechssdata を使用してスパース形式のモデル行列にアクセスできます。spy コマンドを使用して 1 次および 2 次両方のモデル行列のスパース性を可視化することもできます。

また、任意の非 FRD モデルを sparss または mechss オブジェクトにそれぞれ変換することもできます。逆に、full コマンドを使用すると、スパースモデルを密度の高いストレージ ss に変換できます。大規模なスパースモデルでは、密度の高いストレージに変換するとメモリ使用量が大きくなりパフォーマンスが低下するため推奨されません。

スパースモデルの組み合わせ

信号ベースの接続

モデルの相互接続にリストされている標準の信号ベースの接続はすべて、スパースモデルオブジェクトでサポートされています。信号を使用する相互接続モデルでは、制御システム用のモデルを作成できます。制御システムは、フィードバック構成で接続されているプラントやコントローラーなど、複数の相互接続されたコンポーネントを含んでいるブロック線図として概念化できます。モデル演算または相互接続コマンドを使用して、これらの各コンポーネントのモデルを 1 つのモデルに統合し、ブロック線図で全体を表現します。

sparss モデルの接続の直列接続とフィードバック接続では、内部信号を取り除くと、A 行列に望ましくない非スパース要素が発生します。これを回避するために、A のスパース性に影響しない信号のみが取り除かれ、残りの信号は状態ベクトルに追加されます。通常、これにより、A 行列に次の図に示すブロック矢印構造がある相互接続の微分代数方程式 (DAE) モデルが生成されます。

ここでは、各対角ブロックはスパースモデルのサブコンポーネントです。最後の行と列は Interface グループと Signal グループを組み合わせてコンポーネント間のすべての結合と接続を取得します。

2 次 mechss モデルオブジェクトに対しても同様の規則が適用されます。showStateInfo を使用して、コンポーネント、インターフェイス、および信号に対する状態ベクトルの分割の概要を出力できます。

xsort コマンドを使用して、状態の分割に基づいて状態を並べ替えます。並べ替えられたモデルでは、すべてのコンポーネントが最初に表示され、インターフェイスが続き、その後にすべての内部信号の単一グループが続きます。

例については、sparss および mechss のリファレンスページを参照してください。

物理結合

interface コマンドを使用して mechss モデルのコンポーネント間の物理的結合を指定できます。interface は "デュアルアセンブリ" 形式を使用し、自由度 (DoF) のグローバルセットが保持され、インターフェイスにおける結合制約を記述することで、アセンブリが作成されます。サブ構造体 S₁ の DOF N₁ とサブ構造体 S₂ の DOF N₂ 間の剛結合の場合は、次が含まれます。

変位の一致: q(N₁) = q(N₂)
作用/反作用の原理: g(N₁) + g(N₂) = 0。ここで、g(N₁) は S₁ に S₂ が及ぼす力、および g(N₂) は S₂ に S₁ が及ぼす力です。

これらの制約は、次の方程式に集約できます。

$\begin{array}{l} M \ddot{q} + C \dot{q} + K q = B u + g, H q = 0, g = - H^{T} λ, \\ y = F q + G \dot{q} + D u, \end{array}$

ここで、g はインターフェイスの力のベクトル、H は適切な DOF 選択行列です。

interface コマンドを使用して DOF 間の非剛結合を指定することもできます。

詳細については、interface およびモデルコンポーネントの剛体アセンブリを参照してください。

タイプが異なるモデルの組み合わせ

次の優先順位ルールは、タイプが異なるモデルを組み合わせるときに適用されます。

スパースモデルを FRD モデルと組み合わせると frd モデルオブジェクトが生成される
スパースモデルを tf、ss、および zpk などの任意の非 FRD モデルと組み合わせるとスパースモデルオブジェクトが生成される
sparss モデルと mechss モデルを組み合わせると mechss モデルオブジェクトが生成される
現在、スパースモデルを一般化モデルまたは不確かさを持つモデルである genss および uss (Robust Control Toolbox) と組み合わせることはできない

例については、mechss のリファレンスページを参照してください。

時間領域解析

時間領域および周波数領域解析にリストされている標準の時間領域解析関数はすべて、sparss モデルオブジェクトと mechss モデルオブジェクトでサポートされています。

スパースモデルに対する時間領域応答関数を使用する場合は、最終時間または時間ベクトルを指定する必要があります。以下に例を示します。

tf = 10;
step(sys,tf)

t = 0:0.001:1;
initial(sys,x0,t)

時間応答関数は、特に大規模なスパースモデル用に開発されたカスタムの固定ステップ DAE ソルバーに依存します。DAE ソルバータイプを構成し、sparss モデルオブジェクトと mechss モデルオブジェクトの SolverOptions プロパティを使用することで並列計算を有効化できます。各入力チャネルの応答は並列でシミュレーションされるため、並列計算は多入力モデルに対する step または impulse 応答シミュレーションを加速するために使用できます。並列計算の有効化には Parallel Computing Toolbox™ ライセンスが必要です。

使用可能なソルバーオプションを、次の表に示します。

DAE ソルバー	説明	使用方法
`'trbdf2'`[2]	精度が `o(h^2)` の固定ステップソルバー。ここで、`h` はステップサイズです。`'trbdf2'` は `'hht'` より 50% 計算量が多くなります。これは `sparss` モデルと `mechss` モデル両方に対する既定の DAE ソルバーオプションです。	`sparss` モデルと `mechss` モデルの両方で使用できます。
`'trbdf3'`	精度が `o(h^3)` の固定ステップソルバー。`'trbdf3'` は `'trbdf2'` より 50% 計算量が多くなります。	`sparss` モデルと `mechss` モデルの両方で使用できます。
`'hht'`[1]	精度が `o(h^2)` の固定ステップソルバー。ここで、`h` はステップサイズです。`'hht'` は最速ですが、システムのインパルス応答など、初期の加速度が高いと問題が発生する場合があります。	`mechss` モデルでのみ使用可能です。

並列計算とソルバーの選択を有効にするには、次の構文を使用します。

sys.SolverOptions.UseParallel = true;
sys.SolverOptions.DAESolver = 'trbdf3';

例については、片持ちビームの線形解析を参照してください。

周波数領域解析

時間領域および周波数領域解析にリストされている標準の周波数領域解析関数はすべて、sparss モデルオブジェクトと mechss モデルオブジェクトでサポートされています。周波数応答計算の場合、並列計算を有効にすると、使用可能な複数のワーカー間で一連の周波数が分散され、応答計算が高速化されます。並列計算の有効化には Parallel Computing Toolbox ライセンスが必要です。

スパースモデルに対する周波数応答関数を使用する場合は、周波数ベクトルを指定する必要があります。以下に例を示します。

w = logspace(0,8,1000);
bode(sys,w)
sigma(sys,w)

例については、片持ちビームの線形解析を参照してください。

連続と離散の変換

c2d、d2c および d2d を使用して連続時間と離散時間の間の変換を行い、sparss モデルをリサンプリングできます。次の表に、sparss モデルに使用可能な方法を示します。

メソッド	説明	使用方法
`'tustin'`	双一次 (Tustin) 近似	関数 `c2d`、`d2c` および `d2d` で使用できます。
`'damped'`	TRBDF2[2]式に基づいて減衰された (Tustin) 近似。`'tustin'` メソッドとは対照的に、このメソッドでは減衰が無限大に提供されます。つまり、高周波数のダイナミクスはフィルターで除外され、離散モデルで `z = -1` に近いモードの累積に寄与しません (数値的不安定性の原因)。	`c2d` でのみ使用できます。

'tustin' メソッドで c2d、d2c、および d2d を使用して連続時間と離散時間の間の変換を行い、mechss モデルをリサンプリングできます。'tustin' メソッドは Tustin 離散化の 2 次形式を計算します。これは、mechss モデルと同等の 1 次 sparss に Tustin を適用することと等価です。対称性が維持されるため、モーダル近似にはこの形式の方が適しています。 (R2024a 以降)

スパース線形化

Simulink モデルの線形化

Descriptor State-Space (Simulink) または Sparse Second Order ブロックがある場合、Simulink^® モデルからスパース線形化モデルを取得できます。

Sparse Second Order ブロックは常に mechss モデルに線形化されるように構成されます。その結果、線形化されたモデル全体は、このブロックが存在する場合に 2 次スパースモデルになります。
Descriptor State-Space (Simulink) ブロックで [スパースモデルに線形化する] オプションをオンにして sparss モデルを線形化します。これは、Descriptor State-Space (Simulink) で LinearizerToSparse パラメーターを 'on' に設定して実行することもできます。既定では、[スパースモデルに線形化する] ボックスはオフになっていて、LinearizerToSparse は 'off' です。ブロックレベルの LinearizeToSparse 設定は、SCDBlockLinearizationSpecification ブロックパラメーターまたは linearize (Simulink Control Design) に対する blocksub 入力のいずれかを使用してブロックの置換を指定すると無視されます。

線形化ワークフローには、SimulinkControl Design™ ソフトウェアが必要です。

例については、スパース 2 次モデルオブジェクトへの Simulink モデルの線形化を参照してください。

構造モデルまたは熱 PDE モデルの線形化

関数 linearize (Partial Differential Equation Toolbox) を使用して、構造モデルまたは熱 PDE モデルからスパース線形化モデルを取得できます。

構造解析モデルの場合、linearize は mechss モデルを抽出。
熱解析モデルの場合、linearize は sparss モデルを抽出。

linearizeInput (Partial Differential Equation Toolbox) を使用して、PDE モデル内の境界条件、負荷、または内部熱源について線形化モデルの入力を指定します。linearizeOutput (Partial Differential Equation Toolbox) を使用して、セル (3 次元ジオメトリの場合のみ)、面、エッジ、頂点などのジオメトリの領域について線形化モデルの出力を指定します。

線形化ワークフローには、Partial Differential Equation Toolbox™ ソフトウェアが必要です。

例については、片持ちビームの線形解析と音叉の線形解析を参照してください。

モデル次数の低次元化

R2023b 以降

スパース状態空間モデルの低次の近似を計算するには、関数 reducespec を使用します。reducespec は Control System Toolbox™ ソフトウェアにおけるモデル次数の低次元化ワークフローのエントリポイントです。完全なワークフローについては、Task-Based Model Order Reduction Workflowを参照してください。

このソフトウェアでは、以下の手法を使用したスパースモデルの次数の低次元化がサポートされています。

平衡化打ち切り — 寄与度の低い状態を破棄して低次の近似を得る。
モーダル打ち切り — 望ましくないモードを破棄して低次の近似を得る。
固有直行分解 — シミュレーションデータを用いて近似の平衡化打ち切りを実行し、状態ベクトルの主要なモード (主成分) を計算および抽出する。 (R2025a 以降)
周波数応答フィッティング — 低次のモデルを大規模なスパース状態空間モデルの周波数応答データに当てはめる。 (R2025a 以降)
極零点打ち切り — 任意の低周波数帯域でのスパース状態空間モデルの零点-極-ゲインの近似を得る (R2025a 以降)

例については、SparseBalancedTruncation、SparseModalTruncation、ProperOrthogonalDecomposition、FrequencyResponseFitting、および SparseZeroPoleTruncation を参照してください。

さらに、tf および zpk を使用して、特定の低周波数帯域での元のスパースモデルの打ち切られた tf または zpk モデルの近似を得ることができます。 (R2025a 以降)

その他のサポートされる機能

また、以下の機能が現在、スパースモデルでサポートされています。

dcgain を使用した低周波ゲインの計算
evalfr および freqresp を使用した周波数応答の評価
inv を使用したモデルの反転
imp2exp を使用した陰的-陽的な関係の変換
pade を使用したパデ近似
入力、出力、および内部遅延の指定
I/O の選択およびインデックス付け
pole、damp、および isstable を使用した極のサブセット、減衰、固有振動数の計算、および指定の周波数帯域における安定性の確認。 (R2025a 以降)

制限

以下の操作は現在、スパースモデルでサポートされていません。

安定余裕の計算
補償器の設計および調整

full コマンドを使用して小規模なスパースモデルを密度の高いストレージ ss に変換し、前述の操作を実行します。大規模なスパースモデルで密度の高いストレージに変換すると利用可能なメモリが飽和し、パフォーマンスが大幅に低下する場合があるため、推奨されません。

スパース線形化には次の制限があります。

スパース線形化は、調整可能なモデルまたは不確かさをもつモデルに関連するブロックの置換と互換性がありません。genss モデルまたは uss (Robust Control Toolbox) モデルに線形化する場合、[スパースモデルに線形化します] オプションをオフにしなければなりません。
Descriptor State-Space (Simulink) または Sparse Second Order ブロックが存在する場合、スナップショットの線形化は動作しない可能性があります。スナップショットの線形化にはシミュレーションが必要なため、シミュレーション機能は現在 Descriptor State-Space (Simulink) ブロックと Sparse Second Order ブロックに制限されています。

参照

[1] H. Hilber, T. Hughes & R. Taylor. " Improved numerical dissipation for time integration algorithms in structural dynamics." Earthquake Engineering and Structural Dynamics, vol. 5, no. 3, pp. 283-292, 1977.

[2] M. Hosea and L. Shampine. "Analysis and implementation of TR-BDF2." Applied Numerical Mathematics, vol. 20, no. 1-2, pp. 21-37, 1996.

参考