利用系統波德圖與奈氏圖設計分數階PID控制器

國

立

高

雄

應

用

科

技

大

學

電

機

工

程

系

碩

士

班

碩

士

論

文

利用系統波德圖與奈氏圖設計分數階 PID 控制器

Fractional Order PID Controller Design via

Bode Diagram and Nyquist Plot

研 究 生 ： 周永吉

指導教授 ： 李

慶

祥 博士

利用系統波德圖與奈氏圖設計分數階

PID 控制器

Fractional Order PID Controller Design via Bode Diagram and

Nyquist Plot

研 究 生：周 永 吉 Student：Yung-Chi Chou 指導教授：李 慶 祥 Advisor：Ching-Hsiang Lee

國立高雄應用科技大學

電機工程系碩士班

碩士論文

A Thesis

Submitted to

Department of Electrical Engineering

National Kaohsiung University of Applied Sciences

in Partial Fulfillment of the Requirements

for the Degree of Master of Science

in

摘要

本論文藉由在系統奈氏圖中移動頻率響應曲線，使其產生在



G j

( )



 

180

 和G j(



) 1 處，兩個點對點之間的關係式，再藉此取得聯立方程組以求得分數 階控制器參數。本人利用控制器參數調整方法設計分數階PID 控制器，使系統達 到期望的增益邊界與相位邊界，其中系統在期望的交越頻率(crossover frequency) 處，相位曲線是平坦的，而且對增益變化具有強健性。參數設計由粒子群演算法 找出最佳參數，確保得出的控制器參數不會只是區域最佳解。分數階PID 控制器 與傳統整數階PID 控制器相比，分數階控制器的優點為具有五個可調參數，即

,

,

,

p i d

K

K K



_{(積分階數)和}



(微分階數)，從而增加設計的自由度並且可以使系 統獲得較強的強健性。根據數值模擬，本論文的設計方法，可以求出良好的可調 參數。 關鍵字 ─ 分數階微積分；分數階PID 控制器；增益邊界；相位邊界；粒子群演 算法

Abstract

The Nyquist plot of the system is shifted such that a point-to-point relations at

( )

180

G j







 

and G j(



) 1 conditions as well as a system of algebraic equations can be obtained to solve for the parameters of an FOPID controller. The parameters

,

,

,

p i d

K

K K



_{(the integrating order) and}

_

_{(the derivative order) of an FOPID}

controller are properly selected to achieve desired gain margin and phase margin, to obtain a flat phase curve at the crossover frequency, and to increase the robustness property against the plant parameter variations. The searching of the optimal controller parameters is accomplished using PSO (particle swarm optimization) technology. Numerical simulations demonstrated the optimal parameters can be obtained using this method.

Key Words ─ Fractional Calculus, Fractional PID Controller, Gain Margin, Phase Margin, Particle Swarm Optimization (PSO).

誌謝

光陰似箭，在高雄應用科技大學的研究所生涯轉眼間就結束，回想起這兩年間 求學過程的種種點滴，在此衷心感謝在此曾經鼓勵與支持我的師長、同儕與親友 們。 在此論文完成之際，首先要感謝我的指導教授李慶祥博士，給予我許多寶貴的 建議與指導，使能順利完成論文，讓我這段求學過程中受益良多。同時，感謝在 論文口試期間，高雄海洋科技大學凃文福博士與本校李宗恩博士和李慶祥博士在 論文上的指正與建議，以及費心審閱本論文，方能使之更臻完整。 接著，感謝非線性控制實驗室的浩倫學長、至宗學長、京寰學長、亭維學姊、 偉煌、明偉、韶鴻、得益和繼哲，以及智慧自動化實驗室、控制理論與應用實驗 室的明洧學長、翊展學長、忠義學長、信杰學長、若軒學長、建安、順利、智翔、 河仲、品高、彥柔、冠龢、建霖和培峰，他們的陪伴在我碩士生涯中增添了許多 色彩。 最後，感謝一直以來都支持我的家人，有他們的付出與照顧，讓我的求學路途 得以順利，毫無後顧得完成學業。僅以本文獻給我敬愛的家人及所有支持我的人， 願將這份成果和喜悅與各位一同分享。 周永吉 謹誌於 國立高雄應用科技大學 非 線 性 控 制 實 驗 室 中華民國 一百零四 年 七 月

目錄

摘要...i Abstract...ii 誌謝...iii 目錄...iv 圖目錄...vi 表目錄...vii 符號說明...viii 第一章、緒論...1 1.1 文獻回顧與研究動機...1 1.2 論文綱要...5 第二章、預備知識...6 2.1 系統波德圖簡介...6 2.2 系統奈氏圖簡介...6 2.3 分數階微積分定義...7 2.3.1 Grünwald–Letnikov 分數階微積分定義...8 2.3.2 Riemann–Liouville 分數階微積分定義...9 2.3.3 Caputo 分數階微積分定義...10 2.4 分數階微積分定義之間的關係...11 2.6 Ziegler–Nichols 參數調整方法...13

3.3 粒子群演算法...18 3.3.1 線性慣性權重...21 3.3.3 線性更新公式...23 3.4 分數階控制器：調整...24 3.4.1 分數階 PID 控制器：Yeroglu 調整法...25 3.4.2 Yeroglu 調整法的缺陷...28 3.5 改善對分數階 PID 控制器調整方法...29 第四章、例題模擬與討論...35 4.1 UAV 飛行姿態控制模型...35 4.2 使用 PSO 搜尋初始控制器...37 4.2.1 慣性權重的改善...37 4.2.2 粒子的優劣...38 4.3 例題模擬與系統響應...41 第五章、結論與未來研究方向...46 5.1 結論...46 5.2 未來展望...47 參考文獻...48

圖目錄

圖 3-1 控制器調整法示意圖 14 圖 3-2 20 圖 3-3 Yeroglu 調整法示意圖 25 圖 3-4 Deepak 調整法示意圖 29 圖 3-5 改良調整法示意圖 30 圖 4-1 系統波德圖(不好的粒子) 38 圖 4-2 40 圖 4-3 單位步階響應之比較 42 圖 4-4 系統開路波德圖(傳統控制器) 44 圖 4-5 系統開路波德圖(FOPI 和 FOPID 比較) 45

表目錄

表 2-1 PID 控制器參數調整方法 13

符號說明

p

K

控制器比例增益 i

K

_{控制器積分增益} d

K

_{控制器微分增益}  控制器積分階次



_{控制器微分階次} 0 p

K

初始控制器之比例增益 0 i

K

_{初始控制器之積分增益} 0 d

K

_{初始控制器之微分增益} p

K

期望控制器之比例增益 i

K

_{期望控制器之積分增益} d

K

_{期望控制器之微分增益}  _{期望控制器之積分階次}  _{期望控制器之微分階次}

u

T

振盪週期(臨界穩定) m



_相位邊界 m

g

_增益邊界 cp



相位交越頻率 cg



增益交越頻率

第一章、 緒論

1.1 文獻回顧與研究動機

使用非整數階的微分與積分的這個概念並非是新的發想。最早提到分數階微積 分的想法可以追溯到1695 年 Leibniz 給 L’Hospital 的信件[1]，至今這個話題已經 有300 多年的歷史了。分數階微積分理論的研究與發展主要在 19 世紀，由 1832 年Liouville [2]、1953 年 Riemann [3]和 1864 年 Holmgren [4]，還有 1730 年 Euler [5]、1772 年 Lagrange [6]等等，甚至還有人更早作出貢獻。

Liouville 將指數函數擴展到級數中，並如同 q 為正整數一般逐項定義此級數第

q 階的導數。Riemann 提出了針對定積分不同的定義，且可以應用於非整數指數

的冪級數。Grünwald [7]提到 Liouville 的方法受到限制，因為他的微分定義的出 發點為差商(difference quotient)的極限和 q 階導數的定積分公式。Krug [8]通過對 一般導數做柯西(Cauchy)積分公式動作，得出從 Liouville 的定義當中，Riemann 定積分必須被解釋為具有有限的下限，其中沒有可以區別_{下限值的存在。} 在 20 世紀已經有許多著名的貢獻做出分數微積分運算子(differintegral)的理論和 應用。Weyl [9]、Hardy [10]、Hardy & Littlewood [11-13]、Kober [14]以及 Kuttner [15]等人的研究是一些比較特殊但不失一般性的分數階微積分運算子，其性能屬 於Lebesgue 和 Lipschitz 這類的性能。Erdélyi [16-18]和 Osler [19]對任意函數給予 分數階微積分運算子的定義，並且Post [20]用差商定義廣義微分運算子 f D

 

，其

[35] ； 化 學 物 理 學 方 面 如 Somorjai & Bishop [36] ； 一 般 的 運 輸 問 題 (transport problems)如 Oldham[37]以及 Oldham & Spanier [38]。

分數階微積分的發展為各學科提供新的理論基礎，如在冶金、化工、電力和機 械等工業過程中都有應用。分數階微積分對於複雜的系統提供了更完善的數學模 型，在物理、生物工程、控制理論等方面也有許多應用。隨著工業的發展，對於 一些實際模型的建立有了更高的要求，分數階微積分的引入使數學模型變得更加 準確。 分數階控制系統既可以應用到整數階系統，也更適合分數階受控系統模型中。 若在複雜動態系統，使用分數階微積分方程建模要比整數階系統模型更加準確， 尤其是在物理、生物醫學等方面，分數階模型可以準確描述動態系統的特性。一 般對分數階模型而言須使用分數階控制器才能起到很好的控制效果。分數階控制 器增加了可調參數，其控制效果遠遠好過整數階控制器。分數階控制器不僅適用 於分數階系統模型，而對於整數階系統模型也能充分展現其優越性。

由於分數階微積分理論被應用到研究控制領域的時間不是很長，其研究方向包 括以下三個方面：其一，分數微積分逼近法的研究；其二，用分數階微分方程描 述較難通過整數階微分方程的控制系統；最後一個是研究分數階微分控制器和積 分控制器的實現。

近年來，分數階微積分已越來越多應用於機械系統、電力和生物工程。分數階 PID (fractional order PID，FOPID) 控制理論是基於傳統的 PID 控制理論的發展， 分數階控制器的設計與研究已經獲得越來越多的關注。分數階PID 控制器的想法 是由Podlubny [39]所提出。與傳統整數階控制器相比，分數階控制器的優點為具 有五個可調參數(

K K K

p

,

i

,

d

, ,

 

)，從而增加設計的自由度，並且可以實現更好 的控制品質。它能有效地抑制雜訊，對於系統模型的不確定性具有較強的強健性。 群智能優化，是指無智能的或具有簡單智能的個體透過群體協作和組織來表現 出群體智能優化行為的特性。大部分圍繞二大分支進行研發：

(1) 蟻群演算法(Ant Algorithms)：於 1992 年 M. Dorigo [40]觀察螞蟻群的搜尋食物 活動啟發而建立的演算機制模式。

(2) 粒 子 群 演 算 法 (Particle Swarm Optimization) ： 於 1995 年 J. Kennedy & R. Eberhart [41]由鳥群捕食行為的研究得到啟發而產生的最佳化演算法。

J. Kennedy & R. Eberhart 由鳥群捕食行為提出一個簡單的模擬： (1) 每個智能體(agent)會朝著最優個體接進。 (2) 每個智能體會記憶當前的最佳個體。 (3) 每個智能體會與鄰近個體共享與最佳個體相關的信息。 粒子群演算法的開發基於一些很簡單的概念，使其僅需要一般數學運算或使用 幾行程式編碼撰寫就可以實現。再加上其演算法整體有許多可以調整、設計的參 數，有許多學者提出方法來改善粒子群演算法，Shi & Eberhart 分別在 1998 年[42] 提出常數慣性權重(inertia weight)的想法，以及在 1999 年[43]提出線性遞減的慣性 權重改善常數慣性權重；Clerc [44]提出收縮係數的概念；Ratnaweera 等人[45]提 出分別對個體學習參數和群體學習參數線性化；Chatterjee & Sisrry [46]提出權重 係數以非線性方式遞減；Ko 等人[47]從非線性權重的概念，衍生出分別對個體學 習參數和群體學習參數非線性化。

1.2 論文綱要

本論文運用分數階控制器有多個可調參數的特性，使系統達到期望性能以及附 有強健性，在設計方法之中結合了粒子群演算法搜尋最佳參數組合，文章結構如 下：第一章緒論，先提到研究動機與歷史回顧；第二章簡單介紹分數階微積分定 義與Ziegler–Nichols 參數調整方法；第三章方法章節中，分別介紹本論文所提出 的設計步驟與推導過程，以及粒子群演算法的使用簡介；第四章前半段為推導例 題的系統模型，將其代入本篇提出的調整方法中，並藉由例題模擬利用粒子群演 算法用來尋找參數所遇到的問題，然後針對問題來改善粒子群演算法的程序。後 半段為比較傳統控制器與分數階控制器對輸出的響應和頻率響應；第五章是本論 文的結論和未來展望。

第二章、 預備知識

2.1 系統波德圖簡介

波德圖是由兩張圖組合而成，其一幅頻圖為分貝值對頻率的變化，另一張相位 圖為相位對頻率的變化。利用波德圖可以看出在不同頻率下，系統增益的大小及 相位，也可以看出增益大小及相位隨頻率變化的趨勢。 在控制器設計條件中，給予對應於波德圖曲線的各種性能指標(如相位邊界、增 益邊界等)可從中簡單得出設計條件。由於設計方法靈活性高且簡單，故在本論文 中選用章節3.1 的設計方法。

2.2 系統奈氏圖簡介

奈氏圖是將系統的頻率響應增益與相位以及座標的方式在複平面中會出曲線， 常在控制系統或信號處理使用，可以用來判斷系統是否穩定。閉路系統可以由其 開路系統的奈氏圖，配合奈氏穩定準則(Nyquist stability criterion)判斷其穩定性。 從頻率響應得出的設計方法可以用在有延遲的系統，或是傳輸函數不是有理函 數的系統，這階系統用其他方法很難分析。故本論文的控制器設計法是可以有效 地應用在這類系統上。

2.3 分數階微積分定義

分數階微積分作為建模工具的使用上擁有很好的效果，並且可以從數學邏輯上 證明系統的有效性。 分數階理論的發展主要在十九世紀，分數階微積分理論上為任意階微積分，也 就是階數可以為任意實數甚至可以為任意複數。分數階微積分的基本運算子為 a

D

t  ，其中_a和

_t

為運算子的上下限，_ 為微積分階次並可為任意複數，本討論為 實數的情況。

 

 

 

 

_{ }

, Re 0 1 , Re 0 , Re 0 a t t a d dt D d              _     



(2-1) 其中Re

 

 為_ 的實部，透過引入分數階運算子a

D

t  ，積分和微分可以被統合在 一起。 目前分數階微積分還沒有被廣泛地應用在生活中，原因之一為到現在為止還沒 有一個統一的定義。在分數階微積分理論發展過程中，有很多種分數階微積分定 義，以下舉出三種分數階微積分定義做簡單的介紹，分別為Grünwald–Letnikov 分數階微積分定義[48-49]、Riemann–Liouville 分數階微積分定義[50-54]和 Caputo 分數階微分定義[55]。

2.3.1 Grünwald–Letnikov 分數階微積分定義

對於任意實數 m，設 m 的整數部分為[m]，對函數 f t

 

的_ 階引入分數階微積 分運算子a

D

t  的定義為

 

 

 





 





 





  ( ) 0 0 ( ) 0 0 lim 1 1 lim 1 t h j a t _h j t h h j D f t h f t jh j j f t jh h j            _   _     _{ }          





(2-2) 其中



 

 







1 2 1 ! ! ! ! j j j j j                  _    (2-3) Grünwald–Letnikov 分數階微積分定義是從尋找 n 階導數和 n 階積分的統一性整 理出來的。首先對於充分連續函數 f t

 

整理出n 階導數



n N



的公式

 

( ) 

 





( ) 0 ₀ lim 1 t n h j n n h _j n f t h f t jh j         _{ }   



(2-4) 其中



 

 







1 2 1 ! ! ! ! n n n n n j n j j j n j           _    (2-5)

接著把(2-5)式中的 n 擴展到任意整數 p。當 p 為負整數n_時，

 

1 n n n j j                (2-6) 再把p 擴展到任意實數 _{便可得到(2-2)式。}

2.3.2 Riemann–Liouville 分數階微積分定義

Riemann–Liouville (RL)對分數階微分的定義如下， 對於任意實數m   ， m N1  m 

 

_

_

 





1 1 m _t a t m _a m f d D f t d m dt _t           



_ (2-7) 分數階積分的定義為

 

_{ }

 





1 1 t a t _a f I f t d t      _    



 (2-8) 其中

t



0 ,





R



RL 分數階微分和積分定義也可以整合到一個表示式中， RL 分數階微積分定義為

 

_

_

 





1 ( ) 1 m t a t _a m f d D f t d m dt _t      _      _{ }    



 (2-9) 其中m   ， m N1  m  。

2.3.3 Caputo 分數階微積分定義

Caputo 分數階微分定義為

 

_

_

 





( 1) 1 1 m t a t _a f D f t d t  







_

   



 (2-10) 其中   ， m Nm   ，0  1。 Caputo 分數階積分定義為

 

_{ }

 





1 1 t a t _a f D f t d t          



_ (2-11)

上述兩式也可以整合成起來，Caputo 分數階微積分定義為

 

_

_

 





( ) 1 1 t m a t _a m f D f t d m _t  







_

    



 (2-12) 其中m   ， m N1  m  。

2.4 分數階微積分定義之間的關係

Riemann–Liouville 分 數 階 微 積 分 定 義 和 Caputo 分 數 階 微 積 分 定 義 都 是 對 Grünwald–Letnikov 分 數 階 微 積 分 定 義 進 行 改 進 。 這 兩 者 之 間 的 差 異 分 別 為 ， Riemann–Liouville 分數階微積分定義是對 f t

 

的正非整數_ 階導數先進行_m 階積分，再進行m 階微分；Caputo 分數階微積分定義則是對正非整數 _階導數先 進行m 階微分，再進行m_{階積分。在階數} _{為副實數和正整數時，兩者之間} 的定義是等價的。 這兩者的定義也有不同用途，引入 Riemann–Liouville 分數階微積分定義可以簡 化分數階導數的計算；引入Caputo 分數階微積分定義可以讓拉氏轉換(Laplace transform)更簡潔，有利於分數階微分方程的討論。

2.5 分數階微積分的性質

就像整數階微積分一樣，分數階微積分也具有基本性質： (1) 解析函數 f t

 

的分數階導0D f tt

 

 對t 和_ 都是解析的。 (2) 當 n_{，n 為整數時，分數階微積分和整數階微積分的值完全一樣，而且}

 

 

0 0D f tt  f t 。 (3) 分數階微積分運算子為線性的，

 

 

 

 

0Dt a f t bg t a D f t0 t b D g t0 t  _{  }  _    (2-13) 其中a、b 為任意常數。 (4) 分數階微積分運算子滿足交換律，並且滿足疊加關係

 

 

( )

 

0Dt 0D f tt 0 Dt 0D f tt 0 Dt f t _  _  _  _       (2-14)

2.6 Ziegler–Nichols 參數調整方法

早期由 J.G. Ziegler 和 N.B. Nichols 提出[56]，故以此為命名。調整方法首先將 積分增益

K

i和微分增益

K

p設為0，接著把比例增益K從0 開始增加，直到輸出成 穩定週期震盪時，增益將達到臨界增益

K

_u，震盪週期為

T

_u。 除了原本 Ziegler–Nichols 參數調整方法外，McCormack [57]提到許多調整方法 是根據ZN 參數調整方法進行改進，得出不同用途的參數調整方法，表 2-1 為文 獻中所整理的調整法。在本篇中是使用原方法ZN 參數調整方法做應用。 PID 控制器的形式為

 

_c 1 1 _d i C s K T s T s    _   _   (2-13) 表 2-1 PID 控制器參數調整方法 調整方法 控制器參數 ZN

K

c



0.6

K T

u

,

i



0.5 ,

T T

u d



0.125

T

u PIAE

K

c



0.7

K T

u

,

i



0.4 ,

T T

u d



0.15

T

u SO–OV

K

c



0.33

K T

u

,

i



0.5 ,

T T

u d



0.33

T

u

0.2

,

0.5 ,

0.33

K



K T



T T



T

第三章、 方法介紹

本方法主要針對不穩定系統做控制，在設計步驟上分為兩個階段，如圖 3-1 所 示，第一階段為設計初始控制器使系統產生實線曲線；第二階段為調整控制器， 把實線調整到虛線。在第一階段中，由於本篇所使用的方法在尋找控制器參數的 自由度不足以直接到達虛線部分，因此需要分為兩階段來做設計。 圖 3-1 控制器調整法示意圖

第二階段為調整，此階段的方法對尋找分數階控制器的參數擁有較高的自由度 ， 但是因為此調整方法運用在不穩定系統上時，必須要有可以使系統穩定的初始控 制器，也就是說要先使系統能畫出實線部分，然後再調整初始控制器使其達到期 望的性能，將實線移動至虛線。

3.1 初始控制器設計規範

從 Wang 等人[60]所提出的分數階 PI 控制器設計方法作為基礎，將其方法推廣 成分數階PID 控制器，使其能套用到本論文初始控制器的設計。在此篇提到，對 於系統的穩定性和強健性，根據開路傳輸函數的相位和增益提出三個設計規範：

假設增益交越頻率(gain crossover frequency)為



c，和相位邊界(phase margin)為 m



(i) 相位邊界規範

[ (

_c

)]

[ (

_c

) (

_c

)]

_m

Arg G j





Arg C j



P j



  

 

(ii) 增益交越頻率的規範 ( _c)_db ( _c) ( _c)_db 0 G j  C j P j 

分數階 PID 控制器的形式如下所示， 0 0( ) 0(1 0 ) i p d K C s K K s s      _(3-1) 其頻率響應如下 0 0 0 0 0 0 0 0 0 ( ) (1 ( ) ( ) ) 1 cos sin 2 2 cos sin 2 2 p i d p i i d d C j K K j K j K K jK K jK                              _  _{ } _{ }            _{ } _{ }     (3-2) 相位角度和增益的推導為，將(3-2)式分為實部與虛部，再分別將它們做相位角 度和增益振幅的運算 0 0 0 0 0 sin( / 2) sin( / 2) [ ( )] arctan 1 cos( / 2) cos( / 2) i d i d K K Arg C j K K                     (3-3) 2 0 0 0 0 1 2 2 0 0 ( ) (1 cos( / 2) cos( / 2)) ( sin( / 2) sin( / 2)) p i d i d C j K K K K K                  _      _ (3-4) 將系統接上分數階 PID 控制器後，其開迴路頻率響應的相位與振幅(amplitude) 分別如下所示， 0 0 0 0 0 [ ( )] [ ( )] [ ( )] sin( / 2) sin( / 2) arctan [ ( )] 1 cos( / 2) cos( / 2) m i d i d

Arg G j Arg C j Arg P j

K K Arg P j K K              _                  (3-5)

0 2 0 0 0 1 2 2 0 0 ( ) ( ) ( ) 1 (1 cos( / 2) cos( / 2)) ( sin( / 2) sin( / 2)) ( ) p i d i d G j C j P j K K K K K P j                       _      _  (3-6) 透過章節 3.1 提到控制器的設計規範 (i)、(ii) 和 (iii)，將上述開迴路頻率響應 的相位與振幅(3-5)和(3-6)式套入後，再給予期望的系統性能，相位邊界



m以及增 益交越頻率



_c，便可得出符合系統性能的控制器參數

K

p0、

K

i0、

K

d0、和



。 其控制器參數是以何種方法得到，在 Wang 等人[60]的文獻中並沒有詳細規定 也沒有提到。因此在本碩論將三個設計規範套入粒子群演算法(particle swarm optimization，PSO)之中，透過粒子群演算法全域搜尋的性質來找出符合設計規範 的控制器參數。

3.3 粒子群演算法

粒子群演算法(particle swarm optimization，PSO)的理論基礎是將粒子(particle)模 擬程鳥類族群中的個體，並賦予每個粒子(個體)擁有記憶性，使其能夠透過粒子 群體中每個粒子之間的互動進而搜尋到最佳解。每個粒子皆擁有各自的位置與速 度，並在移動的過程中會產生經驗與記憶，每當個體移動的同時，其粒子便能依 照自身的經驗與記憶來學習調整移動方向。由於在粒子群演算法中是以多個粒子 同時移動的，同時以自身的經驗和其他粒子提供的經驗進行比對，藉此尋得最佳 解。 在本篇的應用中，每個粒子皆有與之對應的目標函數值 (objective function value)，每個粒子與目標函數值的誤差將作為判斷該粒子優劣的依據，依此決定 自身最佳解(當代最佳解)及群體最佳解(歷代最佳解)的更新動作，利用這樣的機制 不斷更新每個粒子的移動速度，反覆地在解空間中迭代搜尋。每次迭代的過程不 是完全隨機的，如果找到最佳解，將會以此為依據來尋找下一個解。

在 D 維空間中，設有 m 個粒子，第 i 個粒子的位置為xi 



x xi1, i2, , xiD



，第i 個粒子的速度為vi 



v vi1, i2, , viD



。將

x

_i代入目標函數中計算其適應值，衡量

x

_i 的優劣。假設搜尋迄今為止，單一粒子的最佳位置為Pi 



P Pi1, i2, ,PiD



  ，群體粒 子群的最佳位置為Pg  P Pg1, g2, ,PgD   ，參照Chou [61]講義中說明基本粒子群演 算法的粒子速度與位置的更新公式如下









1 1 1 2 2 t t t t t t id id id id gd id V  v c r P x c r P x _(3-7)

x

idt 1

x

idt

v

idt 1 

_

_

 (3-8) 其中i1, 2,_,m和d 1, 2,_,D，

r

1和

r

2為介於0 到 1 之間的隨機亂數，

c

1為個 體學習參數，

c

₂為群體學習參數，



為慣性權重。

圖 3-2 為基本粒子群演算法的運作過程，在本論文中，粒子的速度初值定為 0， 粒子的位置初值給予隨機變數。粒子適應值是把章節3.1 的設計規範(i)、(ii)和(iii) 做為目標函數，將對目標函數的誤差總和作為適應值的計算。在更新完粒子後， 判斷是否滿足終止條件，粒子的誤差總和若低於目標則跳出迴圈。

3.3.1 線性慣性權重

慣性權重關係著粒子搜尋解的範圍，如何設定其值與最佳解的搜尋有相當的關 連性。本篇中所使用的線性慣性權重為Shi & Eberhart [43]所提出的設計方法，他 們將常數的慣性權重以線性遞減的方式取代，

max

min max min

max

( )k iter iter ( )

iter

       _(3-9)

其中

iter

max為最大迭代次數；iter為當前代數也就是第k 代。

這個設計方法可以使粒子隨著迭代的次數增加來調整慣性權重值的大小，在搜 尋前期給予粒子群較大的全域搜尋能力，隨著代數的增加使慣性權重遞減，讓粒 子群在搜尋後期收斂為局部搜尋。在本篇的應用中，此方法能有效地使粒子收斂 至最佳解。

3.3.2 線性學習參數

上述的慣性權重對粒子群的影響為解的搜尋範圍，而學習參數對粒子群所影響 的方面為個體與群體之間的影響力。Ratnaweera 等人[45]提出了時變的個體學習 參數

c

1和群體學習參數

c

2，加入線性學習參數後的PSO 搜尋方式，在演算法最 初進行全域性的搜尋，最後收斂至群體最佳解之中。 在粒子打散後開始搜尋時，給予較大的個體學習參數值和較小的群體學習參數 值，隨著迭代次數的增加，在搜尋的後期給予粒子群較小的個體學習參數值和群 體學習參數值。也就是說，在搜尋的過程中，前期是以個體最佳解的影響較大， 到了後期開始往群體收斂則以群體最佳解的影響較大。 此方法的學習參數設計分別為，個體學習參數

c

1_{從最大值隨著迭代次數以線性} 遞減的方式減至最小值；相對的，群體學習參數

c

2_{則是從最小值以線性遞增的方} 式增加到最大值。 max

1 1min 1max 1min max ( ) iter iter ( ) c k c c c iter      _(3-10) max

2 2max 2min 2max max ( ) iter iter ( ) c k c c c iter      _(3-11)

3.3.3 線性更新公式

為了使 PSO 更有效率地找到最佳參數，本篇將線性公式融入更新公式之中，藉 由線性公式的特性使PSO 在求解初始控制器參數的過程能有效地往最佳解收斂。 將上述線性慣性權重和學習因子(3-9)、(3-10)和(3-11)結合 PSO 基本更新公式(3-7)和(3-8)，便可得出本篇所使用的線性更新公式









1 1 2 2 ( 1) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) i i i i g i V k k V k c k r P k X k c k r P k X k             (3-12)

(

1)

( )

(

1)

i i i

X k

 

X k V k





(3-13)

3.4 分數階控制器：調整

在 Deepak & Kumari [62]文獻中所使用的調整方法改善了 Yeroglu & Tan [63]的 控制器調整方法。原方法(以下稱 Yeroglu 調整法)存在著一些缺陷，其中影響最大 的問題在於此調整法中只調積分階數  和微分階數



來獲得期望的頻率響應，以 及無法運用在不穩定的系統。Deepak & Kumari [62]在改善後的方法中提到，此方 法要運用在不穩定系統中必須要有傳統PID 控制器使系統穩定，以此條件為前提 套入調整方法從中獲得期望的性能。

本篇中，吾人針對不穩定系統的部分進行修改，因為 Deepak & Kumari [62]提 到，要進行調整不穩定系統時，必須先找到一組PID 控制器使系統穩定後再調整 而在這當中沒有明確地提出要用什麼方法找出這組控制器，因此本篇將原本所使 用的傳統PID 控制器變換為分數階 PID 控制器，也就是將第一階段的設計方法結 合此處的調整方法進而得到期望性能的控制器。

3.4.1 分數階 PID 控制器：Yeroglu 調整法

文獻中提出多種分數階 PID 的調整方法[64-65]，在這部分介紹 Yeroglu 調整法。 圖3-3 中的曲線分別為接上控制器前後的頻率響應曲線，藉由調整控制器參數使 系統的響應曲線從實線移動到虛線的位置。控制器參數的設計是由點1 與點 2 的 關係式來進行調整。此方法是將Ziegler–Nichols 和 Åström–Hägglund 調整方法套 入分數階PID 控制器調整方法之中，也就是用來得出

K

p、

K

i和

K

d的初值。其中， p

K

_和

K

_i的初值是使用Ziegler–Nichols 調整方法給定，而

K

d的初值是使用Åström– Hägglund 調整方法給定，然後設定期望達到的相位邊界後，代入經由下述的推導 所得出的兩個非線性方程式進行調整。

1 ( ) j cp c G j e K    _(3-14) 其中， 1 ( ) c m cp K g G j   _(3-15) 和



cp_{為相位交越頻率。假設接上分數階PID 控制器的系統頻率響應曲線會通過} 點2，也就表示點 2 為對應接上控制器的系統頻率響應的相位邊界。再假設接上 控制器的系統增益交越頻率會等於未接上控制器的系統相位交越頻率



cg cp



_。 根據這些假設可將點2 表示為 ( ) ( ) ( ) j pm cp cp G j C j _e   (3-16) 將(3-14)和(3-16)式相除後 ( ) j pm cos sin cp c c pm c pm C j K e K   jK  (3-17)

由於控制器是分數階控制器，所以將控制器先進行整理 2 2 ( ) ( ) ( ) cos( ) cos( ) 2 2 sin( ) sin( ) 2 2 j j i cp p d cp p i cp d cp cp p i cp d cp i cp d cp K C j K K j K K e K e j K K K j K K _ _                                          _  _   (3-18) 把(3-17)和(3-18)式分為實部和虛部後再進行整理





cos cos cos 0

2 2

p i cp d cp c pm

K K _ _K  _ _K  

    (3-19)





sin sin sin 0

2 2 i cp d cp c pm K   K    K   _{ } _{ }      (3-20) 得出這兩式之後便可用來調整控制器的分數階數



 &



，藉此滿足相位邊界 的設計需求。方法整體的運行程序如下： Step.1 設定期望達到的相位邊界值



m Step.2 從Ziegler–Nichols 調整方法獲得

K

p_和

K

i

此篇提到模擬結果表示

K

d對(3-19)和(3-20)式求解參數的影響很大，因此可以微 調

K

d來尋找(3-19)和(3-20)式的最佳解。

3.4.2 Yeroglu 調整法的缺陷

此分數階 PID 控制器調整方法只透過一個點的移動來獲得所需的響應，也就是 透過調整  和



來獲得新的頻率響應，對於這個問題至少2 個點移動到新位置來 改善。 在這個方法中是假定系統擁有一些增益邊界，但是生活中有許多系統是不穩定 的以及沒有增益邊界，對於這些系統要套用到此方法是有難度的，所以後續的調 整方法也對於此問題進行改善。

3.5 改善對分數階 PID 控制器調整方法

上述調整方法主要的缺陷在於無法應用在不穩定系統，但是不穩定系統若可以 使用傳統PID 控制器穩定，就可以藉此調整到期望性能的分數階 PID 控制器。 在 Deepak & Kumari [62]文中提出的改善調整方法分為兩類，分別為針對穩定 系統以及不穩定系統進行探討。圖3-4 為運用在不穩定系統的示意圖，將兩曲線 各別取



G j

( )



 

180

 (A 和_A)和 G j(



) 1 (B 和_B)兩處的關係式，由得出 關係式調整控制器參數，使頻率響應曲線由實線部分移動到虛線部分。

在本篇中，針對不穩定系統的部分行修改，在不穩定系統中使用 Deepak & Kumari 調整方法，必須先找到一組傳統 PID 控制器使系統穩定後再調整，而在這 當中沒有明確地提出要用什麼方法找出這組控制器，因此本篇將原本所使用的傳 統PID 控制器變換為分數階 PID 控制器，如圖 3-5 所示，也就是將第一階段的設 計方法結合此處的調整方法進而得到期望性能的控制器，並且可以增加設計方法 的靈活性與自由度。 圖 3-5 改良調整法示意圖

由於這一階段控制器形式與前一階段不同，故必須進行參數變換的動作 ( ) i p d K C s K K s s      _(3-21) 其中

K

p



K

p0；

K

i



K K

p0 i0；

K

d



K K

p0 d0。 定義期望性能的 FOPID 控制器為 ( ) i p d K C s K K s s            _(3-22) A 點到A點為對應於系統的增益邊界，其公式的描述如下 Point A： ( ) ( ) j cp cp a a G j C j r e  r _(3-23) Point A： 1 1 ( ) ( ) j cp cp a a G j _ C j  _ r e   r _(3-24) 其中

r

_a和

r

_a1是頻率響應曲線在



G j

( )



 

180

處的大小。假設



cp





cp'，進行 以下的推導過程。

由前一階段得到的初始控制器頻域形式為









( ) i cp p d cp cp K C j K K j j         _(3-26) 和達到期望性能的控制器頻域形式為 ( ) ( ) ( ) i cp p d cp cp K C j K K j j               _(3-27) 將(3-25)、(3-26)和(3-27)經過實部和虛部的分類整理後可得出以下兩式，並且可 用來調整控制器。 1 cos₂ cos ₂ cos cos 0 2 2 a p i cp d cp a p i cp d cp r K K K r K K K                       _{ }       _{    }  _   _   (3-28) 1 sin sin 2 2 sin sin 0 2 2 a d cp i cp a d cp i cp r K K r K K                      _       _{   }  _  _   (3-29)

相同地，B 點到B點為對應於系統的相位邊界，公式的描述如下所示 Point B： 1 ( ) ( ) j pm cg cg G j C j e  (3-30) Point B： 2 ( ) ( ) j pm cg cg G j C j  e   _(3-31) 假設



cg





cg後，得出下式 1 2 ( ) ( ) ( ) ( ) pm pm j cg j cg C j e C j e           (3-32)

利用(3-26)、(3-27)和(3-32)，同樣地整理成實數部分和虛數部分後可得以下兩式

1

1

2

2

( cos cos ) cos

2 2

( sin sin ) sin

2 2

( cos cos ) cos

2 2

( sin sin )sin 0

2 2 p i cg d cg pm d cg i cg pm p i cg d cg pm d cg i cg pm K K K K K K K K K K                                                               (3-33) 1 1 2 2

( cos cos )sin

2 2

( sin sin ) cos

2 2

( cos cos )sin

2 2

( sin sin ) cos 0

2 2 p i cg d cg pm d cg i cg pm p i cg d cg pm d cg i cg pm K K K K K K K K K K                                                               (3-34) 在系統的頻率響應裡，透過 A 和 B 兩點的移動，吾人得到了四條方程式 [(3-28), (3-29), (3-33), (3-34)]，運用這組聯立方程式找出最佳參數

K

p、

K

i、

K

d、 和



，並設計出符合期望效能的FOPID 控制器。 此處吾人使用 Matlab 來求解[(3-28), (3-29), (3-33), (3-34)]_{這四條方程式，從} Matlab 內建的函數庫中選用 ’Fsolve’ 來編寫程式。使用此函數必須先給予

K

p、 i

K

、

K

d、和



的初值，

K

p、

K

i和

K

d初值給予的方法吾人選用Ziegler–Nichols 參數調整方法來設計，和



的初值就直接選用由第一階段PSO 得出的參數。

第四章、 例題模擬與討論

對無人飛行載具(unmanned aerial vehicle，UAV)做縱向飛行控制設計，本章節 將Wang 等人[60]的方法所設計出來的分數階 PI 控制器、傳統控制器和用本篇方 法所設計出來的分數階PID 控制器進行比較。在本章節中會討論使用 PSO 尋找初 始控制器所遇到的問題，並且針對問題點提出修改後的PSO。

4.1 UAV 飛行姿態控制模型

當飛機在穩定航行時，假設速率變化和仰角變化為 0，其俯仰角速度根據 Kada & Ghazzawi [66]文獻中的模型為 2 2 17.11 36.06 ( ) 3.981 26.47 s P s s s       (4-1)

舵機模型使用一階慣性環節(first-order inertia link)，其轉移函數如

 

1 0.1 1 G s s     (4-2) 然後俯仰率迴授系統的開路轉移函數為 3 2 171.1 360.6 ( ) 13.981 66.28 264.7 s P s s s s       (4-3) 再接上磁環外圈(outer ring)後的開路轉移函數為

4 3 2 171.1( ) 360.6 ( ) ( ) 13.981( ) 66.28( ) 264.7( ) j P j j j j j













     (4-5) 其相位和振幅分別為 3 4 2 171.1 336.82 13.981 [ ( )] arctan arctan 360.6 100.5 Arg P j            _{ } _{ }    _{ } (4-6) 2 2 4 2 2 3 2 360.6 (171.1 ) ( ) ( 100.5 ) (336.82 13.981 ) P j           (4-7) 將(4-6)和(4-7)式代入(2-5)和(2-6)式便可得出接上控制器後系統的相位和振幅 0 0 0 0 0 3 4 2 [ ( )] [ ( )] [ ( )] sin( / 2) sin( / 2) arctan 1 cos( / 2) cos( / 2) 171.1 336.82 13.981 arctan arctan 360.6 100.5 i d i d

Arg G j Arg C j Arg P j

K K K K                                    _{ } _{ }      (4-8) 0 2 0 0 0 1 2 2 0 0 2 2 4 2 2 3 2 ( ) ( ) ( ) (1 cos( / 2) cos( / 2)) ( sin( / 2) sin( / 2)) 360.6 (171.1 ) ( 100.5 ) (336.82 13.981 ) p i d i d G j C j P j K K K K K                          _      _      (4-9) 把(4-8)和(4-9)式套入章節 3.1 的設計規範(i)、(ii)和(iii)中，運用 PSO 找尋分數階 PID 控制器參數後就可以得出初始控制器。

4.2 使用 PSO 搜尋初始控制器

直接把章節 3.3.3 的線性 PSO 更新公式套用到章節 3.1 的設計規範(i)、(ii)和 (iii)，設計出來的控制器不僅效果不如預期，還會發生由 PSO 得出來的最佳粒子 中參數雖然完全符合設計規範，但是接上系統後的輸出響應與頻率響應所呈現的 效果極差，而且系統會不穩定。 為了因應此問題的發生，吾人將 PSO 做了一些改善，例如慣性權重的分配、粒 子優劣的條件設定等等。

4.2.1 慣性權重的改善

如果將慣性權重



直接套入更新公式，則會使控制器參數其中的分數階次  和



_{超出範圍，分數階次範圍通常設為0 到 2 之間的實數，原因在於使用線性更新} 方法會使_{PSO 前期粒子速度的變化較為巨大，控制器參數當中 和}



與

K

p、

K

i和 d

K

相較之下有範圍的限制，若是直接使用相同的慣性權重就會很容易發生超出範 圍的問題。 因此，在慣性權重的設計方面，在

K

p、

K

i和

K

d以及 和



這兩組參數採用不同 比例的方式進行粒子更新。比例的選用經過測試後，得出在對參數

K

p、

K

i和

K

d 以及 和



之間採用1 比 0.1 時效果顯著，也就是說對參數

K

p、

K

i和

K

d維持原本 的慣性權重，對  和



則採用原本慣性權重的0.1 倍，在這樣的調整之下能夠明

4.2.2 粒子的優劣

在 PSO 的過程中，每個粒子的適應值是利用目標函數計算出來的，從本篇的應 用來看的話，就是把章節3.1 的設計規範(i)、(ii)和(iii)做為目標函數，計算每個粒 子對目標函數的誤差，再選出最佳粒子。 但是這裡卻存在非常大的問題，如果只考慮粒子對目標函數的誤差，會讓有些 粒子對目標函數的誤差表面上看起來很低，但是透過觀察圖4-1 的系統開路頻率 響應曲線就可以發現增益曲線多次通過0 dB，發生這個現象的系統一般都是不穩 定系統，我們將這類粒子分類為不好的粒子。 圖 4-1 系統波德圖(不好的粒子)

從增益曲線可看到系統多次的通過 0 dB，由於這不是一個好的現象，所以必 須加入一些條件使PSO 在粒子優劣的選擇方面進行改善。 因此，吾人將在 PSO 過程中加入一道程序，如圖 4-2 所示，在計算粒子適應值 到更新粒子速度與位置的中間加入判斷條件，如果有粒子的參數會造成系統增益 曲線通過0 dB 兩次以上，吾人將會給予其粒子懲罰，也就是在做適應值的比較時， 不把此粒子一起進行比較，但並不代表此粒子就從PSO 的粒子群中剔除，在進行 粒子更新時此粒子也會一起更新，原因在於若是都把不好的粒子剔除，只留下優 良的粒子會使PSO 受到限制，不能往全域的方向進行搜尋。

4.3 例題模擬與系統響應

選用 UAV 飛行姿態控制模型(4-4)式當作受控系統，對此系統分別設計傳統控 制器、分數階PI 控制器和分數階 PID 控制器做比較，傳統控制器和分數階 PI 控 制器使用Wang 等人[60]所設計的控制器，分數階 PID 控制器則是使用本篇提出 的方法做設計。 這邊給予控制器的設計條件分別為，增益邊界

g

m



10dB

(

r

a1



0.31622

)，相位 邊界

  

m

75

，增益交越頻率

 

c

1.85 rad / s

。 使用傳統控制器設計的控制器為 1 ( ) 1.2 IO C s s   _(4-10) 搭配控制器章節 3.1 的設計規範所設計的分數階 PI 控制器為

 

0.4409 1.1495 0.5561 PI C s s   _(4-11) 使用改善後的 PSO 找出初始控制器之後，再將其代入章節 3.5 的調整方法，運 用Matlab 內建的函數 ’fslove’ 求解聯立方程[(3-28), (3-29), (3-33), (3-34)]，得出 分數階PID 控制器的最佳參數組合。 其控制器轉移函數為

(a) 0 到 30 秒反應時間

(b) 暫態響應放大圖

從圖 4-3 可以看出傳統控制器與分數階控制器在步階響應上有相當的差異，尤 其在最大超越量明顯大於其餘兩組分數階控制器。從放大圖之中比較兩組分數階 控制器的響應，吾人發現用本篇提出的方法設計出來的分數階PID 控制器比原例 題設計的PI 控制器更快達到穩定區間，而且最大超越量也比較低。

圖 4-4 系統開路波德圖(傳統控制器)

圖 4-4 當中的相位曲線圖可以很明顯地看出，此傳統控制器在設計時並沒有考 慮控制器章節3.1 的設計規範，故使曲線呈現如上圖的起伏現象，這也造成系統 並沒有強健性。

圖 4-5 系統開路波德圖(FOPI 和 FOPID 比較)

上圖為兩組分數階控制器的比較，由於這兩組控制器在設計時都有使用控制器 章節3.1 的設計規範，所以從相位曲線圖可以明顯地比較出與圖 4-4 的差異，在分 數階控制器圖中可以看出這兩組控制器都有賦予系統強健性。

第五章、 結論與未來研究方向

5.1 結論

本篇針對 Wang 等人[60]所提出的方法只能針對相位邊界、增益交越頻率和相 位曲線的變化量做控制器設計這點進行改善。由於Wang 等人[60]的方法無法充 分的運用分數階控制器有5 個參數可以調整的特性，因此本篇融入了 Deepak & Kumari [62]的方法，其方法的特性為可以透過控制器的 5 個參數同時進行調整， 將原本沒辦法加入增益邊界這點進行改善，並且也能夠保有章節3.1 中控制器設 計規範的特性，也就是能夠保持相位曲線在交越頻率的變化量為零的設計規範， 使設計出來的控制器附有強健性。 從單位步階輸出響應的最大超越量作比較，由範例所提供的控制器分別是，傳 統PID 控制器為 1.172，FOPI 控制器為 1.08；由本論文提出的方法所設計出的 FOPID 控制器為 1.075。傳統控制器的最大超越量明顯大於分數階控制器，而從 本論文方法所設計的FOPID 控制器亦能更進一步地降低原例題所提供的 FOPI 控 制器所產生的最大超越量。

5.2 未來展望

舉出一些未來能夠繼續研究的方向與相關的主題： 1. 將本論文的方法擴充到其他系統的應用，如時間延遲系統。 2. 使用其他參數調整方法替換 Ziegler–Nichols 參數調整方法使其在用 Matlab 求 解參數能有更好的效果，或者是使用其他方法取代Matlab 的函數 ’fsolve’ 來 尋找控制器參數。 3. 尋找其他更有效的方法將本篇使用兩個階段做設計這點改善至可以直接設計 出達到期望效能的控制器。

參考文獻

[1] G. W. Leibniz, Letter from Hanover, Germany, September 30, 1695 to G. A. L'Hospital. Leibnizen Mathematische Schrifien, Vol. 2, pp. 301-302, Olms Verlag, Hildesheim, Germany, 1962, First published in 1849.

[2] J. Liouville, " Memoire sur quelques Questions de Geometrie et de Mecanique, et sur un nouveau genre de Calcul pour resoudre ces Questions," J. Ecole Polytech, Section 21, pp. 1-69, 1832.

[3] B. Riemann, "Versuch einer allgemeinen Auffasung der Integration und Differen-Differentiation," The Collected Works of Bernhard Riemann (H. Weber, ed.), 2nd ed. Dover, New York, 1953.

[4] H. J. Holmgren, "Om differentialkalkylen med indices of hvad nature sam heist," Kongliga Svenska Vetenskaps-Akademiens Handlinger 11, pp. 1-83, 1864. [5] L. Euler, "De Progressionibus Transcentibus, sev Quarum Termini Algebraice

Dari Nequeunt," Comment. Acad. Sci. Imperialis Petropolitanae 5, pp. 38-57, 1738.

[6] J. L. Lagrange, " Sur une nouvelle espece de calcul relatif a la differentiation et a l'integration des quantites variables," Oeuvres de Lagrange, Vol. 3, pp. 441-476. Gauthier-Villars, Paris, 1849, First appeared in Nouv. Mem. Acad. Roy. Sci. Belles- Lett. Berlin 3, pp. 185-206, 1772.

[7] A. K. Griinwald, " Ueber ‘begrenzte’ Derivationen und deren Anwendung," Z. Math. Phys. 12, pp. 441-480, 1867.

[8] A. Krug, "Theorie der Derivationen," Akad. Wiss. Wien Denkenschriften, Math. Naturwiss. Kl. 57, pp. 151-228, 1890.

[9] H. Weyl, "Bemerkungen zum Begriff des Differentialquotienten gebrochener Ordnung," Vierteljschr. Naturforsch. Gesellsch. Zurich 62, pp. 296-302, 1971.

[10] G. H. Hardy, "On Some Properties of Integrals of Fractional Order," Messenger Math. 41, pp. 145-150, 1917.

[11] G. H. Hardy and J. E. Littlewood, "Some Properties of Fractional Integrals," Proc. London Math. Soc. [2], 24, pp. 37-41, 1925.

[12] G. H. Hardy and J. E. Littlewood, "Some Properties of Fractional Integrals, I," Math. Z. 27, pp. 565-606, 1928.

[13] G. H. Hardy and J. E. Littlewood, "Some Properties of Fractional Integrals, II," Math. Z. 34, pp. 403-439, 1932.

[14] H. Kober, " On Fractional Integrals and Derivatives," Quart. J. Math. Oxford Ser. 11, pp. 193-211, 1940.

[15] B. Kuttner, "Some Theorems on Fractional Derivatives," Proc. London Math. Soc. [3] 3, pp. 480-497, 1953.

[16] A. Erdélyi, "Transformation of Hypergeometric Integrals by Means of Fractional Integration by Parts," Quart. J. Math. Oxford Ser. 10, pp. 176-189, 1939.

[17] A. Erdélyi, "On fractional integration and its applications to the theory of Hankel transforms," Quart. J. Math. Oxford Ser. 11, pp. 293, 1940.

[18] A. Erdélyi and staff of the Bateman Manuscript Project, "Tables of Integral Trans- Transforms," Vol. 2, pp. 181-214. McGraw-Hill, New York, 1954.

[19] T. J. Osler, "Leibniz Rule for Fractional Derivatives Generalized and an Application to Infinite Series," SIAM J. Appl. Math. 16, pp. 658-674, 1970. [20] E. L. Post, "Generalized Differentiation," Trans. Amer. Math. Soc. 32, pp.

[23] A. Erdelyi, "Axially Symmetric Potentials and Fractional Integration," SIAM J. Appl. Math. 13, pp. 216-228, 1965.

[24] T. P. Higgins, "The Use of Fractional Integral Operators for Solving Nonhomo-geneous Differential Equations," Document DI-82-0677, Boeing Sci. Res. Lab., Seattle, Washington, 1967.

[25] G. W. Scott Blair, "The role of psychophysics in rheology," J. Colloid Sci. 2, pp. 21, 1947.

[26] T. D. Shermergor, "On the use of fractional differentiation operators for describing the hereditary properties of materials," Z. Prikl. Mech. i Tekhn. Fiz. 6, pp. 118, 1966.

[27] G. W. Scott Blair, "The role of psychophysics in rheology," J. Colloid Sci. 2, pp. 21, 1947.

[28] G. W. Scott Blair, "Measurements of Mind and Matter," Dennis Dobson, London, 1950.

[29] G. W. Scott Blair, " Some aspects of the search for invariants," British J. Philos. Sci. 1, pp. 230, 1950.

[30] G. W. Scott Blair, and J. E. Caffyn, "An application of the theory of quasi-properties to the treatment of anomalous stress-strain relations," Philos. Mag. 40, pp. 80, 1949.

[31] A. Graham, G. W. Scott Blair, and R. F. J. Withers, "A methodological problem in rheology," British J. Philos. Sci. 11, pp. 265, 1961.

[32] V. A. Belavin, R. Sh. Nigmatullin, A. I. Miroshnikov, and N. K. Lutskaya, "Fractional differentiation of oscillographic polarograms by means of an

[34] K. B. Oldham, and J. Spanier, "The replacement of Fick's laws by a formulation involving semidifferentiation," J. Electroanal. Chem. Interfacial Electrochem. 26, pp. 331, 1970.

[35] M. Grenness, and K. B. Oldham, " Semiintegral electroanalysis: theory and verification," Anal. Chem. 44, pp. 1121, 1972.

[36] R. L. Somorjai, and D. M. Bishop, "Integral-transformation trial functions of the fractional-integral class." Phys. Rev. A 1, pp. 1013, 1970.

[37] K. B. Oldham, "Diffusive transport to planar, cylindrical, and spherical electrodes," J. Electroanal. Chem. Interfacial Electrochem. 41, pp. 351, 1973. [38] K. B. Oldham, and J. Spanier, "A general solution of the diffusion equation for

semiinfinite geometries," J. Math. Anal. Appl. 39, pp. 655, 1972.

[39] I. Podlubny, "Fractional-order systems and PI D controllers," IEEE Transactions on Automatic Control, volume. 44, issue 1, pp 208-214, 1999. [40] M. Dorigo, "Optimization, Learning and Natural Algorithms," PhD Dissertation,

Politecnico di Milano, 1992.

[41] J. Kennedy, and R. Eberhart, "Particle Swarm Optimization," Proc. of the IEEE International Conference on Neural Networks, N. J. IEEE Press, pp. 1942-1948, 1995.

[42] Y. Shi, and R. Eberhart, "A Modified Particle Swarm Optimizer," Procedings of the 1998 IEEE International Conference on Evolutionary Computation, pp. 69-73, 1998.

[43] Y. Shi, and R. Eberhart, "Empirical Study of Particle Swarm Optimization," Proceedings of the 1999 Congress on Evolutionary Computation, Vol.3, pp.

[45] A. Ratnaweera, S. K. Halgamuge, and H. C. Watson, "Self-Organizing Hierarchical Particle Swarm Optimizer With Time-Varying Acceleration Coefficients," IEEE Transactions on Evolutionary Computation, Vol. 8(3), pp.240-255, 2004.

[46] A. Chatterjee, and P. Siarry, "Nonlinear inertia weight variation for dynamic adaptation in particle swarm optimization," Computers and Operations Research, Vol. 33, pp. 859-871, 2006.

[47] C. N. Ko, Y. P. Chang, and C. J. Wu, "An orthogonal-array-based particle swarm optimizer with nonlinear time-varying evolution," Applied Mathematics and Computation, Vol. 191, pp. 272-279, 2007.

[48] K. Oldham, and J. Spanier, "The Fractional Calculus," Hardcover: 234 pages. Publisher: Academic Press, pp. 45-60. 1974.

[49] M. D. Ortigueira, and F. Coito, "From Differences to Derivatives," Fractional Calculus and Applied Analysis. 7(4), pp. 459-471, 2004.

[50] M. Riesz, "L'intégrale de Riemann-Liouville et le problème de Cauchy," Acta Mathematica 81 (1), pp. 1-223, 1949.

[51] E. Hille and R. S. Phillips, "Functional analysis and semi-groups," Revised Edition, Amer. Math. Soc. 1957.

[52] K. S. Miller, and B. Ross, "An Introduction to the Fractional Calculus and Fractional Differential Equations," John Wiley & Sons, 1993.

[53] P.I. Lizorkin, "Fractional integration and differentiation in M. Hazewinkel, "Encyclopedia of Mathematics," Springer Verlag, Berlin, Heidelberg, New York, 2001.

[56] J. G. Ziegler and N. B. Nichols, "Optimum settings for automatic controllers," Trans. Amer. Soc. Mech. Eng., vol. 64, pp. 759–768, 1942.

[57] A. S. McCormack, and K. R. Godfrey, "Rule-Based Autotuning Based on Frequency Domain Identification," IEEE Transactions on Control Systems Technology, Vol. 6, No. 1, pp. 43-60, 1998.

[58] D. W. Pessen, "A new look at PID-controller tuning," Trans. Amer. Soc. Mech. Eng., J. Dynamic Syst., Meas., Contr., vol. 116, pp. 553–557, 1994.

[59] D. E. Seborg, T. F. Edgar, and D. A. Mellichamp, "Process Dynamics and Control." New York, Wiley, 1989.

[60] C. Y. Wang, J. F. Li, Y. D. Liu and M. Q. Li, "Research on UAV Attitude Control Based on Fractional Order Proportional Integral Controllers," IEEE Proceedings of the 33rd Chinese Control Conference, Nanjing, China, pp. 3552-3557, 2014.

[61] J. H. Chou, "Optimization approaches," Lecture Notes and Handouts at the Institute of Electrical Engineering (National Kaohsiung University of Applied Sciences), Taiwan, part 3, pp. 1-42, 2013.

[62] D. V. Deepak, and S. U. Kumari, "Modified method of tuning for fractional PID controllers," Power Signals Control and Computations (EPSCICON), 2014 International Conference on IEEE, pp. 1-6, 2014.

[63] C. Yeroglu, and N. Tan, "Note on fractional-order proportional integral-differential controller design," IET Control Theory Application, Vol. 5, issue 17, pp. 1978-1989, 2011.

[64] C. A. monje, Y. Q. Chen, D. Xue, B. M. Vinagre and V. Feliu, "Fractional-order Systems and Control : Fundamentals and Applications," Advances in Industrial Control, Springer, 2010.

[65] D. Xue, Y. Q. Chen and D. P. Atherton, "Linear Feedback Control Analysis and Design with MATLAB, " Advances in Design and Control, Siam, 2007.