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