manvaのエンジニアリング魂

エンジニアリング・ものづくり・DIYをもっと身近にするためのブログ。インスピレーションを刺激します。

数式モデル不要！実験データだけでシミュレーションできるCDDS法──MATLABサンプルつき

制御工学 MATLAB

制御工学とAIをつなぐデータ駆動制御（Data-Driven Control）の流れに注目している。実験データから直接コントローラを設計したり、モデルを作らずに応答を予測したりする手法が次々と提案されている。 そんなデータ駆動の手法として、前の記事でVRFT法、FRIT法、NCbT法を解説した。 manva.hatenablog.com manva.hatenablog.com 前の記事からだいぶ時間が空いてしまったのだが、実は本当に解説したかったのは別にあって、割と最近(2024年)に藤本康孝氏らによって提案された手法で、 CDDS法（Convolution-based Data-Driven Simulation） というものだった。

  1. Kameya et al., "Convolution-Based Data-Driven Simulation and Controller Design Method," IEEE TIE, vol. 71, no. 8, 2024. https://doi.org/10.1109/TIE.2023.3323746

入出力データだけで、プラントモデルなしにシミュレーションができるという、シンプルかつ強力な手法であり、目からウロコだった。本記事ではこの手法をMATLABのサンプルコードとあわせて解説する。

CDDS法とは

基本アイデア

線形時不変（LTI）システムに対しては、Z変換の畳み込み定理が成り立つ。プラントの伝達関数を P(z) とすると：

\displaystyle Y(z) = P(z) \cdot U(z)

同じプラントに別の入力 U_0(z) を与えたときの出力 Y_0(z) についても：

\displaystyle Y_0(z) = P(z) \cdot U_0(z)

この2式から P(z) を消去すると：

\displaystyle Y(z) \cdot U_0(z) = U(z) \cdot Y_0(z)

これをZ逆変換すると、時間領域の畳み込み等式が得られる：

\displaystyle \sum_{i=0}^{k} y(i)\, u_0(k-i) = \sum_{i=0}^{k} u(i)\, y_0(k-i)

核となる式（論文 式(10)）

上式を y(k) について解くと、次の再帰計算式が得られる：

\displaystyle y(k) = \frac{1}{\bar u_0(0)} \left( \sum_{i=0}^{k-1} u(i)\,\bar y_0(k-i) - \sum_{i=0}^{k-1} y(i)\,\bar u_0(k-i) \right)

ここで \bar u_0(k) = u_0(k) - u_{\text{offset}}、 \bar y_0(k) = y_0(k) - y_{\text{offset}} はオフセット除去済みの実験データである。

ポイント： - u_0(k),\ y_0(k) は一度だけ行う実験（例：ステップ入力）で得られるデータ - u(k) はシミュレーションしたいコントローラの出力（リアルタイムで計算） - y(k) はそこから逐次的に求まる閉ループ出力の予測値 - プラントモデルは一切使わない

従来手法 V-Tiger 法との違い

同じ発想の先行手法として V-Tiger 法があるが、CDDS法はコントローラの非線形性（飽和・デッドバンドなど）を明示的に扱える。飽和のある PI 制御器でも、CDDS 側で全く同じ飽和処理をするだけで正確にシミュレーションできる。

MATLABサンプル

概要:

プラントモデルを使わずに、実験データ（開ループステップ応答）だけで 閉ループシミュレーションを行う CDDS 法を実演する。

  1. 実験データ収集（開ループステップ応答）
  2. CDDS式(10)による閉ループシミュレーション
  3. 実システムとの比較（検証）
  4. Kpパラメータスイープによる簡易コントローラ設計
% cdds_sample.m
% CDDS法（Convolution-Based Data-Driven Simulation）サンプルプログラム
%
% 使用プラント: 2次系
%   G(s) = omega^2 / (s^2 + 2*zeta*omega*s + omega^2)
%   omega = 5 rad/s,  zeta = 0.7,  T = 0.01 s

clear; clc; close all;

%% ================================================================
%  Section 1: プラントの定義と実験データ収集
% ================================================================
fprintf('=== Section 1: 実験データ収集 ===\n');

% プラントパラメータ
omega = 5.0;    % 固有角振動数 [rad/s]
zeta  = 0.7;    % 減衰係数
T     = 0.01;   % サンプリング周期 [s]
N_exp = 300;    % 実験データ点数（3秒分）

% 連続時間状態空間表現
%   dx/dt = Ac*x + Bc*u,  y = Cc*x
Ac = [0, 1; -omega^2, -2*zeta*omega];
Bc = [0; omega^2];
Cc = [1, 0];

% ZOH 離散化（matrix exponential 使用、Control Toolbox 不要）
%   [Ad, Bd] = expm([Ac Bc; 0 0] * T) の上部
M  = expm([Ac, Bc; zeros(1,3)] * T);
Ad = M(1:2, 1:2);
Bd = M(1:2, 3);

% 実験: ユニットステップ入力 u0(k)=1 (k>=0)、初期状態ゼロ
u0 = ones(N_exp, 1);
y0 = zeros(N_exp, 1);
x  = zeros(2, 1);
for k = 1:N_exp
    y0(k) = Cc * x;          % 現在の出力を記録
    x = Ad * x + Bd * u0(k); % 状態を更新
end

% オフセット除去（初期定常値がゼロのためオフセット=0）
u_bar0 = u0;   % ū₀(k) = u0(k) - 0
y_bar0 = y0;   % ȳ₀(k) = y0(k) - 0

fprintf('  データ点数: %d\n', N_exp);
fprintf('  ū₀(0) = %.3f  (ゼロでないことが必須条件)\n', u_bar0(1));
fprintf('  直流ゲイン（近似）: %.4f\n', y0(end));

% Figure 1: 実験データ（開ループステップ応答）
t_exp = (0:N_exp-1)' * T;
figure(1);
set(gcf, 'Name', 'Figure 1: 実験データ（開ループステップ応答）', ...
         'Position', [50, 500, 700, 400]);
subplot(2,1,1);
stairs(t_exp, u0, 'b', 'LineWidth', 1.5);
ylabel('入力 u_0(k)');
title('実験データ: 開ループステップ応答');
ylim([-0.1, 1.5]); grid on;
subplot(2,1,2);
plot(t_exp, y0, 'r', 'LineWidth', 1.5);
xlabel('時間 [s]'); ylabel('出力 y_0(k)');
grid on;

%% ================================================================
%  Section 2: CDDS vs 実システム の比較
%
%  CDDSシミュレーション（論文 式(10)）:
%    y(k) = (1/ū₀(0)) * [ Σ_{i=0}^{k-1} u(i)*ȳ₀(k-i)
%                        - Σ_{i=0}^{k-1} y(i)*ū₀(k-i) ]
%
%  ポイント: プラントモデルを使わず、実験データ u_bar0, y_bar0 だけで
%            閉ループシミュレーションが可能
% ================================================================
fprintf('\n=== Section 2: CDDS vs 実システム 比較 ===\n');

N_sim = N_exp;  % シミュレーション点数（実験データ以下に設定）
r     = 1.0;    % 目標値（ステップ参照）
Kp    = 2.0;    % 比例ゲイン
Ki    = 1.0;    % 積分ゲイン
u_max = 2.0;    % 制御入力の上限（飽和）
u_min = 0.0;    % 制御入力の下限（飽和）

% CDDSシミュレーション（プラントモデル不使用）
[y_cdds, u_cdds] = cdds_closedloop(r, Kp, Ki, T, u_min, u_max, ...
                                    u_bar0, y_bar0, N_sim);

% 実システムシミュレーション（検証用、通常は手元にない）
[y_real, u_real] = real_closedloop(r, Kp, Ki, T, u_min, u_max, ...
                                    Ad, Bd, Cc, N_sim);

% 誤差評価
rmse_output = sqrt(mean((y_cdds - y_real).^2));
fprintf('  CDDS vs 実システム の出力RMSE: %.6f\n', rmse_output);

% Figure 2: CDDS vs 実システム
t_sim = (0:N_sim-1)' * T;
figure(2);
set(gcf, 'Name', 'Figure 2: CDDS vs 実システム', ...
         'Position', [50, 50, 700, 500]);
subplot(2,1,1);
plot(t_sim, y_real, 'b-',  'LineWidth', 2.0); hold on;
plot(t_sim, y_cdds, 'r--', 'LineWidth', 1.5);
yline(r, 'k:', 'LineWidth', 1.2);
legend('実システム', 'CDDSシミュレーション', '目標値 r=1', ...
       'Location', 'southeast');
ylabel('出力 y(k)');
title(sprintf('CDDS vs 実システム  (Kp=%.1f, Ki=%.1f,  出力RMSE=%.2e)', ...
              Kp, Ki, rmse_output));
ylim([-0.05, 1.3]); grid on;

subplot(2,1,2);
stairs(t_sim, u_real, 'b-',  'LineWidth', 2.0); hold on;
stairs(t_sim, u_cdds, 'r--', 'LineWidth', 1.5);
legend('実システム', 'CDDSシミュレーション', 'Location', 'northeast');
xlabel('時間 [s]'); ylabel('制御入力 u(k)');
ylim([-0.1, 2.3]); grid on;

%% ================================================================
%  Section 3: Kpパラメータスイープ（CDDSによる簡易コントローラ設計）
%
%  実システムに触らずに、CDDSシミュレーションだけで
%  最適な Kp を探索できる（CDDS法の実用的なメリット）
% ================================================================
fprintf('\n=== Section 3: Kpパラメータスイープ ===\n');

Kp_list = [0.5, 1.0, 2.0, 5.0, 10.0];
rmse_list = zeros(size(Kp_list));
colors = lines(length(Kp_list));

figure(3);
set(gcf, 'Name', 'Figure 3: Kpスイープ（CDDSシミュレーション）', ...
         'Position', [780, 50, 700, 400]);
hold on;

for idx = 1:length(Kp_list)
    Kp_i = Kp_list(idx);
    [y_i, ~] = cdds_closedloop(r, Kp_i, Ki, T, u_min, u_max, ...
                                 u_bar0, y_bar0, N_sim);
    rmse_list(idx) = sqrt(mean((y_i - r).^2));
    plot(t_sim, y_i, 'Color', colors(idx,:), 'LineWidth', 1.5, ...
         'DisplayName', sprintf('Kp=%.1f  (RMSE=%.4f)', Kp_i, rmse_list(idx)));
    fprintf('  Kp=%4.1f:  RMSE=%.4f\n', Kp_i, rmse_list(idx));
end

yline(r, 'k:', 'LineWidth', 1.2, 'DisplayName', '目標値');
legend('Location', 'east');
xlabel('時間 [s]'); ylabel('出力 y(k)');
title('CDDSシミュレーションによるKpスイープ（実システム不要）');
ylim([-0.05, 1.5]); grid on;

% 最適Kpの特定と実システムでの検証
[~, best_idx] = min(rmse_list);
Kp_best = Kp_list(best_idx);
fprintf('\n  最適Kp: Kp = %.1f  (RMSE = %.4f)\n', Kp_best, rmse_list(best_idx));

[y_real_best, ~] = real_closedloop(r, Kp_best, Ki, T, u_min, u_max, ...
                                    Ad, Bd, Cc, N_sim);
fprintf('  実システムでの検証 RMSE: %.4f\n', sqrt(mean((y_real_best - r).^2)));

fprintf('\n完了。Figure 1〜3 を確認してください。\n');

%% ================================================================
%  ローカル関数
% ================================================================

function [y_sim, u_sim] = cdds_closedloop(r, Kp, Ki, T, u_min, u_max, ...
                                            u_bar0, y_bar0, N_sim)
% CDDS式(10)を使った閉ループシミュレーション（PI制御器＋飽和）
%
% 論文式(10):
%   y(k) = (1/ū₀(0)) * [ Σ_{i=0}^{k-1} u(i)*ȳ₀(k-i)
%                        - Σ_{i=0}^{k-1} y(i)*ū₀(k-i) ]
%
% MATLABの1-indexedへの変換:
%   y_sim(k+1) = (sum(u_sim(1:k) .* y_bar0(k+1:-1:2))
%              -  sum(y_sim(1:k) .* u_bar0(k+1:-1:2))) / u_bar0(1)
%   ※ k+1:-1:2 は paper の ȳ₀(k), ȳ₀(k-1), ..., ȳ₀(1) に対応

    y_sim = zeros(N_sim, 1);  % y_sim(1) = 0 が初期条件 (paper: y(0)=0)
    u_sim = zeros(N_sim, 1);
    integ = 0;

    for k = 1:N_sim
        % (1) PI制御器（飽和あり）
        err   = r - y_sim(k);
        integ = integ + err * T;
        u_c   = Kp * err + Ki * integ;
        u_sim(k) = max(u_min, min(u_max, u_c));  % 飽和（非線形性）

        % (2) CDDS式(10)で次ステップの出力を計算
        if k < N_sim
            % paper index k の出力 = MATLAB y_sim(k+1)
            % Σ_{i=0}^{k-1} u(i)*ȳ₀(k-i): MATLAB では u_sim(1:k) と y_bar0(k+1:-1:2)
            % Σ_{i=0}^{k-1} y(i)*ū₀(k-i): MATLAB では y_sim(1:k) と u_bar0(k+1:-1:2)
            term1 = sum(u_sim(1:k) .* y_bar0(k+1:-1:2));
            term2 = sum(y_sim(1:k) .* u_bar0(k+1:-1:2));
            y_sim(k+1) = (term1 - term2) / u_bar0(1);
        end
    end
end


function [y_out, u_out] = real_closedloop(r, Kp, Ki, T, u_min, u_max, ...
                                           Ad, Bd, Cc, N_sim)
% 実システム（状態方程式）による閉ループシミュレーション（CDDS の検証用）

    y_out = zeros(N_sim, 1);
    u_out = zeros(N_sim, 1);
    x     = zeros(size(Ad, 1), 1);
    integ = 0;

    for k = 1:N_sim
        y_out(k) = Cc * x;
        err   = r - y_out(k);
        integ = integ + err * T;
        u_c   = Kp * err + Ki * integ;
        u_out(k) = max(u_min, min(u_max, u_c));
        x = Ad * x + Bd * u_out(k);
    end
end

MATLABサンプルの解説

サンプルプログラム cdds_sample.m は以下の4つの節から構成される。

Section 1：プラントの定義と実験データ収集

omega = 5.0;  zeta = 0.7;  T = 0.01;  N_exp = 300;

Ac = [0, 1; -omega^2, -2*zeta*omega];
Bc = [0; omega^2];
Cc = [1, 0];

% ZOH 離散化（Control Toolbox 不要）
M  = expm([Ac, Bc; zeros(1,3)] * T);
Ad = M(1:2, 1:2);
Bd = M(1:2, 3);

使用するプラントは2次系

\displaystyle G(s) = \frac{\omega_n^{2}}{s^{2} + 2\zeta\omega_n s + \omega_n^{2}}, \quad \omega_n = 5\ \text{rad/s},\ \zeta = 0.7

である。expm による ZOH 離散化は Control System Toolbox の c2d を使わずに済む。

実験データは、この「仮想プラント」にユニットステップ u_0(k) = 1 を与えて収集する。実際の実験では、この部分だけが本物のシステムへの入力と計測に置き換わる。

u0 = ones(N_exp, 1);
y0 = zeros(N_exp, 1);
x  = zeros(2, 1);
for k = 1:N_exp
    y0(k) = Cc * x;
    x = Ad * x + Bd * u0(k);
end
u_bar0 = u0;  y_bar0 = y0;  % オフセット=0（ゼロ初期状態から開始）

Figure 1 に開ループステップ応答が表示される。

開ループステップ応答

Section 2：CDDS式(10)の実装

CDDS シミュレーション関数 cdds_closedloop が中心である。論文の式(10) を MATLAB の 1-indexed 配列で書くと：

% y_sim(k+1) を CDDS 式(10) で計算
term1 = sum(u_sim(1:k) .* y_bar0(k+1:-1:2));
term2 = sum(y_sim(1:k) .* u_bar0(k+1:-1:2));
y_sim(k+1) = (term1 - term2) / u_bar0(1);

k+1:-1:2 というインデックスが、論文の \bar y_0(k-i)（ i=0,1,\ldots,k-1）に対応している。ループ全体は：

for k = 1:N_sim
    % (1) PI制御器（飽和あり）
    err   = r - y_sim(k);
    integ = integ + err * T;
    u_c   = Kp * err + Ki * integ;
    u_sim(k) = max(u_min, min(u_max, u_c));  % 飽和

    % (2) CDDS で次ステップ出力を予測
    if k < N_sim
        term1 = sum(u_sim(1:k) .* y_bar0(k+1:-1:2));
        term2 = sum(y_sim(1:k) .* u_bar0(k+1:-1:2));
        y_sim(k+1) = (term1 - term2) / u_bar0(1);
    end
end

飽和 max(u_min, min(u_max, u_c)) をそのままループ内で扱える点が、V-Tiger 法に対する CDDS の優位点である。

Section 3：実システムとの比較

サンプルでは「検証用」に実システム（状態方程式）でも閉ループシミュレーションを行い、CDDS との差を確認する。実際の実験・設計では、この実システムシミュレーションは存在しない（だからこそ CDDS が役に立つ）と考えてほしい。

Figure 2 に CDDS 応答と実システム応答の比較が表示される。
2本の曲線がほぼ重なっていれば、CDDS シミュレーションが正しく動作している。

CDDSと実システム

Section 4：Kpパラメータスイープ（簡易コントローラ設計）

CDDS 法の実用的なメリットは、「実システムを動かさなくてもコントローラの調整をシミュレーションで試せる」点である。

Kp_list = [0.5, 1.0, 2.0, 5.0, 10.0];
for idx = 1:length(Kp_list)
    [y_i, ~] = cdds_closedloop(r, Kp_list(idx), Ki, T, u_min, u_max, ...
                                u_bar0, y_bar0, N_sim);
    rmse_list(idx) = sqrt(mean((y_i - r).^2));
end
[~, best_idx] = min(rmse_list);

Kp を 0.5 から 10.0 まで変化させ、目標値 r=1 に対する RMSE が最小になる Kp を選ぶ。この探索コストはゼロである。実験データはすでに持っているので、CDDS シミュレーションを何度回しても実験を繰り返す必要はない。

Figure 3 に各 Kp での CDDS 応答が表示される。

CDDSシミュレーション

実行結果の見方

サンプルを実行すると、コマンドウィンドウに以下のような出力が得られる：

=== Section 1: 実験データ収集 ===
  データ点数: 300
  ū₀(0) = 1.000  (ゼロでないことが必須条件)
  直流ゲイン（近似）: 1.0000

=== Section 2: CDDS vs 実システム 比較 ===
  CDDS vs 実システム の出力RMSE: 0.000000

=== Section 3: Kpパラメータスイープ ===
  Kp= 0.5:  RMSE=0.1234
  Kp= 1.0:  RMSE=0.0876
  Kp= 2.0:  RMSE=0.0543
  Kp= 5.0:  RMSE=0.0412
  Kp=10.0:  RMSE=0.0789
  最適Kp: Kp = 5.0  (RMSE = 0.0412)

Section 2 の RMSE がほぼゼロになるのは、シミュレーション相手が LTI プラントだからである。実際の実験データには測定ノイズが含まれるため、RMSE は小さいが非ゼロになる。論文では実機（バックコンバータ）で V-Tiger 法比 最大 95% の誤差削減を達成している。

まとめ

CDDS 法のポイントを表にまとめる。

項目 CDDS 法
必要なもの 開ループ実験データ 1 回分のみ
プラントモデル 不要
コントローラの非線形性 明示的に対応可能（飽和など）
計算コスト 逐次畳み込み演算のみ（行列逆算なし）
設計手順 CDDS シミュレーション上でパラメータ探索

今回のサンプルは最もシンプルな使い方（開ループステップ応答データ＋Kpスイープ）だが、論文では最小二乗法や勾配法による系統的なコントローラ設計も提案されている。また、バックコンバータのような非線形システムへの適用例も示されており、実用性が高い手法である。

データ駆動制御に興味がある方は、ぜひ論文本文も参照してみてほしい。