charter(JJ1PID)

MATLABで非線形曲線近似する

非線形曲線近似する最低限のソースコードを載せました

Case.1 減衰振動のデータの近似

 y = 12.5e^{-0.2t}\sin(t)で減衰振動する物体の時刻tでの位置yの計測値、tdataとydataを考えます


MATLAB 非線形近似 減衰振動

%  tdata: 0から20まで100等分した値
tdata = linspace(0, 20, 100);

% ydata: 計測して得られた物体の位置(を想定)
% ここでは理想のモデル式に乱数のノイズを足しています
ydata = 12.5*exp(-0.2*tdata) .* sin(tdata) + 0.3*randn(size(tdata));


% 係数ベクトルtと、tdataを入力とする関数を作成
% 考えているモデル式の係数の部分を係数で置き換えます
fun = @(t, tdata) t(1)*exp(t(2)*tdata) .* sin(t(3)*tdata);

% 係数ベクトルxの各要素の初期値(なんとなく予想できる値)
% 近ければOKです
t0 = [12, -0.1, 1.1];

% 各係数が出力されます
t = lsqcurvefit(fun, t0, tdata, ydata);
% t=  [12.4747   -0.1979    0.9979]だったとします

% 計測値を黒の丸で、近似曲線を青の実線でプロットします
% plot(tdata, ydata, 'ko' ,tdata, fun(t, tdata), 'b-');
% なら計測値と近似曲線を一気にプロット出来ますが、近似曲線がtdataの値ごとの点を結んでいるためカクカクになります
% そこで、近似曲線の式を別に新しく定義し、重ねます
plot(tdata, ydata, 'ko');
hold on;
f = @(t)  12.4747* exp(-0.1979* t) .* sin(0.9979* t);
fplot(f, [0 20], 'b-');

% 以下、包絡線のプロットを付け加えています
f1 =@(t) 12.4747 * exp(-0.1979*t);
f2 =@(t) -12.4747 * exp(-0.1979*t);
fplot(f1, [0 20]);
fplot(f2, [0 20]);
hold off;


あとはプロパティインスペクターでグリッド入れたり、ラベルやタイトルを入れます


Case.2 LPFの利得の周波数特性

利得の周波数特性が G = \frac{1000}{j\omega+1000}と表されるLPFについて、その計測値の近似曲線を求めます
dB単位では G = 20\log_{10}{\frac{1000}{j\omega+1000}}\,\rm[dB]となります


LPF 利得 周波数特性

% 10から10^6までログスケールで等間隔な100個の点ω [rad/s]
omega = logspace(1,6,100)

% Gdata: 計測して得られた利得(を想定)
% ここでは理想のモデル式に乱数のノイズを足しています
Gdata = 20*log10(abs(1000./(i*omega+1000))) + randn(100,1)


G = @(k,omega) 20*log10(abs(k(1)./(i*omega + k(2))))


% 実際にフィッティングするときは大雑把でいいです
k0 = [990,1020]
k = lsqcurvefit(G,k0,omega,Gdata)


% 片対数グラフを作成
% 今回は近似曲線が緩やかで直線部分が長いのでカクカクしません
semilogx(omega,Gdata,'ko',omega,G(k,omega),'b-')

Case.3 インパルス電圧印加時の火花電圧の推定

棒ー平板電極にインパルス電圧を印加したときの火花電圧を推定します
直流や交流のときの火花電圧を知りたいのであれば、少しずつ印加電圧を上げていき火花放電が起きたときの値を使えばいいですが、インパルス電圧を印加する場合は印加電圧を変えてそのときの火花放電率を計測して推定する必要があります
その方法には補間法や昇降法がありますが、ここでは火花放電率Pは印加電圧Vに関する正規累積分布に従うとして、正規累積分布にフィッティングし、その結果から火花放電率と標準偏差を推定してみます


MATLAB 非線形近似 火花放電率 不平等電界

% V: 印加した電圧のベクトル
V = [74.121; 88.9452; 93.8866; 96.3573; 97.59265; 98.828; 101.2987; 103.7694; 106.2401; 108.7108; 111.1815; 113.6522; 123.535] 
% P: その電圧での火花放電率  
P = [0; 0; 0.2; 0.2; 0.2; 0.1; 0.4; 0.7; 0.7; 0.8; 1; 1; 1]

% 係数ベクトルkと、入力vの関数を作成
% normcdf(x,m,v)は平均値m, 標準偏差vの正規累積分布のxのときの値を出力します
% k(1)が平均値, k(2)が標準偏差を表します
p = @(k,vdata) normcdf(vdata,k(1),k(2)); 

% 係数ベクトルxの各要素の初期値(なんとなく予想できる値)
k0 = [100 10];

% 各係数が出力されます
k = lsqcurvefit(p, k0, V, P)
% k = [102.5146  5.7996]だったとします
% これは平均値が102.5146 , 標準偏差が5.7996の累積正規分布に近似できたことを意味します

% 計測値と近似曲線を重ねます
plot(V, P, 'o');
hold on;
% 近似曲線の式を作り直します
% ↑キーでモデル式定義したときの式に遡り、@(k,vdata)を@(vdata)に変えるだけです
p = @(vdata) normcdf(vdata,k(1),k(2)); 
fplot(f);


% テキストボックスに数式を入れるとき数式に翻訳してくれないときがあるので、先にコマンドで書き込みます
% 数式を書き込めたらあとは数式をダブルクリックして編集できます
text(90,0.9,'$$ P=\frac{1}{2}\rm{erfc}\left(-\frac{x-102.5}{\sqrt{2}\times 5.8}\right) $$','Interpreter','latex');

% あとはプロパティインスペクターでグリッド入れたり、ラベルやタイトルを入れます

hold off;


フィッティングの結果、計測した条件での火花電圧は102.5±5.8 kVと推定できました。


この方法は昇降法や補間法より正確ですが、100回以上印加する必要がありかなり大変ですし、印加してるうちに供試物がストレスで火花電圧が変動する可能性もあります
そのため、より少ない回数で火花電圧を推定する補間法や昇降法が一般に使われます


グラフ中のerfcとは相補誤差関数で、標準正規分布の累積分布関数を意味します


参考
正規累積分布関数 - MATLAB normcdf - MathWorks 日本
誤差関数とは - simonの開発日記

おわりに

上の例では標本値を作るところから始めていますが、実際にはデータをテーブル形式でインポートし(EPLとする)、次のようにすればいいです

xData = EPL.x
yData = EPL.y

(そのままk = lsqcurvefit(p, k0, EPL.V, EPL.P)としてももちろん良いです。その方が余計に変数が増えないので筆者はよくそうします。)


もっと詳しく知りたい方はこちらへ
非線形曲線近似 (データ適合) 問題を最小二乗近似的に解く - MATLAB lsqcurvefit - MathWorks 日本


線形・多項式モデルの近似はこちらへ
www.charter-blog.com