遺伝的アルゴリズムは有名ですが、実数を扱う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世代(数分かかる)くらいで評価値が収束する感じですかね。

GA評価値

位置応答(ステップ応答)も、それなりに安定した応答になりました。

GA_ステップ




今回は短いですがこのへんで。
 
このエントリーをはてなブックマークに追加

コメント

コメントフォーム
評価する
  • 1
  • 2
  • 3
  • 4
  • 5
  • リセット
  • 1
  • 2
  • 3
  • 4
  • 5
  • リセット