概要
イナータンスを実験同定する必要があり，lsqcurvefitでカーブフィットさせるコードを書いたのでその備忘録。
やったこと
上図の様な1自由度マスバネダンパモデルでのイナータンスを導出する。
運動方程式を書くと
Mx''+Cx'+Kx=Fe^{j\omega t}
フーリエ変換して
{-\omega^2}MX(\omega)+{j\omega}CX(\omega)+KX(\omega)=F
式変形してイナータンスを求める。
\frac{{-\omega^2}X(\omega)}{F} = \frac{1}{M+\frac{C}{j\omega}+\frac{K}{-\omega^2}}\\
コード
入力1：$x$（$[M,C,K]$）
入力2：$xdata$（n行1列の周波数データ）
出力：y（n行1列のイナータンス）
の関数dof1_Faを作成
function y = dof1_Fa(x,xdata)
% 1自由度イナータンスモデル作成
% x(1) = M
% x(2) = C
% x(3) = K
w = 2*pi*xdata;
A = x(1) + x(2)./(1i*w) + x(3)./(-w.^2);
X = A.\1;
y = 20 * log10(abs(X));
end
初期値$[M,C,K]=[1,50,20000]$としてlsqcurvefitでカーブフィットを実施。
%% 実験データの読み込み
xdata = [17 24 31 36 41 42 43 45 48 53 56]'; % frequency
ydata = [-25 -18 -9 -4 5 9 10 8 3 -1 -3]'; % inertance
%% 1自由度イナータンスモデル呼び出し
fun = @dof1_Fa;
%% 開始点を指定してモデルの当てはめ
x = lsqcurvefit(fun,[1,50,20000],xdata,ydata);
%% 実験データと当てはめた曲線をプロット
freq = linspace(xdata(1),xdata(end));
plot(xdata,ydata,'ko',freq,fun(x,freq),'b-')
xlabel('Frequency [Hz]')
ylabel('Inertance [dB]')
legend('Data','Fitted Curve')
title('Data and Fitted Curve')
初期値はえいやで決めており，値によっては以下警告が出る。
lsqcurvefit は関数評価の制限、options.MaxFunctionEvaluations = 3.000000e+02
を超過したため、停止しました。
結果
x(1) = 3.34899555583242;
x(2) = 85.2344191122035;
x(3) = 247375.437493596;
で
実験データ○に対し，同定結果は青線となった。
