粒子群演算法尋找比例控制器P最佳參數

為了驗證粒子群演算法在定位控制的應用之可行性，所有應用到粒子群演算法的相關參數、係數與作法同前述。

茲僅以比例控制器(P-controller)做驗證，至於比例-積分控制(Proportional-Integrate control, PI control)、比例-微分控制(Proportional-Differential control, PD control)有興趣的研究者可以參考。

茲選擇受控對象的轉移函數為:
y(s)/u(s)=(s+1)/(s^2+2s+1)...................(1)

此受控對象的步階響應(Step response)為下圖所示:

(圖一)

由圖一中可得知受控對象的響應，上升時間大約5秒(0~100% level)，安定時間約6秒。
換句話說，此受控對象加入一個穩定的電源，從開始到穩定，大概需要6秒的時間。

為了使受控對象的輸出響應可以獲得比較滿意的結果(例如:上升時間低於1秒，或安定時間能
縮短)，因此設計一個良好的控制系統是非常重要的。
PID控制器對於受控對象系統模型不明確仍有極佳的控制效能，因此廣泛使用於工業控制。然而PID雖然結構簡單，但如果參數設定不佳，仍有可能使系統產生不佳的控制效果，甚至造成系統的不穩定。

廣泛地從相關文獻可略而得知:
比例控制器(P-Controller):有效縮短響應時間、存在穩態誤差
積分控制器(I-Controller):增加系統的型式、改善響應的穩態誤差、但易使系統產生不穩定
微分控制器(D-Controller):改善響應時間、增加系統穩定能力

模擬實驗中，我們初期先撒出200個particles，每一個particle代表kp值(數值介於5~20之間)，其參數設定請參酌程式碼。每一個kp帶入副函式(PSO_position_control.m)，求得每一個kp所產生的響應結果，將適應函數的值(本例採error平方的加總結果)記錄之，並比較每一組的優劣，最終找出最佳的結果。

[程式碼]

clear all
clc

%初始化
N=200;         %partical 個數
Period=10;   % 迭代次數
Vmax=1;

kp=15*rand(1,N)+5;            %產生N個數的kp隨機值，範圍介於+5~+20
P_best=kp;                    %個體最優值，初始狀態取kp相同值
G_best=25;
V=zeros(1,N);
C1=2;
C2=2;
W=0.8;

for i=1:Period
i
for j=1:N
    V(2,j)=W*V(1,j)+C1*rand(1)*(P_best(1,j)-kp(1,j))+C2*rand(1)*(G_best-kp(1,j));
    if(abs(V(2,j)>Vmax))
        V(2,j)=sign(V(2,j))*Vmax;
    end
    kp(2,j)=kp(1,j)+V(2,j);
    P_best_fitness=PSO_position_control(P_best(1,j));
    current_particle_fitness=PSO_position_control(kp(2,j));
    if(current_particle_fitness<P_best_fitness)
        P_best(1,j)=kp(2,j);
    end
end
V(1,:)=V(2,:);
kp(1,:)=kp(2,:);
for j=1:N
    P_best_fitness=PSO_position_control(P_best(1,j));
    G_best_fitness=PSO_position_control(G_best);
    if(P_best_fitness<G_best_fitness)
        G_best=P_best(1,j);
    end
end
end

G_best

[副函式PSO_position_control.m]

function fitness=PSO_position_control(kp)
%系統轉移函數  sys_open(s)=(s+1)/(s^2+2s+1)
num=[1 1];
den=[1 2 1];
sys_open=tf(num,den);

%轉移函數轉換為狀態空間方程式
[a,b,c,d]=tf2ss(num,den);
sys_open=ss(a,b,c,d);

%初始設定
Ts=0.01;       %sampling time
x=zeros(2,1);  %state space variables is set to zero
u=0;
R=10;          %reference value
y=0;           %current state space output value, assume as velocity
position=0;    %current position
cnt=1;
err=0;         %error between reference position value and current one.
sum_err=0;

for t=0:Ts:10
 dot_x=a*x+b*u;
 y=c*x+d*u;

 position=position+y*Ts;

 err=R-position;
 sum_err=sum_err+0.5*(err)^2;
 u=kp*err;

 x=dot_x;
 plot_position(cnt)=position;
 cnt=cnt+1;
end

fitness=sum_err;

[執行結果]
經由最佳化搜尋，找到kp=29.0621可有良好的效能，因此將此參數帶入觀察，其響應結果如下圖所示：

（圖二）

其上升時間、安定時間都在1秒之內，效果卓著。

求解最佳化的問題-使用粒子群演算法,part4

最佳化問題，求解一維函數F(x)=x^2的極值

粒子群演算法相關參數初值定義如下：

particle 個數N=200

迭代次數period=1000次

Vmax=5粒子飛行速度最大值

X=10*rand(1,N)-5，隨機產生N個-5~+5之間的數

V＝zeros(1,N)，所有粒子飛行速度初值均為0

C1=C2=2

W=0.8

G_best=1000，寫大一點，否則會一下子就收斂

P_best=X ，一開始，個體的最佳值即一開始當前的位置

[程式碼]

%一維函數F(x)=x^2，給定範圍找最小值

clear all
clc

%初始化
N=200;         %partical 個數
Period=1000;   % 迭代次數
Vmax=5;

X=10*rand(1,N)-5;            %產生N個數的隨機值，範圍介於-5~+5
P_best=X;                    %個體最優值，初始狀態取X相同值
G_best=1000;
V=zeros(1,N);
C1=2;
C2=2;
W=0.8;

cnt=1;
for i=1:Period
    for j=1:N
        V(2,j)=W*V(1,j)+C1*rand(1)*(P_best(1,j)-X(1,j))+C2*rand(1)*(G_best-X(1,j));
        if(abs(V(2,j)>Vmax))
            V(2,j)=sign(V(2,j))*Vmax;
        end
        X(2,j)=X(1,j)+V(2,j);
        P_best_fitness=fitness(P_best(1,j));
        current_particle_fitness=fitness(X(2,j));
        if(current_particle_fitness<P_best_fitness)
            P_best(1,j)=X(2,j);
        end
    end
    V(1,:)=V(2,:);
    X(1,:)=X(2,:);
    for j=1:N
        P_best_fitness=fitness(P_best(1,j));
        G_best_fitness=fitness(G_best);
        if(P_best_fitness<G_best_fitness)
            G_best=P_best(1,j);
        end
    end
    G_best_dot(cnt)=G_best;
    cnt=cnt+1;
end

plot(G_best_dot)
G_best

副函數fitness:
function fit=fitness(X)
    fit=X^2;
end
 
[執行結果]
從圖中可看出最佳值收斂，且經計算得到G_best=-3.9680e-021
換句話說，電腦分析當-3.9680e-021可得函數F(-3.9680e-021)的最小值。當然，迭代的分析，必然存在誤差，但其結果仍令人滿意，與期望值0極為接近。

求解最佳化的問題-使用粒子群演算法,part3

接續前一篇討論…

粒子群演算法，在運算上非常的簡單，首先先初始撒下粒子之個數、粒子飛行速度。

演算法：

1.每一粒子飛行速度與方向遷就於粒子本身及群體，定義飛行速度：

Vi(k+1)=W*Vi(k)+C1*rand()*(P_best_i-Xi(k))+C2*rand()*(G_best-Xi(k))............(1)

其中:

W是權重常數，考慮前一刻飛行速度，對下一刻的飛行速度的影響程度

Vi(k+1)：第i個粒子，第k+1時刻的飛行速度

Vi(k)：第i個粒子，第k時刻的飛行速度

C1,C2：學習係數，當前位置與最佳位置，影響下一時刻飛行速度的權重

rand()：隨機函數，產生0~1之間的小數

P_best_i：第i個粒子，所走路徑中自始至今之最佳值

G_best：群體粒子之最佳值

Xi(k)：第i個粒子，當前所在位置

2.定義粒子下一時刻的位置：

Xi(k+1)=Xi(k)+Vi(k+1)...........................(2)

其中：

Xi(k+1)：第i個粒子，在第k+1時刻的位置

Xi(k)：第i個粒子，在第k時刻的位置

(2)式，速度‧時間＝Vi(k+1)△t=位置

3.第i粒子最佳位置取代原則：

需定義適應函數(fitness function)，若fit(Xi(k+1))優於fit(P_best_i)，則P_best_i=Xi(k+1)

4.群體最佳值的修正時機：

若fit(P_best_i)優於fit(G_best)，則取代當前G_best地位

求解最佳化的問題-使用粒子群演算法,part2

粒子群演算法仿效鳥類飛行追求卓越的生命歷程，很像「雁行理論」。

粒子群演算法(Particle Swarm Optimization, PSO)，其精神為「一群隨機粒子分佈於一個區域範圍，每一個體粒子搜尋鄰近範圍的最優，個體與個體之間相互比較（取適合度fitness function）並擇其優，經過反覆迭代的過程，找出總體的最佳值。」

以前述的例子來解釋，假設一個一維函數F(x)=x^2..............(1)

首先，先散佈一群(N個)隨機的粒子，每一個粒子都代表此函數的X,例如X1,X2,X3....,XN

每一個X帶入(1)都會有一個相對應的函數值，我們的目的是:

找出X等於多少時，函數有最小值？

在每一次迭代的過程，個體(particle)在一個環境探索(exploit) 的歷程中，反覆與舊有的經驗（個體自始至今的最佳）比較，如果當前的位置優於舊有的經驗，最佳位置則以新值取代。

簡單的說，一隻鳥覓食，假設一個環境，甲處有一隻蟲、乙處有兩隻蟲、丙處有五隻蟲，以鳥的觀點來說，這三處最好的位置，應該是丙處。而鳥在飛翔的過程，假設會先經過甲、再到乙，最終到達丙處，鳥在觀察的過程中會發現乙優於甲、丙又優於乙，所以最佳的位置會不斷被取代。

而一群個體，在一個環境探索的歷程，每一個個體都會找到當前最好的位置，但以一個群體來說，整體最優的位置在哪？

因此，在粒子群演算法中討論到，如果一群個體相互比較，發現某一個體找到的位置是群體中所有的最佳，那麼所有的個體，將追隨這個總體最佳的位置前進。

簡單地說，如果一群鳥在一個環境中，每一隻鳥都會找到不錯的覓食點，假設其中有一隻鳥挖到寶，發現在某一個地方有無限量供應的蟲（酷斃了！），那麼所有的鳥將追隨這個位置並往這個方向飛，而放棄他自己找到還不錯的覓食點。

所以，當你有空停下腳步觀看天上飛翔的鳥，你會發現：一群鳥會跟著其中一隻往某個方向飛…

（待續）

求解最佳化的問題-使用粒子群演算法,part1

研究所時，跟隨恩師學習各種控制技術，印象中最為深刻地莫過於探討最佳化的問題。

有關最佳化的研究，有許多相關文獻運用不同的技術或方法求解最佳化，而最終的目的無非是獲得最佳的結果。至於最佳化的目的，主要為克服人為的經驗誤差，因此求助於微電腦控制，找尋參數的最佳值。

在暴力搜尋法中，遺傳演算法(Genetic Algorithm) 是最經典的技術，另外，仿生物型態的，還包含有蟻行演算法（Ants Algorithm）、及粒子群演算法（Particle Swarm Optimization,仿鳥類生態）。

再談最佳化的問題前，我們先瞭解一個最簡單的問題：

已知一個一維變數的函數F(x)=X^2 (X平方)..............(1)

若X的範圍從-5到+5，我們如果用畫圖的方式，可得以下的圖形：

從圖可以看出，此函數為拋物線函數，且最小值為F(0)=0，即：當X=0時有最小值0

當然，如果用微分來計算，我們可知：

dF(x)/dx=2x....................(2)

令dF(x)/dx=0 可得極值，由（2）可知2x=0，所以x=0............(3)

把(3)的結果帶回(1)，所以F(x)=0

數位I/O與邏輯分析器的常用名詞及定義

數位I/O與邏輯分析器的常用名詞及定義
http://zone.ni.com/devzone/cda/tut/p/id/6247