買書捐殘盟

2011年7月28日 星期四

粒子群演算法尋找比例控制器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秒之內,效果卓著。

2011年7月25日 星期一

求解最佳化的問題-使用粒子群演算法,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

2011年7月15日 星期五