【実用ガイド】制御工学でよく使うMATLABコマンド一覧

制御工学を学んでいると、システムの解析や設計に様々なコマンドを使います。この記事では、MATLABの Control System Toolbox でよく使われるコマンドとその使い方を実例付きで紹介します。特に伝達関数の解析に焦点を当てています。

目次

1. 伝達関数の定義と基本操作
2. 極とゼロ点の解析
3. 時間応答の解析
4. 周波数応答の解析
5. システムの変換
6. 制御系の設計と解析
7. サンプルスクリプト

伝達関数の定義と基本操作

伝達関数の作成

% 方法1: 分子・分母の多項式係数を使って定義
num = [1 2];         % 分子の係数 (s+2)
den = [1 3 5];       % 分母の係数 (s^2+3s+5)
G1 = tf(num, den);   % G(s) = (s+2)/(s^2+3s+5)

% 方法2: sをシンボリックに使って定義
s = tf('s');
G2 = (s+2)/(s^2+3*s+5);

% 方法3: zero-pole-gain形式で定義
z = [-2];            % ゼロ点
p = [-1.5+1.66i, -1.5-1.66i]; % 極
k = 1;               % ゲイン
G3 = zpk(z, p, k);   % G(s) = k(s-z1)(s-z2).../(s-p1)(s-p2)...

% 作成した伝達関数の表示
G1

伝達関数の演算

% 伝達関数の演算
G4 = G1 + G2;        % 並列接続 (加算)
G5 = G1 * G2;        % 直列接続 (乗算)
G6 = feedback(G1, G2); % フィードバック接続 G1/(1+G1*G2)

% 伝達関数の特性表示
tf(G1)               % 伝達関数形式で表示
zpk(G1)              % ゼロ・極・ゲイン形式で表示

極とゼロ点の解析

極とゼロ点は、システムの動特性を理解する上で非常に重要です。

極とゼロ点の取得

% 極とゼロ点の取得
p = pole(G1)         % G1の極を取得
z = zero(G1)         % G1のゼロ点を取得

% 極とゼロ点のプロット
figure;
pzmap(G1);           % 極・ゼロマップをプロット
grid on;
title('極・ゼロマップ');

% 複数のシステムの極・ゼロマップを比較
figure;
pzmap(G1, 'r', G2, 'b');
grid on;
legend('System G1', 'System G2');

極とゼロ点に関する追加情報

% 全ての極の実部が負か確認（安定性チェック）
stable = isstable(G1)

% 極の実部と虚部を個別に取得
poles = pole(G1);
real_parts = real(poles)
imag_parts = imag(poles)

% 減衰比と固有角周波数を計算
[wn, zeta] = damp(G1)

% ゼロと極の個数の確認
order = order(G1)           % システムの次数（極の数）
n_zeros = length(zero(G1))  % ゼロの数

時間応答の解析

システムの時間応答は、ステップ入力やインパルス入力などに対する出力の時間変化を示します。

ステップ応答とインパルス応答

% ステップ応答のプロット
figure;
step(G1);
grid on;
title('ステップ応答');

% 複数のシステムのステップ応答比較
figure;
step(G1, 'r', G2, 'b');
grid on;
legend('System G1', 'System G2');

% インパルス応答のプロット
figure;
impulse(G1);
grid on;
title('インパルス応答');

% カスタム時間範囲でのプロット
t = 0:0.01:10;      % 0から10秒まで0.01秒刻み
figure;
step(G1, t);
grid on;

時間応答の特性解析

% ステップ応答の特性（立ち上がり時間、オーバーシュートなど）
info = stepinfo(G1)

% 自由応答（初期状態からの応答）
[y, t] = initial(G1, [1; 0], 0:0.01:5);  % 初期状態[1; 0]からの応答
figure;
plot(t, y);
grid on;
title('自由応答');

任意入力に対する応答

% 任意入力に対する応答
t = 0:0.01:10;
u = sin(t);         % 正弦波入力
[y, t] = lsim(G1, u, t);  % 線形シミュレーション
figure;
subplot(2,1,1);
plot(t, u);
grid on;
title('入力信号');
subplot(2,1,2);
plot(t, y);
grid on;
title('出力応答');

周波数応答の解析

周波数応答は、システムが様々な周波数の正弦波入力に対してどのように応答するかを示します。

ボード線図

% ボード線図（ゲイン特性と位相特性）
figure;
bode(G1);
grid on;

% 周波数範囲を指定したボード線図
w = logspace(-2, 2, 100);  % 10^-2から10^2までの対数スケール
figure;
bode(G1, w);
grid on;

% 複数システムのボード線図比較
figure;
bode(G1, 'r', G2, 'b');
grid on;
legend('System G1', 'System G2');

ナイキスト線図とニコルス線図

% ナイキスト線図
figure;
nyquist(G1);
grid on;
title('ナイキスト線図');

% ニコルス線図
figure;
nichols(G1);
grid on;
title('ニコルス線図');

周波数応答データの取得

% 特定の周波数での応答を取得
w = [0.1, 1, 10];    % 周波数点
[mag, phase] = bode(G1, w);

% マグニチュードをdBに変換
mag_db = 20*log10(mag(:,:))

% 位相を度数に変換（すでに度数で出力されるが、確認のため）
phase_deg = phase(:,:)

% 周波数応答データの詳細な取得
[mag, phase, wout] = bode(G1, w);  % woutは実際に使用された周波数

% ゲイン余裕と位相余裕
[Gm, Pm, Wcg, Wcp] = margin(G1);
disp(['ゲイン余裕: ' num2str(Gm) ' dB']);
disp(['位相余裕: ' num2str(Pm) ' deg']);
disp(['位相交差周波数: ' num2str(Wcg) ' rad/s']);
disp(['ゲイン交差周波数: ' num2str(Wcp) ' rad/s']);

システムの変換

制御系の解析では、伝達関数、状態空間表現、零点-極-ゲイン表現など異なる表現を使い分けることがあります。

表現形式の変換

% 伝達関数から状態空間表現へ
[A, B, C, D] = tf2ss(num, den);
sys_ss = ss(A, B, C, D);

% 状態空間表現から伝達関数へ
[num, den] = ss2tf(A, B, C, D);
sys_tf = tf(num, den);

% 伝達関数からzero-pole-gain表現へ
sys_zpk = zpk(G1);

% zero-pole-gainから伝達関数へ
sys_tf2 = tf(sys_zpk);

離散時間システムへの変換

% 連続時間システムから離散時間システムへ
Ts = 0.1;  % サンプリング時間
G_d = c2d(G1, Ts, 'zoh');  % ゼロ次ホールド法
G_d2 = c2d(G1, Ts, 'tustin');  % 双一次変換法

% 離散時間システムのステップ応答
figure;
step(G_d);
grid on;
title('離散時間システムのステップ応答');

制御系の設計と解析

PID制御器の設計

% PID制御器の設計
Kp = 1;    % 比例ゲイン
Ki = 0.5;  % 積分ゲイン
Kd = 0.2;  % 微分ゲイン

% PID制御器の伝達関数
C = pid(Kp, Ki, Kd);

% プラントと制御器を組み合わせたオープンループ伝達関数
G_ol = C * G1;

% クローズドループ伝達関数
G_cl = feedback(G_ol, 1);

% ステップ応答の確認
figure;
step(G_cl);
grid on;
title('PID制御系のステップ応答');

ルートローカス（根軌跡）

% 根軌跡のプロット
figure;
rlocus(G1);
grid on;
title('根軌跡');

% 特定のゲインでの極の位置を取得
K = 2.5;
poles = rlocus(G1, K);

LQR制御器設計（状態空間表現用）

% 状態空間表現でのシステム
A = [-3 1; -5 -2];
B = [0; 1];
C = [1 0];
D = 0;
sys = ss(A, B, C, D);

% LQR制御器設計
Q = eye(2);  % 状態重み行列
R = 1;       % 入力重み行列
K = lqr(sys, Q, R);

% 閉ループシステム
sys_cl = ss(A-B*K, B, C, D);

% ステップ応答
figure;
step(sys_cl);
grid on;
title('LQR制御系のステップ応答');

サンプルスクリプト

以下に、一連の解析を行うサンプルスクリプトを示します。このスクリプトは、二次系の伝達関数を定義し、極・ゼロ点解析、時間応答解析、周波数応答解析を行います。

% サンプルスクリプト：二次系の解析
clear;
clc;
close all;

% パラメータ定義
wn = 1;         % 固有角周波数
zeta_values = [0.2, 0.5, 0.7, 1.0, 1.2]; % 減衰比
colors = {'r', 'g', 'b', 'm', 'k'};      % プロット用の色

% 図の準備
figure('Position', [100, 100, 1200, 800]);

% 各減衰比に対して解析
for i = 1:length(zeta_values)
    zeta = zeta_values(i);
    
    % 伝達関数の定義
    s = tf('s');
    G = wn^2 / (s^2 + 2*zeta*wn*s + wn^2);
    
    % 極とゼロ点
    p = pole(G);
    z = zero(G);
    
    % ステップ応答
    subplot(2, 2, 1);
    step(G, 10);
    hold on;
    grid on;
    title('ステップ応答');
    xlabel('時間 [s]');
    ylabel('振幅');
    
    % 極・ゼロマップ
    subplot(2, 2, 2);
    pzmap(G);
    hold on;
    grid on;
    axis([-2 0.5 -2 2]);
    title('極・ゼロマップ');
    
    % ボード線図（ゲイン）
    subplot(2, 2, 3);
    w = logspace(-1, 1, 100);
    [mag, phase] = bode(G, w);
    mag_db = 20*log10(squeeze(mag));
    semilogx(w, mag_db, colors{i});
    hold on;
    grid on;
    title('周波数応答（ゲイン）');
    xlabel('周波数 [rad/s]');
    ylabel('ゲイン [dB]');
    
    % ボード線図（位相）
    subplot(2, 2, 4);
    phase_deg = squeeze(phase);
    semilogx(w, phase_deg, colors{i});
    hold on;
    grid on;
    title('周波数応答（位相）');
    xlabel('周波数 [rad/s]');
    ylabel('位相 [deg]');
    
    % ステップ応答の特性
    info = stepinfo(G);
    fprintf('減衰比 ζ = %.1f の場合:\n', zeta);
    fprintf('  極: %.3f + %.3fi, %.3f + %.3fi\n', real(p(1)), imag(p(1)), real(p(2)), imag(p(2)));
    fprintf('  立ち上がり時間: %.3f 秒\n', info.RiseTime);
    fprintf('  整定時間: %.3f 秒\n', info.SettlingTime);
    fprintf('  オーバーシュート: %.2f %%\n\n', info.Overshoot);
end

% 凡例の追加
subplot(2, 2, 1);
legend('\zeta = 0.2', '\zeta = 0.5', '\zeta = 0.7', '\zeta = 1.0', '\zeta = 1.2');
subplot(2, 2, 3);
legend('\zeta = 0.2', '\zeta = 0.5', '\zeta = 0.7', '\zeta = 1.0', '\zeta = 1.2');
subplot(2, 2, 4);
legend('\zeta = 0.2', '\zeta = 0.5', '\zeta = 0.7', '\zeta = 1.0', '\zeta = 1.2');

このスクリプトを実行すると、異なる減衰比を持つ二次系の特性を比較するグラフと、各システムの詳細情報（極の位置、立ち上がり時間、整定時間、オーバーシュートなど）がコマンドウィンドウに表示されます。

まとめ

この記事では、MATLABを使った制御系解析の基本的なコマンドとその使い方を紹介しました。これらのコマンドを使いこなすことで、制御系の特性を詳細に解析し、制御器の設計・評価を効率的に行うことができます。特に、極とゼロ点の解析、時間応答と周波数応答の解析は、制御工学において非常に重要なスキルです。

より高度な制御系設計（ロバスト制御、モデル予測制御など）や特殊なシステム（非線形系、時変系など）の解析には、追加のツールボックスやより専門的なコマンドが必要になることもありますが、この記事で紹介したコマンドは、制御工学を学ぶ上での基礎となるものです。

基本的にはMTLABが公式で出しているドキュメントを確認することを推奨します。

「MATLA グラフ出力」　「MATLAB 極」

などで検索すれば公式ドキュメントが見つかります。最近ではchatGPTなどに「matlabで○○したいから教えて」と言えばコードが返ってきたりもします。