“使用K-Mean均值算法進(jìn)行聚類只是一種判斷干擾因素的手段，不是最終目的”

近期有幾個(gè)特別熱衷于調(diào)參的小伙伴問我，“那些大神/老師總說通過階躍響應(yīng)曲線來調(diào)PID參數(shù)并不正確，是這樣嗎”。

我難以揣測這背后到底有著多少對于數(shù)學(xué)的惡意，也無法猜測這其中摻雜了多少不懂裝懂的尷尬。但是既然大家對這些工具的惡意這么大，索性今天也就不計(jì)字?jǐn)?shù)，給大家粗略解釋一下PID-Analyzer和PIDtoolbox的階躍響應(yīng)形狀不同的原因，以及其背后的反卷積原理。

01

—

為什么同一個(gè)LOG會(huì)出現(xiàn)截然不同的結(jié)果？

簡單來說，這是因?yàn)?strong>PID-Analyzer直接輸出全部階躍響應(yīng)幀的平均結(jié)果，而PIDtoolbox啟用了根據(jù)穩(wěn)定時(shí)長篩選階躍響應(yīng)幀的算法來控制分析數(shù)據(jù)的質(zhì)量。

在PIDtoolbox內(nèi)看起來還算正常的log，在PID-analyzer里則完全超調(diào)查閱代碼可知，PID-analyzer在進(jìn)行維納反卷積之后直接對全部階躍響應(yīng)幀進(jìn)行了求平均，并進(jìn)行了圖像繪制。def?wiener_deconvolution(self,?input,?output,?cutfreq):??????#?input/output?are?two-dimensional?輸入輸出信號(hào)均為二維信號(hào)??pad?=?1024?-?(len(input[0])?%?1024)????????????????????????#?padding?to?power?of?2,?increases?transform?speed?將?padding設(shè)置為2的冪次數(shù)以提高變換速度??input?=?np.pad(input,?[[0,0],[0,pad]],?mode='constant')????#?指定輸入信號(hào)及模式，下為輸出信號(hào) output = np.pad(output, [[0, 0], [0, pad]], mode='constant') H = np.fft.fft(input, axis=-1)??G?=?np.fft.fft(output,axis=-1)??freq?=?np.abs(np.fft.fftfreq(len(input[0]),?self.dt)) sn?=?self.to_mask(np.clip(np.abs(freq),?cutfreq-1e-9,?cutfreq))??len_lpf=np.sum(np.ones_like(sn)-sn)??sn=self.to_mask(gaussian_filter1d(sn,len_lpf/6.))??sn=?10.*(-sn+1.+1e-9)???????#?+1e-9?to?prohibit?0/0?situations??Hcon?=?np.conj(H)??deconvolved_sm?=?np.real(np.fft.ifft(G?*?Hcon?/?(H?*?Hcon?+?1./sn),axis=-1))??#?此為反卷積核心公式??return?deconvolved_sm??#?將階躍響應(yīng)結(jié)果返回輸出def?stack_response(self,?stacks,?window):??#?繪制階躍響應(yīng)曲線 inp?=?stacks['input']?*?window??outp?=?stacks['gyro']?*?window??thr?=?stacks['throttle']?*?window??deconvolved_sm?=?self.wiener_deconvolution(inp,?outp,?self.cutfreq)[:,?:self.rlen]??delta_resp?=?deconvolved_sm.cumsum(axis=1)??max_thr?=?np.abs(np.abs(thr)).max(axis=1)??avr_in?=?np.abs(np.abs(inp)).mean(axis=1)??max_in?=?np.max(np.abs(inp),?axis=1)??avr_t?=?stacks['time'].mean(axis=1)??return?delta_resp,?avr_t,?avr_in,?max_in,?max_thr

摘自PID-Analyzer_v0.52.py

而PIDtoolbox在完成反卷積之后，使用stepinfo()函數(shù)對階躍響應(yīng)質(zhì)量進(jìn)行了提取，并隨后以穩(wěn)定時(shí)間為指標(biāo)對所有階躍響應(yīng)幀進(jìn)行了篩選，以完成質(zhì)量控制過程。如下所示。j?=?0；rateHigh =?0；for?i=1:size(SPseg,1) waitbar(i/size(SPseg,1),hw,['computing?step?response?functions...?']);? a=fft(GYseg(i,:).*hamming(length(GYseg(i,:)))');%?output,?use?hann?or?hamming?taper b=fft(SPseg(i,:).*hamming(length(SPseg(i,:)))');%?input,?use?hann?or?hamming?taper G=a/length(a); H=b/length(b); Hcon=conj(H);????? imp=real(ifft((G .* Hcon) ./ (H .* Hcon + .0001)))'; % impulse response function, .0001 to avoid divide by 0 impf = smooth(imp, lograte*10); % minor smoothing resptmp(i,:) = cumsum(impf); % integrate impulse resp functions resptmp(i,:)=resptmp(i,:)-resptmp(i,1);??a=stepinfo(resptmp(i,1:wnd));?%?gather?info?about?quality?of?step?resp?function?抓取階躍響應(yīng)質(zhì)量??if?a.SettlingMin>.5?&&?a.SettlingMin<=1?&&?a.SettlingMax>1?&&?a.SettlingMax<2?%Quality?control?質(zhì)量控制過程????j=j+1; stepresponse(j,:)=resptmp(i,1:1+wnd);? if?max(abs(SPseg(i,:)))?>=?rateHighThreshold rateHigh(j,:)=1; else rateHigh(j,:)=0; end end t = 0:1/lograte:StepRespDuration_ms;% time in ms end摘自PIDtoolbox/PTstepcalc.m也就是說，某些情況下，對階躍響應(yīng)幀進(jìn)行質(zhì)量篩選，是導(dǎo)致同一段LOG在兩個(gè)工具內(nèi)跑出完全不一樣結(jié)果的主要原因。

02

—

那PIDtoolbox的結(jié)果還是正確的嗎？

部分是 - 前提是要理解為何我們能計(jì)算出階躍響應(yīng)。已知輸入輸出信號(hào)的時(shí)域曲線，將其變換到復(fù)頻域后，使用(輸出/輸入)便可求得系統(tǒng)的傳遞函數(shù)，再將其變換回時(shí)域即可得沖激響應(yīng)。對其積分即可得到階躍響應(yīng)。

這在數(shù)學(xué)上稱為反卷積。反卷積是毫無對錯(cuò)可言的 - 它本來就是可推導(dǎo)的正確的數(shù)學(xué)定理。問題出現(xiàn)在質(zhì)量控制部分。有一部分人認(rèn)為系統(tǒng)的所有階躍響應(yīng)幀均可反映系統(tǒng)的實(shí)際性能，有一部分人則不這么認(rèn)為。

認(rèn)為所有階躍響應(yīng)幀均可反映系統(tǒng)的實(shí)際性能的人的考慮點(diǎn)是這樣的：系統(tǒng)的PIDF設(shè)計(jì)之初就是用來完成穩(wěn)定問題和抗擾動(dòng)問題，所以由于系統(tǒng)外非線性擾動(dòng)造成的異常響應(yīng)幀也應(yīng)視為有效幀。洗槳、強(qiáng)風(fēng)這些情況下的響應(yīng)幀也應(yīng)被計(jì)算在內(nèi)。

另外一部分人認(rèn)為，系統(tǒng)的設(shè)計(jì)前提是假定系統(tǒng)是局部線性，我們?nèi)粘ｏw行也基本是處于局部線性段，而反卷積也是針對線性系統(tǒng)，所以我們應(yīng)該只關(guān)注于局部線性部分的響應(yīng)性能。而前饋確實(shí)會(huì)給系統(tǒng)帶來一個(gè)相位超前的加速信號(hào)，這個(gè)信號(hào)會(huì)顯著改變系統(tǒng)的階躍響應(yīng)曲線。所以，要分析系統(tǒng)自身的階躍響應(yīng)曲線，確實(shí)不應(yīng)該使用前饋。

無論如何，反卷積算法在局部線性部分確實(shí)是科學(xué)有效的計(jì)算工具。但是，穿越機(jī)是一個(gè)由線性系統(tǒng)和非線性系統(tǒng)相互耦合的復(fù)雜系統(tǒng)，只使用反卷積階躍響應(yīng)曲線法來作為判斷PIDF參數(shù)的唯一判據(jù)顯然是不正確的。

03

—

那么這種算法本身的可靠性如何？

要評估算法自身的可靠性，就要分析全部階躍響應(yīng)幀，就需要找到一種科學(xué)有效的方式來對階躍響應(yīng)幀進(jìn)行歸一化統(tǒng)計(jì)。在數(shù)據(jù)分析領(lǐng)域，使用K-means平均算法來對信息進(jìn)行聚類分析是常見的數(shù)據(jù)分析方法。下面我們便使用K-means對原始階躍響應(yīng)幀進(jìn)行聚類分析。

圖片源自K-means_clustering維基百科第一步，打印出階躍響應(yīng)幀的數(shù)量。然后先對對數(shù)據(jù)進(jìn)行白化(Whiten)處理，之后再調(diào)用Kmeans算法對其聚類，并在在歐幾里得空間里標(biāo)示出數(shù)據(jù)分布與聚類結(jié)果的質(zhì)心位置。step_response[j,:]?=?respprint("We now have %d different step responses" %j)from scipy.cluster.vq import vq,kmeans,whiten# Set the number of clusters# In my test, it seems that the best value is 6GROUP_COUNT = 6# Pass step resonse to featuresfeatures?=?step_response[0:j,:]# Whiten datawhitened?=?whiten(features)# Find N clusters in the datacodebook, distortion = kmeans(whitened, GROUP_COUNT)print("The distortion of the signal is" %distortion)print(distortion)result?=?vq(whitened,?codebook)# Display centroids in Euclidean spaceplt.scatter(whitened[:,0],whitened[:,1])plt.scatter(codebook[:,0],codebook[:,1],c='r')plt.show()然后，標(biāo)示出各聚類結(jié)果對應(yīng)的樣本數(shù)量，這與PIDtoolbox的樣本數(shù)量N的意義基本一致。在對階躍響應(yīng)幀進(jìn)行處理時(shí)，同時(shí)統(tǒng)計(jì)輸入信號(hào)(搖桿輸入信號(hào))，并在最終繪制圖像時(shí)繪制與聚類后的階躍響應(yīng)曲線相對應(yīng)的輸入信號(hào)曲線，以判斷搖桿輸入信號(hào)對階躍響應(yīng)結(jié)果的影響。# Display clustered step response curvesresp_cluster?=?np.zeros((GROUP_COUNT,?segment_length))#print(resp_cluster)setpoint_cluster = np.zeros((GROUP_COUNT, segment_length))#print(setpoint_cluster)member_count = np.zeros(GROUP_COUNT)for i in range(j): codebook_idx = result[0][i] resp_cluster[codebook_idx] = resp_cluster[codebook_idx] + step_response[i] setpoint_cluster[codebook_idx] = setpoint_cluster[codebook_idx] + input_segments[i] member_count[result[0][i]] = member_count[result[0][i]] + 1for i in range(GROUP_COUNT): plt.subplot(GROUP_COUNT,2,i*2+1) plt.plot(resp_cluster[i]/member_count[i]) plt.subplot(GROUP_COUNT,2,i*2+2) plt.plot(setpoint_cluster[i]/member_count[i])print(member_count)通過對一段Betaflight 4.1的花飛LOG進(jìn)行分析，產(chǎn)生如下結(jié)果。其階躍響應(yīng)幀數(shù)量為1069個(gè)，聚類結(jié)果在歐幾里得空間里分布的如下圖所示，藍(lán)點(diǎn)為白化后的分布點(diǎn)，紅點(diǎn)為聚類結(jié)果質(zhì)心。

聚類結(jié)果如下所示，其中左側(cè)為階躍響應(yīng)曲線；右側(cè)為與之對應(yīng)的輸入信號(hào)曲線；橫軸2000代表樣本數(shù)量，實(shí)際表示1秒的時(shí)間；自上到下的樣本數(shù)量分別為469、15、398、6、109、72。

與之對應(yīng)，下面是此條LOG分別在PTB與PID-Analyzer內(nèi)分析出來的階躍響應(yīng)曲線。

可以看到，經(jīng)K-means聚類后，樣本數(shù)量為469與398的階躍響應(yīng)曲線與PTB大致相同。其樣本數(shù)量最多，大致反應(yīng)了飛行時(shí)的絕大多數(shù)情況下的階躍響應(yīng)曲線。另外也可以看到，“奇怪”的階躍響應(yīng)常發(fā)生于猛烈打桿后收桿歸中時(shí)。我們應(yīng)在飛行中盡量避免產(chǎn)生此類打桿操作。

04

—

結(jié)語

盡管K-means算法可以將階躍響應(yīng)進(jìn)行聚類，產(chǎn)生一系列較為科學(xué)的結(jié)果，但仍需注意，K-means并不是一個(gè)穩(wěn)定的算法，同一條LOG上可能會(huì)產(chǎn)生多種聚類結(jié)果。雖然K-means的結(jié)果更加復(fù)雜難以直觀判斷，但是它卻可以將各類階躍反應(yīng)直觀繪制出來。經(jīng)驗(yàn)豐富的技師可以根據(jù)聚類結(jié)果，對PIDF參數(shù)的適配程度有一個(gè)更完整的認(rèn)知；執(zhí)行能力強(qiáng)大的飛手可以根據(jù)聚類結(jié)果對應(yīng)的搖桿輸入情況，調(diào)整自己的打桿習(xí)慣，使飛行時(shí)飛機(jī)的動(dòng)態(tài)范圍總保持在局部線性區(qū)間內(nèi)。

圖源GmyFC很抱歉將大家?guī)胍粋€(gè)更加復(fù)雜的、尚未有確切定論的話題內(nèi)，但我們的本意是打消大家對工具的算法的惡意。希望諸位讀者在讀完本文之后可以完成我們的初始目標(biāo)(打消惡意)。畢竟，真正的數(shù)學(xué)意義上的線性非線性并不是什么順滑不順滑；反卷積算法也不是什么僅僅理論可行實(shí)踐掉鏈子的雞肋算法；真正的PID階躍響應(yīng)結(jié)果也并不會(huì)如“大神”所言“我打個(gè)桿就能改變這條曲線”。但是切莫不可將任何一種算法的結(jié)果直接當(dāng)成唯一的判據(jù)。正如本程序的幕后作者所言，“使用K-Mean均值算法進(jìn)行聚類只是一種判斷干擾因素的手段，不是最終目的”。

---

圖文：市民姜先生

排版：SinCerely

程序：ColinNiu、市民姜先生

審核：ColinNiu、山鷹_SamFisher、SinCerely、輝光管

如需轉(zhuǎn)載請后臺(tái)留言

大神請勿直接抄襲

分享給朋友或朋友圈請隨意哦

---

參考資料：

1.?反卷積和信號(hào)復(fù)原,?鄒謀炎,?國防工業(yè)出版社

2.?K-Mean-Cluster維基百科,?https://en.wikipedia.org/wiki/K-means_clustering

3.?PID-Analyzer源碼

4. PIDtoolbox源碼

5. pyPIDtoolbox源碼

總結(jié)

以上是生活随笔為你收集整理的阶跃函数卷积自己_如何看待大神/老师说“阶跃响应曲线调参并不准确”的观点？...的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。