遺伝的アルゴリズムは有名ですが、実数を扱うPID制御のパラメータ調整では、ビット数を実数に変えるテクニックが必要。
そこで、直接実数を使える実数値遺伝的アルゴリズムを使って、パラメータ調整をしてみました。
コードは以下の通り。
#遺伝的アルゴリズム from control.matlab import * import numpy as np import random import math import warnings import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D %matplotlib inline warnings.simplefilter('ignore')#警告を無視する #制御対象のパラメータ(2慣性、剛体モード+2次までの振動モードとする) kt=10 #剛体モード k0=100 #1次振動モード k1=100 zeta1=0.2 omega1=30*2*3.14 #2次振動モード k2=100 zeta2=0.04 omega2=200*2*3.14 #むだ時間(パデ近似で表現) deadtime_num,deadtime_den=pade(500e-6,4) deadtime_sys=tf(deadtime_num,deadtime_den) #制御対象の伝達関数 P=kt*(tf(k0,[1,0,0])+tf(k1,[1,2*zeta1*omega1,omega1**2])+tf(k2,[1,2*zeta2*omega2,omega2**2]))*deadtime_sys #print(P) fmin=1 fmax=10000 wrange=linspace(fmin*2*3.14,fmax*2*3.14,4000) bode(P,wrange) plt.show() #遺伝的アルゴリズムパラメータ generation=30#世代数 indivisual_num=1000#集団数 crossover_probability=0.3#交叉確率 min_para=0.1 max_para=1000 FB_para=np.zeros((indivisual_num,3))#個体 ISE_min=1e12#評価値 ISE_list=[]#世代ごとの評価値 #BLX_a法による交叉 def BLX_a(patents): alpha=0.5#拡張倍率 gene_num=3#遺伝子数 minP=[0,0,0] maxP=[0,0,0] child=[0,0,0] #親個体をランダム抽出 patent1=patents[random.randint(0,indivisual_num-1)] patent2=patents[random.randint(0,indivisual_num-1)] #子個体の作成 for gene in range(gene_num): minP[gene]=min(patent1[gene],patent2[gene])-alpha*abs(patent1[gene]-patent2[gene]) maxP[gene]=max(patent1[gene],patent2[gene])+alpha*abs(patent1[gene]-patent2[gene]) if minP[gene]<min_para: minP[gene]=min_para if maxP[gene]>max_para: maxP[gene]=max_para child[gene]=minP[gene]+random.random()*(maxP[gene]-minP[gene]) return child for i in range(generation): print("世代{}".format(i)) #開ループ特性による安定性判別 #初期はランダム設定 if i==0: for g in range(indivisual_num): kp=random.uniform(1,1000) kd=random.uniform(0.1,10) ki=random.uniform(10,1000) FB_para[g]=[kp,kd,ki] #以降は交叉もしくはエリート解 else: for g in range(indivisual_num): if random.random()>=crossover_probability: FB_para[g]=BLX_a(FB_para_patent) else: FB_para[g]=FB_para_opt #各個体についてシミュレーション for g in range(indivisual_num): #PID制御器 kp=FB_para[g][0] kd=FB_para[g][1] ki=FB_para[g][2] K=tf([kd,kp,ki],[1,0]) #開ループ特性 Popen=K*P #PID制御を用いたステップ応答 ref=1 Gyr=feedback(P*K,1)#FB制御系の設計 y_FB,t_FB=step(Gyr,np.arange(0,0.1,0.0001))#ステップ応答 ISE=0 #評価値計算 for j in range(int(0.01/0.0001),len(y_FB)): ISE+=(ref-y_FB[j])**2 #print(ISE) #エリート解との比較 if ISE<ISE_min: ISE_min=ISE kp_opt=kp kd_opt=kd ki_opt=ki FB_para_opt=FB_para[g] fig = plt.figure() ax = Axes3D(fig) ax.set_xlabel('kp') ax.set_xlim(0,1000) ax.set_ylabel('kd') ax.set_ylim(0,10) ax.set_zlabel('ki') ax.set_zlim(0,1000) plt.plot(FB_para[:,0],FB_para[:,1],FB_para[:,2],marker="o",linestyle='None') plt.show() print(ISE_min) ISE_list.append(ISE_min) FB_para_patent=FB_para print("kp={},ki={},kd={}".format(kp_opt,ki_opt,kd_opt)) #PID制御器 K=tf([kd_opt,kp_opt,ki_opt],[1,0]) #開ループ特性 Popen=K*P bode(Popen,wrange) plt.show() #安定余裕の確認 #print(margin(Popen))#2番目が位相余裕、3番目がゲイン余裕 print('ゲイン余裕:'+str(margin(Popen)[2])) print('位相余裕:'+str(margin(Popen)[1])) #ナイキスト線図 nyquist_range=linspace(30*2*3.14,500*2*3.14,1000)#図を見やすくするため開始周波数を調整 nyquist(Popen,omega=nyquist_range) plt.show() #感度関数 Psens=1/(1+K*P) bode(Psens,wrange) plt.show() #PID制御を用いたステップ応答 ref=1 Gyr=feedback(P*K,1)#FB制御系の設計 y_FB,t_FB=step(Gyr,np.arange(0,0.1,0.0001))#ステップ応答 #ステップ応答波形 fig,ax=plt.subplots() ax.plot(t_FB,y_FB*ref) ax.axhline(ref,color='k',linewidth=0.5) #評価値推移 plt.figure() plt.plot(range(0,generation),ISE_list)
交叉手法については、一番簡単なBLX-α法を使用。
用いる個体の選別方法については、ある一定の確率で交叉した子個体を使用、それ以外ではエリート会を使用する仕組みになっています。
なので、やや初期エリート解の値に左右されがちなところはあります。
まぁ交叉確率を調整すれば多様性は確保されますが、個体の収束速度とトレードオフになるので、そこらへんは試行錯誤で調整していく必要があります。
0~1000×3次元・1000個体の集団で回せば、およそ20世代(数分かかる)くらいで評価値が収束する感じですかね。
位置応答(ステップ応答)も、それなりに安定した応答になりました。
今回は短いですがこのへんで。
コメント