Dual-monitor engineering workspace displaying Simulink control block diagrams and Bode plots next to an active suspension hardware testbench setup for LQR design

Simulinkマスターガイド:制御システムシミュレーションとLQR設計の実践ガイド

  • 著者: ダウウェイ・ビークルのCEO、ジョニー・リュー
  • 公開日: 2026年4月29日
  • カテゴリ: 制御システム/自動車工学/シミュレーション

著者について

私はジョニー・リウです。 ダウウェイ車両. 私たちはクリーンな動力伝達システムと高度な車両ダイナミクスを開発しています。私は20年以上にわたり、アクティブサスペンション、ドライブバイワイヤステアリング、シャシーの安全性に関する研究開発に携わってきました。私たちの工場では、実際に部品を製造する前に必ず画面上でテストを行います。

Simulinkは、まさにそれを可能にしてくれます。金属を曲げる前に、安全性のテストを実行したり、新しい制御アルゴリズムを試したりできるのです。このガイドでは、動的システムのモデリング、解析、制御に私たちが実際に使用しているワークフローを詳しく解説します。

パート1:制御システム理論

ブロックを動かし始める前に、数学を理解しなければなりません。制御理論では、モデルを主に3つの形式に分類します。

1.1 数理モデリング

動的モデルは、入力が加えられた際に、システムの状態と出力が時間とともにどのように変化するかを示します。

1. 投入産出モデル

  • 微分方程式: これらは、システムの物理法則を時間領域で記述したものです。例としては、質量-ばね系や電気回路などが挙げられます。
  • 伝達関数(TF): これらは、すべての初期条件がゼロであると仮定して、線形時不変(LTI)微分方程式にラプラス変換を適用することによって得られます。$$G(s) = \frac{Y(s)}{U(s)} = \frac{b_m s^m + b_{m-1} s^{m-1} + \dots + b_0}{a_n s^n + a_{n-1} s^{n-1} + \dots + a_0}$$
  • 周波数特性: これらは、伝達関数中の$s$と$j\omega$を入れ替えることで得られます。これにより、システムがさまざまな速度の正弦波入力にどのように応答するかがわかります。

2. 状態空間モデル

現代の制御理論では、状態空間モデルを使用します。システムを隠された箱のように扱う代わりに、状態と呼ばれる内部変数を追跡します。これらの状態は、任意の時点でのシステムのエネルギーと挙動を定義します。$$\dot{x}(t) = Ax(t) + Bu(t)$$$$y(t) = Cx(t) + Du(t)$$

どこ:

  • $x(t)$は状態ベクトルです。
  • $u(t)$は入力ベクトルです。
  • $y(t)$は出力ベクトルです。
  • $A$はシステム行列(内部ダイナミクスを支配する行列)である。
  • $B$は入力行列です(入力が状態にどのように影響するかを示します)。
  • $C$は出力行列です(状態を測定された出力にマッピングします)。
  • $D$は直接フィードスルー行列です。

伝達関数と状態空間モデルの比較は以下のとおりです。

特徴伝達関数(TF)状態空間(SS)
システム互換性単入力単出力(SISO)および線形時不変(LTI)構成でのみ動作します。マルチ入力マルチ出力(MIMO)、非線形、および時変構成に対応します。
独自性個性的 任意の入力出力ペアに対して。非一意. 異なる状態を選択することで、行列が変わります。
完全のみ表示 制御可能 部品。極零点相殺により、隠れたダイナミクスが隠される。すべてを表示します。 制御可能なものと制御不可能なものの両方 部分、さらに観測可能な部分と観測不可能な部分。
計算周波数領域代数を使用します。コンピュータおよび線形代数ソルバーに適しています。

3. ブロック図

物理システムを接続するために、以下の3つのシンプルな構成を採用しています。

  • シリーズ(カスケード): 伝達関数は次のように乗算されます:$G_{\text{total}}(s) = G_1(s) \cdot G_2(s)$
  • 平行: 伝達関数は加算される:$G_{\text{total}}(s) = G_1(s) + G_2(s)$
  • フィードバック: 閉ループフィードバック $H(s)$ の式は次のとおりです。$$G_{\text{closed}}(s) = \frac{G(s)}{1 \mp G(s)H(s)}$$

システムパフォーマンス指標

ステップテストを実行する際、過渡応答性能は4つの数値で評価します。

  • オーバーシュート($M_p$): ピーク値を最終的な定常値と比較した値をパーセンテージで表すと次のようになります。$$M_p = \frac{y(t_p) – y(\infty)}{y(\infty)} \times 100\%$$
  • 沈降時間($t_s$): システムが安定し、最終値の小さな誤差範囲(通常は2%または5%)内に収まるまでの時間。
  • ピークタイム($t_p$): 最初のオーバーシュートピークに到達するまでの時間。
  • 立ち上がり時間($t_r$): 最終値の10%から90%(または0%から100%)に達するまでにかかる時間。

1.2 分析方法

コントローラーを構築する前に、システムの安定性をテストする必要があります。

1. 時間領域解析

システムに標準テスト信号を適用します。これには以下が含まれます。 単位インパルス単位ステップユニットランプ単位加速度、 または 単位正弦波

  • 安定性ルール: LTIシステムは、そのすべての極(分母の根)の実部が負である場合に安定である。極が1つでもs平面の右側に位置すると、システムは崩壊する。
  • MATLABツール:
    • zpk: ゼロ極ゲインシステムを定義します。
    • ルーツ分母の根を求めて、極の位置を確認します。

2. 根の軌跡分析

ゲイン(K)を0から∞まで変化させたときに、閉ループ極がs平面内でどのように移動するかをプロットします。

  • プロット上のどの点においても、以下の2つのルールが一致しなければならない。
    • マグニチュード条件: $|KG(s)H(s)| = 1$
    • 位相条件: $\angle G(s)H(s) = (2k+1)\pi$
  • MATLABツール: 使用 rlocus(sys) 線を引くには、常に 軸は等しい 直後に配置することで、水平軸と垂直軸の目盛りが同じになります。これにより、角度が正しく表示されます。

3. 周波数領域解析

システムが様々な周波数の定常正弦波をどのように処理するかを検証します。

  • プロット:
    • ナイキスト線図: 周波数 ($\omega$) が $-\infty$ から $+\infty$ まで変化するときの $G(j\omega)$ の極座標プロット。
    • ボーデ線図: 2つの別々のグラフ。上のグラフはデシベル単位の振幅($20\log_{10}|G(j\omega)|$)を示し、下のグラフは位相角を​​度単位で示しています。
    • ニコルズ・プロット: これにより、2つのボード線図が1つの画面に表示されます。周波数は線上の隠れたマーカーとして使用され、デシベル値が位相角に対して直接プロットされます。
  • MATLABツール: 使用 ナイキスト()ボード線()、 そして ニコルス()
  • 主要数値: 共振ピーク($M_r$)、共振周波数($\omega_r$)、帯域幅、およびゼロ周波数ゲイン。
  • 安定性チェック(ナイキストの定理): ナイキスト線が点 $(-1, j0)$ を反時計回りに $N$ 回周回する場合、システムは安定していると言えます。ここで、$N$ は不安定な開ループ極の数 ($P$) に等しくなります。
    • 利益率($G_M$): システムが不安定になる前に、どれだけのゲインを追加できるか。
    • 位相余裕($\Phi_M$): システムを不安定にするために必要な追加の位相遅延。
    • MATLABコマンド: margin(sys) これらの安全マージンを自動的に計算します。

4. 状態空間解析

行列は、多数の入力と出力を持つ複雑なシステムを研究するために用いられます。

  • 標準様式: 変換行列($x = Pz$)を用いることで、単一の伝達関数を異なる状態空間レイアウトにマッピングできます。主な構成は以下の4つです。
    1. 制御可能な正準形式 (状態フィードバック制御器の設計に役立ちます。)
    2. 観測可能な正準形式 (状態推定器の構築に役立ちます。)
    3. 対角線正準形 (対角線に沿った非結合状態。極が明確に区別できる場合に使用。)
    4. ジョーダン正準形式 (極が繰り返される場合に使用)
    • MATLABコマンド: ss()tf2ss()zp2ss()キャノン()、 そして ヨルダン()
  • リャプノフ安定性: システムが安定であるのは、時間微分 $\dot{V}(x)$ が常に負となる正定値エネルギー関数 $V(x)$ を見つけることができる場合である。線形システムの場合、 リアプノフ方程式: $$A^TP + PA = -Q$$ 任意の正定値行列 $Q$ を選択して一意の正定値行列 $P$ を得ると、システムは安定します。
    • MATLABコマンド: lyap(A, Q) または lyap2() この方程式を解いてください。

それでは、Simulinkでこれらのシステムを構築するために使用するブロックを見ていきましょう。

2.1 標準ブロックライブラリ

1. 連続サブライブラリ

ここでは、連続時間モデルを構築します。

  • 状態空間ブロック: $\dot{x} = Ax+Bu$ と $y=Cx+Du$ を計算します。 行列A、B、C、およびD開始状態ベクトル $x_0$ と共に。
  • 転送Fcnブロック: あなたは 分子(num) そして 分母(den) $s$ の降べき乗としてのベクトル。
  • ゼロ極ゲインブロック: 配列を入力します ゼロ(z)極(p)、ゲイン定数 (K)
  • PIDコントローラーブロック(1自由度および2自由度): これらのブロックは、比例項、積分項、および微分項を計算します。
    • パラメータ: コントローラの種類、形式(並列または理想)、時間領域、係数($P、I、D$)、フィルタ定数($N$)を設定します。開始状態を設定したり、上限と下限をオンにしたりできます。 飽和限界 積分器のワインドアップを停止するには、[データ型] タブで境界を設定できます。デフォルトでは、制限チェックなしの継承ルールが使用されます。
    • 2自由度PID制御と1自由度PID制御の違い: 標準的な1自由度コントローラは、誤差$e = r – y$のみを考慮します。 2自由度PID制御 コントローラは、コマンドの追跡方法とノイズの除去方法を分離します。比例経路と微分経路に重み値 $b$ と $c$ を使用します。$$u(t) = P \cdot (b \cdot r – y) + I \cdot \int (r – y) dt + D \cdot \frac{d}{dt}(c \cdot r – y)$$ これにより、設定値を素早く変更したときの「微分キック」が抑制され、アクチュエータの動きがスムーズになります。

2. 不連続性サブライブラリ

物理的なハードウェアは決して完璧ではありません。私たちはこれらのブロックを使用して、現実世界の限界をモデル化します。

  • バックラッシュブロック: ギアの遊びをシミュレートします。デッドバンド幅と開始出力を設定できます。
  • デッドゾーンブロック: 入力値が上限値または下限値を超えるまで、ブロックの出力はゼロのままです。
  • 飽和ブロック: 出力を制限します。上限と下限を設定します。常にオンにしてください。 ゼロ交差検出 つまり、ソルバーは信号が制限値に達する正確な瞬間を見つけ出す。

3. 個別サブライブラリ

これらのブロックは、デジタルプロセッサやECUをモデル化するために使用します。

  • ユニット遅延ブロック($z^{-1}$): 入力値を1クロックステップの間保持し、その後出力します。
  • ゼロオーダーホールド(ZOH)ブロック: 連続信号をサンプリングし、設定されたサンプリング時間($T_s$)の間、その値を一定に保ちます。

ヒント:ソルバー設定がシミュレーションに与える影響

を使用するとどうなるかを見てみましょう ユニット遅延 ブロック:

  • 可変ステップソルバー(ode45など)を使用する場合: ソルバーのステップサイズは動的に変化する可能性があります。これにより、離散ブロック内で補間誤差が発生し、本来連続的ではない信号が連続的に見えてしまうことがあります。
  • 固定ステップソルバーを使用する場合: ソルバーを次のように設定します 固定ステップ サンプル時間を1.0秒に設定することで、すべてが同期して動作します。ユニット遅延ブロックは値を正確に1秒間保持し、スコープ上にきれいで離散的な階段状のグラフを作成します。

2.2 LTIシステムブロックの使用

MATLABワークスペースに線形モデルがある場合は、以下の方法でSimulinkに取り込むことができます。 LTIシステムブロック (制御システムツールボックス内にあります。)

  • 設定方法:
    1. 直接入力してください: ブロックパラメータボックスにTFコマンドを直接入力します(例: tf([1 1]、[1 5 1]))
    2. ワークスペース変数を使用する: まず、MATLABスクリプトでモデルを定義します。G1 = tf([1 1]、[1 5 1]); 次に、 G1 LTIブロックに組み込みます。Dowway Vehicleでは、このスクリプト方式を使用してパラメータを整理しています。LTIモデルが状態空間形式である場合にのみ、開始状態を設定できることに注意してください。

パート3:ステップバイステップのプロジェクト

それでは、実用的なシステムを2つ構築し、運用してみましょう。

プロジェクト1:伝達関数安定性解析

連続伝達関数を解析し、ステップ入力および周波数掃引に対する挙動を調べます。

当社のプラントモデル:

$$G(s) = \frac{25(s + 17)}{(s + 25)(s + 37)} = \frac{25s + 425}{s^2 + 62s + 925}$$

MATLABスクリプト:

このスクリプトを実行すると、ボード線図、ナイキスト線図、ニコルス線図が生成されます。

% Project 1: Frequency Response Analysis
clear; clc; close all;

% Define numerator and denominator
num = 25 * [1 17];                 % 25s + 425
den = conv([1 25], [1 37]);        % s^2 + 62s + 925

% Create Transfer Function
G = tf(num, den);
disp('System Transfer Function:');
printsys(num, den);

% Draw Bode Plot
figure(1)
bode(G)
grid on;
title('System Bode Diagram');

% Draw Nyquist Plot
figure(2)
nyquist(G)
axis equal % Keeps the scale of real and imaginary axes equal
grid on;
title('System Nyquist Plot');

% Draw Nichols Plot
figure(3)
nichols(G)
grid on;
title('System Nichols Chart');
  1. Simulinkで空のモデルを開きます。
  2. 追加する ステップ ブロックから 情報源 ライブラリ。ステップ時間を$0$に設定します。
  3. 追加する 送金FCN ブロックから 連続 ライブラリ。「分子係数」を [25 425] そして「分母係数」 [1 62 925]
  4. 追加する 範囲 ブロックから シンク 図書館。
  5. ブロックをつなげる: [Step]->[Transfer Fcn]->[Scope]
  6. シミュレーションを1.0秒間実行してください。オシロスコープには、振動のない安定した滑らかなステップ応答が表示されます。

プロジェクト2:倒立振子のLQR状態フィードバック制御

倒立振子は不安定で制御が難しい。本稿では、その挙動を線形化し、特性を検証し、最適な設計を行う。 線形二次レギュレータ(LQR) バランスを保つため。

私たちの州ベクトル:

$$x = \begin{bmatrix} x & \theta & \dot{x} & \dot{\theta} \end{bmatrix}^T$$

ここで、$x$は台車の位置、$\theta$は振り子の角度、$\dot{x}$は台車の速度、$\dot{\theta}$は角速度です。

線形化行列(開ループ不安定):

$$A = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 0 & -0.524 & -3.129 & 0 \\ 0 & -15.332 & 9.012 & -0.562 \end{bmatrix}, \quad B = \begin{bmatrix} 0 \\ 0 \\ 81.23 \\ -120.1 \end{bmatrix}$$$$C = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 1 & 0 & $$ \begin{bmatrix} 0 \\ 0 \end{bmatrix}$$

これは MIMO システム。入力は力($u$)です。出力は台車の位置($x$)と振り子の傾斜角($\theta$)です。

MATLAB最適化スクリプト:

このコードを実行して、LQRフィードバック行列$K$を計算してください。

% Project 2: Inverted Pendulum LQR Design
clear; clc; close all;

% System Matrices
A = [0,       0,      1,      0;
     0,       0,      0,      1;
     0,  -0.524, -3.129,      0;
     0, -15.332,  9.012, -0.562];

B = [0; 
     0; 
 81.23; 
-120.1];

C = [1, 0, 0, 0;    % Output 1: Cart Position
     0, 1, 0, 0];   % Output 2: Pendulum Angle

D = [0; 
     0];

% Create state-space model
sys = ss(A, B, C, D);

% Calculate eigenvalues
E = eig(A);
fprintf('Open-Loop Eigenvalues:\n');
disp(E);

% Test Controllability
OC = ctrb(A, B);
N_C = rank(OC);
fprintf('Controllability Rank: %d\n', N_C);

% Test Observability
OB = obsv(A, C);
N_B = rank(OB);
fprintf('Observability Rank: %d\n', N_B);

% LQR Tuning Matrices
Q = diag([10, 0.01, 0.01, 0.1]); 
R = 45;
N = 0; 

% Compute feedback gain matrix K
[K, S, e] = lqr(sys, Q, R, N);
fprintf('Calculated LQR Feedback Gain Matrix K:\n');
disp(K);

% Create Closed-Loop System (u = -Kx)
sysT = ss(A - B*K, B, C, D);

% Plot closed-loop step response
figure(1)
step(sysT)
grid on;
legend('Cart Position (x)', 'Pendulum Angle (theta)');
title('Closed-Loop LQR Step Response');

開ループ固有値は、元のシステムが不安定であることを示しています。可制御性行列と可観測性行列はどちらもランクが4です。これは、システムが完全に可制御かつ可観測であることを意味し、LQRコントローラを構築できます。

Simulinkでこのモデルを実行するには:

  1. 追加する 状態空間 連続ライブラリからのブロック。
  2. パラメータを設定する:
    • セット AA - B*K
    • セット BB
    • セット CC
    • セット DD. (これらの変数がワークスペースで有効になるように、まずMATLABスクリプトを実行してください。)
  3. 接続する ステップ 入力ポートへのブロック。
  4. 接続する 範囲 出力ポートへのブロック。
  5. シミュレーションを実行してください。スコープには2本の線が表示されます。台車の位置はスムーズに安定し、振り子の角度は短い制御された揺れの後、0度(垂直)に戻ります。

ジョニーの安全制御システムチェックリスト

最後に、Simulinkで制御ループを構築するための私個人のチェックリストを以下に示します。

  • まず制御可能性をテストする: 複雑なプラントのゲイン調整に何時間も費やす前に、 ctrb MATLABの場合。ランクが低すぎると、どのアルゴリズムを使ってもシステムを安定させることはできません。
  • 適切なソルバーを選択してください: 非線形ブロックを持つシステムの場合 飽和 または 反発、 オンにする ゼロ交差検出. 可変ステップソルバーを使用する(例: ode45)は迅速なテストに使用しますが、実際のハードウェア向けにコードをコンパイルする準備ができたら、固定ステップソルバーに切り替えてください。
  • 変数を一箇所にまとめておく: Simulinkのパラメータボックスに生の数値を直接入力しないでください。すべての情報は中央のMATLAB初期化スクリプトに記述してください。こうすることで、作業が整理され、見やすくなります。

チューニングの問題やソルバーのエラーが発生した場合は、お知らせください。下のコメント欄にご記入いただければ、一緒に話し合いましょう!

コメントする

メールアドレスが公開されることはありません。 が付いている欄は必須項目です

Need a Quote or Have Questions?

Please fill out the form below, our engineers will contact you within 24 hours.