周波数特性の測定データから二次系の伝達関数パラメータを定めるMATLABプログラム
伝達関数の分子が定数のみの場合とsの1次式の場合とは，"% Initial solution / 初期解" の所で x0 のどちらをコメントアウトするかによって変えられる。 $ \frac{b_{1}s+b_{0}}{s^{2}+a_{1}s+a_{0}} ($ b_{1}=0または$ b_{1}\neq0)
分母多項式から，固有角周波数$ \omega_{n} = \sqrt{a_0}と減衰定数$ \zeta = \frac{a_{1}}{2\omega_{n}}を計算する。 ゲイン特性と位相特性の誤差それぞれに異なる重み係数 c1，c2 をかけられる。より測定誤差の小さい方に大きな重みをかけるとよい。
周波数特性（ボード線図）の測定データは，プログラム内に書くようになっているが，csvやExcelファイルから読み込めるようにすると便利であると思われる。 code:MATLAB
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ボード線図測定データから2次系の伝達関数のパラメータを定めるプログラム
% Program for determining the parameters of the secondary system
% transfer function from Bode diagram measurement data
%
% 上智大学 電気電子工学実験III サーボモータ課題用
%
% Available for 2nd order denominator and constant and 1st order numerator
% 分母が2次，分子が定数のみまたは1次の伝達関数に対応
%
% Required toolboxes: optimization toolbox
%
% 2025-10-30 作成
% by M. Miyatake @Sophia University / 宮武 昌史 @上智大学
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
% Initialization / 初期化
%%%%%%%%%%%%%%%%%%%%%%
clear all;
%%%%%%%%%%%%%%%%%%%%%%
% Measured data / 測定データ (Hz, dB, deg)
%%%%%%%%%%%%%%%%%%%%%%
m_gain = 10.^(m_dB/20); % dB -> ratio
%%%%%%%%%%%%%%%%%%%%%%
% Initial solution / 初期解
%%%%%%%%%%%%%%%%%%%%%%
% Using it if the denominator is constant.
% 分子が定数項のみならこちらを使う
% Using it if the denominator contains the first order term of s.
% 分子がsの1次の項を含むならこちらを使う
%%%%%%%%%%%%%%%%%%%%%%
% Deriving parameters of the transfer function to minimize the sum of squared error
% 二乗和誤差を最小にする伝達関数パラメータを求める
%%%%%%%%%%%%%%%%%%%%%%
x, err, exitflag = fmincon(@(x)quad_err(x,m_freq,m_dB,m_gain,m_angle),x0,[],[],[],[],[],[], ... [],optimset('Algorithm', 'sqp', 'TolX', 1e-14, ...
'TolFun', 1e-14, 'MaxFunEvals', 200000, 'MaxIter', 100000 ));
%%%%%%%%%%%%%%%%%%%%%%
% Extracting the optimal parameters / 最適パラメータの抽出
%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%
% Drawing the bode diagram / ボード線図の描画
%%%%%%%%%%%%%%%%%%%%%%
num = b1 b0; % Numerator coefficients / 伝達関数の分子 den = a2 a1 a0; % Denominator coefficients / 伝達関数の分母 minp = floor(log10(min(m_freq))-0.01); % Lowest frequency / 周波数下限 (10^minp)
maxp = ceil(log10(max(m_freq))); % Highest frequency / 周波数上限 (10^maxp)
n = 200; % 周波数分割数
freq = logspace(minp, maxp, n); % Theoretical plot points / 理論プロット点
omega = 2*pi*freq;
sys = tf(num, den); % Transfer function / 伝達関数
gain = squeeze(gain);
gaindB = 20*log10(gain); % Gain in dB
phase = squeeze(phase); % Phase in deg
phase = phase - ceil(max(phase)/360)*360;
subplot(2,1,1); % Gain characteristics / ゲイン特性
semilogx(freq, gaindB, '-', m_freq, m_dB, '+');
xlabel('Frequency (Hz)');
ylabel('Gain (dB)');
grid( 'on');
legend('Calculated', 'Measured');
subplot(2,1,2); % Phase characteristics / 位相特性
semilogx(freq, phase, '-', m_freq, m_angle, '+');
xlabel('Frequency (Hz)');
ylabel('Phase (deg)');
grid('on');
legend('Calculated', 'Measured');
%%%%%%%%%%%%%%%%%%%%%%
% Display the transfer function parameters / 伝達関数パラメータ表示
%%%%%%%%%%%%%%%%%%%%%%
disp('Transfer Function Parameters:');
disp(b1: ', num2str(b1));
disp(b0: ', num2str(b0));
disp(a2: ', num2str(a2));
disp(a1: ', num2str(a1));
disp(a0: ', num2str(a0));
disp(' ');
%%%%%%%%%%%%%%%%%%%%%%
% Display the system characteristics / システム特性表示
%%%%%%%%%%%%%%%%%%%%%%
omegan = sqrt(a0); % Natural angular frequency / 固有角周波数
zeta = a1/2/omegan; % Damping factor / 減衰定数
disp('System Parameters:');
disp(Sum of squared error: ', num2str(err));
disp(Natural angular frequency: ', num2str(omegan),' rad/s'); disp(Damping factor: ', num2str(zeta));
disp(' ');
sys
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% quad_err
% Calculating the sum of squared error between theoretical and measured
% values in frequency characteristics
% 周波数特性の測定値と理論値の二乗誤差を計算
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function err = quad_err(x,m_freq,m_dB,m_gain,m_angle) %%%%%%%%%%%%%%%%%%%%%%
% Constants / 定数
%%%%%%%%%%%%%%%%%%%%%%
c1 = 1e-2; % Weighting coefficient for gain_dB / ゲイン特性の誤差に対する重み係数
c2 = 1e-4; % Weighting coefficient for angle / 位相特性の誤差に対する重み係数
%%%%%%%%%%%%%%%%%%%%%%
% Parameters / パラメータ
%%%%%%%%%%%%%%%%%%%%%%
% Transfer function
%%%%%%%%%%%%%%%%%%%%%%
% Calculation / 計算
%%%%%%%%%%%%%%%%%%%%%%
% Measured transfer function
m_trans = m_gain .* exp(1j * m_angle*pi/180);
% Calculated transfer function
c_trans = (b1*1j*2*pi*m_freq + b0) ./ (-a2*(2*pi*m_freq).^2 + a1*1j*2*pi*m_freq + a0);
c_gain = abs(c_trans);
c_dB = 20*log10(c_gain);
c_angle = 180/pi*angle(c_trans);
%%%%%%%%%%%%%%%%%%%%%%
% Error evaluation / 誤差評価
%%%%%%%%%%%%%%%%%%%%%%
% Squared error
err = c1*(m_dB-c_dB)*(m_dB-c_dB)' + c2*(m_angle-c_angle)*(m_angle-c_angle)';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% COEFF Mapping from variable x to coefficient of transfer func
% 変数 x と伝達関数パラメータの関連付けをする関数
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
if numel(x) == 3
b1 = 0; % Temporary ignored
else
b1 = x(4);
end
b0 = x(3);
a2 = 1; % Set it as 1
a1 = x(2);
a0 = x(1);
end
実行結果
Transfer Function Parameters:
b1: 0
b0: 2193.843
a2: 1
a1: 45.2179
a0: 2196.2001
System Parameters:
Sum of squared error: 0.18821
Natural angular frequency: 46.8636 rad/s
Damping factor: 0.48244
sys =
2194
--------------------
s^2 + 45.22 s + 2196
https://scrapbox.io/files/6902c79192909b068b4f07bc.png