Разработка и исследование метода грамматической эволюции для структурно-параметрического синтеза системы управления динамическим объектом

def main():

n=6

LenP=40

LenPPL=72

NIter=32

start_t=time.time()

PPL=InitialPop(LenP, LenPPL, n)

p,func=ga(PPL, Criterion, NIter, n, LenP, LenPPL)

stop_t=time.time()-start_t

print p, func

print stop_t

x1,x2,u,t=ForGraphic(p[:-2], p[-2], p[-1])

PhaseTrajectory(x1,x2)

GraphicU(u,t)

if __name__=="__main__":

#cProfile.run('main()')

main()

Листинг модуля GenetikAlg:

from ParetoDomination import Pareto

from random import random,randint

from math import ceil, sqrt

###########################################

def bin2dec(b): # переводит двоичный код в десятичное число

s=''

for i in b: s+=str(i)

return int(s,2)

def merger(ppl,lenp,lenppl): # объединяет подмножества в хромосоме

for i in xrange(lenppl):

x=[]

for j in xrange(lenp+2): x+=ppl[i][j]

ppl[i]=x

return ppl

def dem(ch,n,lenp): # разделяет хромосому на подмножества

ii=0

chi=[0]*(lenp+2)

while (ii+n)<=len(ch):

for i in xrange(lenp+2):

b=ch[ii:ii+n]

ii+=n

chi[i]=b

return chi

def InitialPop(lenp,lenppl,n): # гнерирует начальную популяцию

ppl=[[[randint(0,1) for ii in xrange(n)]

for j in xrange(lenp+2)] for i in xrange(lenppl)]

return ppl

def Tournament(ppl,fit,lenppl): # турнирный отбор

parents=[0]*lenppl

for i in xrange(lenppl):

i1=randint(0,lenppl-1)

i2=randint(0,lenppl-1)

if i1==i2:

while i1==i2: i2=randint(0,lenppl-1)

win=max(fit[i1],fit[i2])

parents[i]=ppl[fit.index(win)]

return parents

def Selection(ppl,fit,lenppl): # отбор родителей со случайной вероятностью

parents=[0]*lenppl

for i in xrange(0,lenppl-1,2):

p = random()

i1 = randint(0,lenppl-1)

i2 = randint(0,lenppl-1)

while fit[i1]<p: i1 = randint(0,lenppl-1)

while fit[i2]<p: i2 = randint(0,lenppl-1)

parents[i]=ppl[i1]

parents[i+1]=ppl[i2]

return parents

def Crossover(parent,lenppl,ch,m,CrossType): # скрещивание

i1=ch

i2=ch+1

if CrossType=='Single-point':

n=randint(1,len(parent[i1])-(2+2*m))

child1=parent[i1][:n]+parent[i2][n:-2*m]+parent[i1][-2*m:]

child2=parent[i2][:n]+parent[i1][n:-2*m]+parent[i2][-2*m:]

if CrossType=='Two-points':

lenpar1=len(parent[i1])

n=randint(1,lenpar1-2)

n1=randint(1,lenpar1-2)

if n==n1:

while n==n1: n1=randint(1,lenppl-2)

if n<n1:

child1=parent[i1][:n]+parent[i2][n:n1]+parent[i1][n1:]

child2=parent[i2][:n]+parent[i1][n:n1]+parent[i2][n1:]

if n1<n:

child1=parent[i1][:n1]+parent[i2][n1:n]+parent[i1][n:]

child2=parent[i2][:n1]+parent[i1][n1:n]+parent[i2][n:]

return child1,child2

def Mutation(genom): # мутация

for i in xrange(len(genom)):

P=random()

if P<=1.0/(3.0+i):

if genom[i]==1: genom[i]=0

else: genom[i]=1

return genom

def ga(P,fun,Nitr,n,lenp,lenppl): # генетический алгоритм

itr=0

obj=Pareto()

func=[fun([bin2dec(P[j][i]) for i in xrange(lenp)],

bin2dec(P[j][-2]) , bin2dec(P[j][-1]) ) for j in xrange(lenppl)]

rank=obj.assign_ranks(func)

fit=obj.set_fitness_from_ranks(rank,lenppl)

while itr<Nitr:

itr+=1

print itr

par=Selection(P,fit,lenppl)

par=merger(par,lenp,lenppl)

Q=[]

for i in xrange(lenppl/2):

ch1,ch2=Crossover(par,lenppl,i*2,n,CrossType='Single-point')

ch1=Mutation(ch1)

ch2=Mutation(ch2)

Q.append(dem(ch1,n,lenp))

Q.append(dem(ch2,n,lenp))

PQ=P+Q

f2=[fun([bin2dec(PQ[j][i]) for i in xrange(lenp)],

bin2dec(PQ[j][-2]), bin2dec(PQ[j][-1])) for j in xrange(2*lenppl)]

rank2=obj.assign_ranks(f2)

fit=obj.set_fitness_from_ranks(rank2,2*lenppl)

fitz=zip(fit,range(2*lenppl))

fitz.sort()

fit=[fitz[i][0] for i in xrange(2*lenppl)]

P=[PQ[fitz[i][1]] for i in xrange(2*lenppl)]

fit=fit[lenppl:]

P=P[lenppl:]

best=[bin2dec(i) for i in P[-1]]

q=f2[fitz[-1][1]][:]

return best,q

Листинг модуля grammer_evolution:

from math import *

def gramevo(p):

expr=['(<expr>)<op>(<expr>)','<fun1>(<expr>)',

'<fun2>(<expr>,<expr>)','<var>']

fun1=['sqr','sign','atan','minus','sqrtX','sin','step','plfun1','plfun2']

fun2=['divX']

op=['+','-','*']

var=['xk[0]','xk[1]','q1','q2']

fun='<expr>'

bool1=1

v=0

while v!=-1:

for i in range(len(p)):

v=fun.find('<')

if fun[v+1:v+5]=='expr':

if bool1==1:

m=p[i] % len(expr)

fun = fun[:v] + expr[m] + fun[v+6:]

else:

fun = fun[:v] + '<var>' + fun[v+6:]

elif fun[v+1]=='f' and fun[v+4]=='1':

m=p[i] % len(fun1)

fun = fun[:v] + fun1[m] + fun[v+6:]

elif fun[v+1]=='f' and fun[v+4]=='2':

m=p[i] % len(fun2)

fun = fun[:v] + fun2[m] + fun[v+6:]

elif fun[v+1]=='o':

m=p[i] % len(op)

fun = fun[:v] + op[m] + fun[v+4:]

elif fun[v+1]=='v':

m=p[i] % len(var)

fun = fun[:v] + var[m] + fun[v+5:]

elif v==-1:

break

bool1=0

return fun

Листинг модуля Prod_laws:

from math import sqrt, exp

def minus(x):

return -(x)

def sqrtX(x):

return sign(x)*sqrt(abs(x))

def sqr(x):

return x**2

def sign(x):

if x>0 or x==0: result=1

else: result=-1

return result

def step(x):

if x>0 or x==0: y=1

else: y=0

return y

def plfun1(x):

try:

result=round(1/(1+exp(-x)),6)

except:

if sign(x)==-1:

result=0

elif sign(x)==1:

result=1

return result

def plfun2(x):

try:

result=round((1-exp(-x))/(1+exp(-x)),6)

except:

if sign(x)==-1:

result=0

elif sign(x)==1:

result=1

return result

def divX(x, y):

try:

z = x / y

except:

z = sign(x) * sign(y) * 10e+8

return z

Листинг модуля RK4:

def rk4(func,dt,time,x,u):

#time=0

#x = x1[:]

n=len(x[0])

f=func(time,x[-1],u)

x[-1]

k1=[dt*f[i] for i in xrange(n)]

f=func(time+(dt/2),[x[-1][i]+(k1[i]/2) for i in xrange(n)],u)

k2=[dt*f[i] for i in xrange(n)]

f=func(time+(dt/2),[x[-1][i]+(k2[i]/2) for i in xrange(n)],u)

k3=[dt*f[i] for i in xrange(n)]

f=func(time+dt,[x[-1][i]+k3[i] for i in xrange(n)],u)

k4=[dt*f[i] for i in xrange(n)]

return [x[-1][i]+(1.0/6.0)*(k1[i]+2*k2[i]+2*k3[i]+k4[i]) for i in xrange(n)]

Листинг модуля ParetoDomination:

class Pareto:

def dominate(self, F1, F2):

d = False

for i in range(len(F1[1])):

if F1[1][i] > F2[1][i]:

d = True

elif F2[1][i] > F1[1][i]:

return False

return d

def non_dominated_front(self, Group):

Front = []

for G in Group:

Front.append(G)

tmpFront = Front[:]

for C in Front:

if G != C:

if self.dominate(G, C):

tmpFront.remove(G)

break

elif self.dominate(C, G):

tmpFront.remove(C)

Front = tmpFront[:]

return Front

def assign_ranks(self, functionals):

tmp_functionals = [(i, functionals[i]) for i in range(len(functionals))]

ranks = []

i = 0

while tmp_functionals != []:

ranks.append(self.non_dominated_front(tmp_functionals))

for R in ranks[i]:

tmp_functionals.remove(R)

i += 1

# print ranks

return ranks

def set_fitness_from_ranks(self, ranks, pop_size):

Fitness = [0.0 for i in range(pop_size)]

cr = 0.0

for R in ranks:

for e in R:

Fitness[e[0]] = 1.0 / (1.0 + cr)

cr += 1

return Fitness

Листинг модуля ParetoDomination:

from grammar_evolution import gramevo

from prod_laws import *

from RK4 import rk4

from Euler2method import Euler2

from math import *

import matplotlib.pyplot as plt

from operator import *

#############################################

LenSys=2

x0=[1,1]

xf=[0.0,0.0]

tf=10.0

dt=0.01

#############################################

def Sys1(time,x,u):

y = [x[1], u]

for i in xrange(2):

if y[i]>=1e8: y[i]=sign(y[i])*1e8

return y

def Duffing(time,x,u):

a=1

B=1

y=[x[1],u-B*x[0]-(a*(x[0]**3))]

for i in xrange(2):

if y[i]>=1e8: y[i]=sign(y[i])*1e8

return y

def norma(x,xf):

s=0.0

for i in xrange(len(xf)): s+=(xf[i]-x[i])**2

return sqrt(s)

def Criterion(p, qx1, qx2):

q1=qx1

q2=qx2

x=[[x0[i] for i in xrange(LenSys)]]

t=0

miss_err=0

flag=False

Ft = tf

Fm = 1e+8

while (t<=tf):

xk = x[-1][:]

try:

u=eval(gramevo(p))

if u>1.0: u=1.0

if u<-1.0: u=-1.0

except:

miss_err=1e+8

print xk, q1, q2

print gramevo(p)

print eval(gramevo(p))

u = 1.0

t+=dt

x.append(rk4(Duffing,dt,t,x,u))

s=norma(x[-1],xf)

if s<=1e-3:

if flag == False:

Ft = t

flag = True

else:

Ft = t

flag = False

Fm = s

return [Ft, Fm+miss_err]

def ForGraphic(p, qx1, qx2):

q1=qx1

q2=qx2

x1=[]

x2=[]

ug=[]

tg=[]

x=[[x0[i] for i in xrange(LenSys)]]

t=0

while (t<=tf):

xk = x[-1][:]

u=eval(gramevo(p))

if u>1.0: u=1.0

if u<-1.0: u=-1.0

t+=dt

x.append(rk4(Duffing,dt,t,x,u))

# s=norma(x[-1],xf)

ug+=[u]

tg+=[t]

for e in x:

x1.append(e[0])

x2.append(e[1])

return x1,x2,ug,tg

def PhaseTrajectory(x1,x2):

plt.plot(x1,x2,label='X2(X1)')

plt.title('Phase trajectory')

plt.xlabel('X1')

plt.ylabel('X2')

plt.legend(loc='best')

plt.grid(True)

plt.show()

def GraphicU(u,t):

plt.plot(t,u,label='U(t)')

plt.title('Control')

plt.xlabel('t')

plt.ylabel('U')

plt.legend(loc='best')

plt.grid(True)

plt.show()

def GraphicX1(x1,t):

plt.plot(t,x1,label='X1(t)')

plt.title('Graphic X1')

plt.xlabel('t')

plt.ylabel('X1')

plt.legend(loc='best')

plt.grid(True)

plt.show()

def GraphicX2(x2,t):

plt.plot(t,x2,label='X2(t)')

plt.title('Graphic X2')

plt.xlabel('t')

plt.ylabel('X2')

plt.legend(loc='best')

plt.grid(True)

plt.show()

Размещено на Allbest.ru