買書捐殘盟

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秒之內,效果卓著。

13 則留言:

  1. 搜尋參數最佳化,搜尋花費時間挺久的,所以程式碼我只做10次迭代,原本想做1000次,可能要跑好幾個小時,所以如果有興趣的人,想迭代比較精緻一點的,也許可以把迭代次數調高一點。

    回覆刪除
  2. 薛老師您好,

    謝謝你分享這麼棒的知識,原來雁行理論也可以運用在最佳化,雖然以前有學過控制理論,但是出社會後就一直停留在二十年前的舊知識,從老師深入淺出的說明,學會一種新方法。

    如果希望縮短計算時間,可以降低 particles ,個人認為收斂狀況 = particles x loops,也就是 200 particles x 10 loops = 10 particles x 200 loops。

    回覆刪除
  3. Dear Bridan:
    particles 數量少,可減少計算時間,卻影響精度,所以無法顧全。但優化的理論,確實有助於改善人為的經驗誤差。

    回覆刪除
  4. particles數量多或少會有甚麼結果?我把它比喻成尋找迷途的羔羊。
    一個農夫走失了一隻羊,在一個開放、未知或偌大的空間,假設只有這個農夫在尋找,找到的機率,相較於10個人一起找,我相信10個人尋羊的速度會比一個人快。如果是一群更多的人(超過10個),那麼效果更加顯著。
    所以電線桿如果看到有人張貼尋狗啟事,無非是希望藉由眾人的力量來協尋,比較有機會找到。

    回覆刪除
  5. 薛老師您好,

    如果目標物是不動的,兩百個人找十天,跟十個人找兩百天, 效果應該是一樣的。但是目標物是會動 的,那越多人一起找是比較好的,不然走失的羊, 拖得越久羊入虎口的機會越高,凶多吉少。

    可能韌體寫多了,對系統資源限制很在意,記憶體少少的系統, 我就只能派十個人去找。^_^

    回覆刪除
  6. PSO會影響執行速度,除了粒子個數,其次便是迭代迴圈,可再探討。

    回覆刪除
  7. 薛老師:
    你好
    小弟最近正在接觸pso相關的演算法
    想要請問pso演算法的"限制條件"的使用方向

    限制條件的副程式應該放在哪邊才好
    如果沒有達到限制條件
    是要重新rand?
    還是有別的方向?

    謝謝!!

    回覆刪除
  8. particles 如果有boundry value,例如我的x的邊界值在[-a,a],超過左右邊界,建議試調一下飛行速度V,可以改小一點,有時會因為飛太快,找到的值超出邊界。

    如果超出邊界,怎辦?你可以選擇放棄該點,那一點也許是更好的fitness沒錯,但既然超出範圍,就不是我們要找的目標對象。你可以在程式增加IF條件去設定他。

    回覆刪除
  9. 你好, 請問不知系統轉移函數的時候, 是否可以用 STEP INPUT 的響應曲線, 套用這個 PSO 的方法呢.

    另外, 需要考慮 DEAD TIME 嗎

    回覆刪除
  10. xiaolaba你好:
    關於你的問題,我的回覆是:
    如果系統數學模型未知,一般便是利用step response或是pulse response觀察它的output response,以便求其近似的數學模型。

    PID控制器對於未知系統模型的控制,有極佳的控制效能,在諸多文獻已有證實。本例中探討的,便是利用pso去找最佳的p,i,d參數來達最佳的控制效能
    (performance)。
    謝謝你的指教!歡迎一起討論!

    回覆刪除
  11. 作者已經移除這則留言。

    回覆刪除
  12. 您好 請問一下如果要找pid的其他參數 也是先寫出個別的副函式 然後直接套入主程式碼下去計算就對了嗎?
    還有一點就是sum_err=sum_err+0.5*(err)^2; 這段大概是什麼意思呢?

    真的很不好意思 希望麻煩大師幫我解惑

    回覆刪除
  13. 最適化參數,將其建構為副函式,取其最佳化參數後再帶回主程式去執行。
    sum_err=sum_err+0.5*(err)^2;
    是為了每次將誤差反饋,至於乘以0.5,你也可以乘以一個定值(0~1)。
    err^2是取誤差的數值永為正,你也可以取絕對值,在誤差反饋都可以接受。

    回覆刪除