Help us understand the problem. What is going on with this article?

MATLABで簡単システム同定(簡易)

データ取得

まずは適切な入力を作成して入出力のデータを得る。

Chirp信号作成

簡単に作れて,かつそこそこの同定をするならChirp信号を用いる。
simulinkのブロックはなんか挙動が思ったように行かなかったので,Mfunctionで作成。

Linearな信号で良ければ1行で作成可能。
下記wikipediaの記事を参照。
https://ja.wikipedia.org/wiki/%E3%83%81%E3%83%A3%E3%83%BC%E3%83%97%E4%BF%A1%E5%8F%B7

  • u: clockから得た時間
  • f0: 最小の周波数
  • f1: 最大の周波数
  • Tchirp: 一度のsweep時間
  • phi0: 初期位相(0でよさそう)
sin(phi0+2*pi*(f0*rem(u,Tchirp)+(f1-f0)/2/Tchirp*rem(u,Tchirp)^2))

こんなかんじ。rem関数を使って周期的な動作を可能にしているのがミソ。

解析

そんでもって得たデータをin,outとする。

以下のがその一例。
(https://qiita-image-store.s3.amazonaws.com/0/188863/1243447c-53d8-b6a2-4eee-6a1d2d42dfa3.png)

周波数応答の同定

freqは当然周波数だが,winsizeは得られたデータをどのサイズで窓関数に打ち込むかを与える。
感覚的にはTchirp*freqで良さそうだが,間違ってたら指摘してください。

[PXX,FREQ] = tfestimate(in,out,winsize,[],[],freq);
Txx = frd(PXX,FREQ*2*pi,1/freq);

2行目でFREQに2*piをかけているのはHzからradに変換する必要があるため。どこかパラメータ設定すればこんなことしなくて済んだ気もする。

ボード線図を引いて見ればキレイに表示されていることと思う。

Bode(Txx)

また,Coherenceを取得するのもコマンド一つで済む。

[COHER,cFREQ] = mscohere(in,out,winsize,[],[],freq);
semilogx(cFREQ,COHER);

coherenceが小さいほどあんまり信用できないデータということになる。

お手軽フィッティング&システム同定

invfreqsという関数を用いることで簡単にフィッティングができる。学生はmatlabに課金すべき。

コマンドの使い方は以下の通りで先程tfestimateで生成したPXXとFREQを用いる。(必要に応じてradに変換すること!)

[num den] = invfreqs(PXX,FREQ*2*pi,nd,dd,weight,iter)
  • 3,4番目のnd,ddは分子と分母の字数で1次/2次なら,1,2と綴る。![octave-online-line-30.png]
  • weightは重み関数である。Coherenceから作成するとスムース。
  • iterはイテレーション回数。大体10回やっておけばまぁまぁ近い値になる。

重み関数作り方例

周波数帯がminf~maxfの間の情報のみかつ,Coherenceが0.95以上の周波数応答を用いてfittingするとする。
こんなかんじに書ける。簡単でしょう。

weight = zeros(size(COHER));
weight(COHER>0.95) = 1;
weight(cFREQ >= fmax)= 0;
weight(cFREQ <= fmin)= 0;

Sample Code(バグあり)

今matlabが消えたPCでやってるので後でちゃんと内容を確認する。あとで。

%% simulation set
close all
freq = 1000; %1000hz

%% 0. Define system
% system (ground truth) 
sys = tf([0.6],[0.02,10]);
figure;bode(sys);

%% 1. read data and set in/out data example

% make chirp signals
time = 0:1/freq:10;
CHP=chirp(time,0.5,2,50);
out=lsim(sys,CHP,time);
in=CHP;

%%load('result.data');
%time = result(:,1);
%in = result(:,2);
%out = result(:,3);

figure(1); plot(time,in,time,out);
xlabel('time[s]')
ylabel('in[A]/out[rad/s]')
legend('In','Out')
grid on;

%% 2. Get bode plot with tfestimate
windowsize = 2000;% this is not nessesary
[PXX,FREQ] = tfestimate(in,out,windowsize,[],[],freq);
Txx = frd(PXX,FREQ,1/freq);

figure(2); 
%opt = bodeoptions('cstprefs');
%opt.PhaseWrapping = 'on';
%opt.PhaseWrappingBranch = -180;
%bode(Txx,opt);
bode(Txx);
title('Freq responce')

%% 3. Get coherence mscoherence
[COHER,cFREQ] = mscohere(in,out,windowsize,[],[],freq);

figure(3); semilogx(cFREQ,COHER);
title('coherent')
xlabel('Freq [Hz]')
%% 4. fitting1
% Chirp Sweep Range
fmax = 50;
fmin = 0.5;

% choose weight from coherence
weight = zeros(size(COHER));
weight(COHER>0.95) = 1;
% use only 0.5--50 Hz
weight(cFREQ >= fmax)= 0;
weight(cFREQ <= fmin)= 0;

% estimate transfer function !
[num den]=invfreqs(PXX,2*pi*FREQ,0,1,weight,10);

% figure;bode(tf(num,den));
fitsys = tf(num,den);

figure(4);bode(sys,fitsys);
legend('ground truth','Estimated')
xlim([fmin fmax])

Why do not you register as a user and use Qiita more conveniently?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
Comments
Sign up for free and join this conversation.
If you already have a Qiita account
Why do not you register as a user and use Qiita more conveniently?
You need to log in to use this function. Qiita can be used more conveniently after logging in.
You seem to be reading articles frequently this month. Qiita can be used more conveniently after logging in.
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away