In [1]:
import random as rnd
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import pandas as pd
from pandas import DataFrame
%pylab inline
In [31]:
# Gera o nome do conjunto
def conj_nome(n, cont):
if cont < n:
return 'S' + str(n - cont)
elif cont == n:
return 'CE'
else:
return 'B' + str(cont - n)
#Função de pertinência triangular
def triangular(a,b,c,xx):
#print('Triangular',a,b,c,xx)
if(xx < a):
return float(0)
elif(xx >= a and xx < b):
return float((xx-a)/(b-a))
elif(xx >= b and xx <= c):
return float((c-xx)/(c-b))
else:
return float(0)
#Gera uma chave única a partir do antecedente da regra
def antecedente(dados,colunas,regra):
key = ''
for c in range(colunas - 1):
col = dados.columns[c]
key += regra[col]
return key
#Erro quadrático médio
def SSE(xx,y):
erro = 0
for c in range(0,len(xx)):
tmp = xx[c]-y[c]
if isnan(tmp):
tmp = 1
erro = erro + (tmp)**2
return erro / len(xx)
#Imprime uma base de regras em formato amigável
def print_rule(regras,atributo_out):
for regra in regras:
out = ''
for atr in regra.keys():
if atr == atributo_out or atr == 'grau':
continue
if len(out) > 0:
out += ' E '
out += atr + ' é ' + regra[atr]
out = 'SE ' + out + ' ENTÃO ' + atributo_out + ' é ' + regra[atributo_out] + ' com grau ' + str(regra['grau'])
print(out)
#Plota o gráfico dos conjuntos 2D
def plot_conj2D(x,y,conjuntos,nomes):
minx = min(x)
maxx = max(x)
incx = (maxx - minx)/20
miny = min(y)
maxy = max(y)
incy = (maxy - miny)/20
plt.scatter(x,y)
labels = []
x = []
for subconj in conjuntos[nomes[0]].keys():
vals = conjuntos[nomes[0]][subconj]
x.append(vals[1])
labels.append(subconj)
plt.axvline(vals[1],label=subconj)
plt.plot([vals[0], vals[1], vals[2]],[miny,miny+incy,miny])
plt.xticks(x, labels, rotation='vertical')
labels = []
y = []
for subconj in conjuntos[nomes[1]].keys():
vals = conjuntos[nomes[1]][subconj]
y.append(vals[1])
labels.append(subconj)
plt.axhline(vals[1],label=subconj)
plt.plot([minx,minx+incx,minx],[vals[0], vals[1], vals[2]])
plt.yticks(y, labels, rotation='horizontal')
plt.xlabel(nomes[0])
plt.ylabel(nomes[1])
def plot_conj3D(x,y,z,conjuntos,nomes):
fig = plt.figure()
ax = Axes3D(fig)
minx = min(x)
maxx = max(x)
incx = (maxx - minx)/20
miny = min(y)
maxy = max(y)
incy = (maxy - miny)/20
minz = min(z)
maxz = max(z)
incz = (maxz - minz)/20
ax.scatter3D(x,y,z)
labels = []
x = []
for subconj in conjuntos[nomes[0]].keys():
vals = conjuntos[nomes[0]][subconj]
x.append(vals[1])
labels.append(subconj)
#ax.axvline(vals[1],label=subconj,zdir='x')
ax.plot([vals[0], vals[1], vals[2]],[miny,miny+incy,miny],zs=0,zdir='x')
plt.xticks(x, labels, rotation='vertical')
labels = []
y = []
for subconj in conjuntos[nomes[1]].keys():
vals = conjuntos[nomes[1]][subconj]
y.append(vals[1])
labels.append(subconj)
#ax.axhline(vals[1],label=subconj,zdir='y')
#ax.plot([minx,minx+incx,minx],[vals[0], vals[1], vals[2]],zdir='y')
plt.yticks(y, labels, rotation='horizontal')
labels = []
z = []
for subconj in conjuntos[nomes[2]].keys():
vals = conjuntos[nomes[2]][subconj]
z.append(vals[1])
labels.append(subconj)
#ax.axhline(vals[1],label=subconj,zdir='z')
#ax.plot([minz,minz+incz,minz],[vals[0], vals[1], vals[2]],zdir='z')
print(labels)
ax.set_zticks(z, labels)
plt.xlabel(nomes[0])
plt.ylabel(nomes[1])
ax.set_zlabel(nomes[2])
In [3]:
def fuzzy_rule(dados,particoes):
retorno = []
linhas,colunas = dados.shape
# ETAPA 1 - Particionamento do Espaço
conjuntos = dict()
for c in range(0,colunas):
col = dados.columns[c]
ma = max(dados[col])
mi = min(dados[col])
k = 2*particoes[c] + 1
tam = (ma - mi)/k
tam2 = tam/2
subconjunto = dict()
cont = 0
for part in np.arange(mi,ma,tam):
subconjunto[conj_nome(particoes[c],cont)] = [part - tam2, part + tam2, part + tam + tam2]
cont += 1
conjuntos[dados.columns[c]] = subconjunto
#return conjuntos
# ETAPA 2 - Regras
# ETAPA 2.1 - Pertinência dos conjuntos
pertinencias = []
for inst in range(0,linhas):
instperts = dict()
for c in range(colunas):
col = dados.columns[c]
conj = conjuntos[col]
perts = dict()
for subconj in conj.keys():
val = conj[subconj]
tmp = triangular(val[0],val[1],val[2],dados[col][inst])
perts[subconj] = tmp
instperts[col] = perts
pertinencias.append(instperts)
#return pertinencias
# ETAPA 2.1 - Regras a partir das pertinências
regras = []
for inst in range(linhas):
instregra = dict()
for c in range(colunas):
col = dados.columns[c]
conj = conjuntos[col]
maxp = 0
maxc = ''
for subconj in conj.keys():
if pertinencias[inst][col][subconj] > maxp: #Escolhe o subconj com maior pertinência
maxp = pertinencias[inst][col][subconj]
maxc = subconj
instregra[col] = maxc
regras.append(instregra)
#return regras
# ETAPA 3 - Graus das regras
# ETAPA 3.1 - Cálculo dos graus das regras
inst = 0
for regra in regras:
grau = 1.0
for c in range(colunas-1):
col = dados.columns[c]
grau = grau * pertinencias[inst][col][regra[col]]
inst += 1
regra['grau'] = grau
# ETAPA 3.2 - Escolha das regras de maior grau
nregras = dict()
for regra in regras:
ant = antecedente(dados,colunas,regra)
if ant in nregras:
if nregras[ant]['grau'] < regra['grau']:
nregras[ant] = regra
else:
nregras[ant] = regra
regras = []
for kreg in nregras.keys():
regra = nregras[kreg]
regras.append(regra)
retorno = []
retorno.append(conjuntos)
retorno.append(regras)
return retorno
def defuzzy_rule(instancia,atributos_in,atributo_out,conjuntos,regras):
#Etapa 5: Deffuzificação
ativacoes = []
centroides = []
for regra in regras:
ativ_regra = 1
for atr in atributos_in:
subconj = regra[atr]
conj = conjuntos[atr][subconj]
ativ = triangular(conj[0],conj[1],conj[2],instancia[atr])
ativ_regra *= ativ
centroides.append(conjuntos[atributo_out][regra[atributo_out]][1])
ativacoes.append(ativ_regra)
#print('Instancia:',instancia)
#print('Ativacoes:',ativacoes)
#print('Centroides:',centroides)
centroides_ponderadas = sum([float(centroides[k])*float(ativacoes[k]) for k in range(len(regras))])
ativacoes = sum([float(ativacoes[k]) for k in range(len(regras))])
return centroides_ponderadas / ativacoes
In [4]:
# Prepara o conjunto de treinamento
df = DataFrame([float(xx) for xx in np.arange(-1,1,0.1)],columns=['x'])
df['y'] = [xx*xx for xx in df['x']]
plot(df['x'],df['y'])
Out[4]:
In [5]:
#Extrai os conjuntos e as regras
conj,reg = fuzzy_rule(df,[20,20])
plot_conj2D(df['x'],df['y'],conj,['x','y'])
print_rule(reg,'y')
In [6]:
# Testa a aproximação
zz = []
for k in df['x']:
zz.append(defuzzy_rule({'x':k},['x'],'y',conj,reg))
plot(df['x'],zz)
Out[6]:
In [32]:
SSE(df['y'],zz)
Out[32]:
In [7]:
def sinc(xx):
return math.sin(xx)/xx
In [8]:
# Prepara o conjunto de treinamento
df2 = DataFrame(np.arange(0.1,2*math.pi,0.1),columns=['x'])
df2['y'] = [float(sinc(xx)) for xx in df2['x']]
plot(df2['x'],df2['y'])
Out[8]:
In [9]:
#Extrai os conjuntos e as regras
conj2,reg2 = fuzzy_rule(df2,[30,30])
print_rule(reg2,'y')
plot_conj2D(df2['x'],df2['y'],conj2,['x','y'])
In [10]:
# Testa a aproximação
zz = []
for k in df2['x']:
zz.append(defuzzy_rule({'x':k},['x'],'y',conj2,reg2))
plot(df2['x'],zz)
Out[10]:
In [48]:
#Erro quadrático médio
SSE(df2['y'],zz)
Out[48]:
In [11]:
# Prepara o conjunto de treinamento
df4 = DataFrame([[0,0,0]],columns=['x','y','z'])
X = []
Y = []
Z = []
for x in np.arange(-2*math.pi,2*math.pi,0.1):
tmpX = []
tmpY = []
tmpZ = []
for y in np.arange(-2*math.pi,2*math.pi,0.1):
z = sinc(x+y)
df4 = df4.append(DataFrame([[x,y,z]],columns=['x','y','z']))
tmpX.append(x)
tmpY.append(y)
tmpZ.append(z)
X.append(tmpX)
Y.append(tmpY)
Z.append(tmpZ)
df4.index = range(0,len(df4['x']))
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(X, Y, Z)
Out[11]:
In [12]:
#Gera os conjuntos e regras
conj4,reg4 = fuzzy_rule(df4,[3,3,3])
plot_conj3D(df4['x'],df4['y'],df4['z'],conj4,['x','y','z'])
In [43]:
#Testa a aproximação
ZZ = []
zz = [0]
for x in np.arange(-2*math.pi,2*math.pi,0.1):
tmpZ = []
for y in np.arange(-2*math.pi,2*math.pi,0.1):
z = defuzzy_rule({'x':x, 'y': y},['x','y'],'z',conj4,reg4)
zz.append(z)
tmpZ.append(z)
ZZ.append(tmpZ)
fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(X, Y, ZZ)
Out[43]:
In [47]:
#Erro quadrático Médio
SSE(df4['z'],zz)
Out[47]:
In [14]:
table1 =[[1.0, 0.0, -19],
[1.95,9.37,-17.95],
[2.88,18.23,-16.9],
[3.79, 26.59, -15.85],
[4.65, 34.44, -14.80],
[5.45, 41.78, -13.74],
[6.18,48.60, -12.7],
[7.48, 54.91,-11.65],
[7.99,60.71,-10.6],
[8.72,65.99,-9.55],
[9.01,70.75,-8.5],
[9.28,74.98,-7.45],
[9.46, 78.70,-6.4],
[9.59,81.90,-5.34],
[9.72,84.57,-4.3],
[9.81, 86.72,-3.25],
[9.88,88.34,-2.2],
[9.91,89.44,0.0]
]
In [15]:
table4 = [[7.0,0.0,-40.0],
[7.76, 18.75,-38.5],
[8.51,36.88,-37.0],
[9.15,54.4,-35.5],
[9.62,71.28,-34.0],
[9.88,87.51,-3.0],
[9.93,89.9,0]
]
In [16]:
table10= [[13.0,180.0,40.0],
[12.23,161.25,38.5],
[11.49,143.11,37.0],
[10.85,125.60,35.5],
[10.38,108.72,34.0],
[10.11,92.49,5.50],
[10.08,98.80,0.0]
]
In [17]:
#Cria o conjunto de treinamento
df5 = DataFrame(table1,columns=['x','phi','theta'])
df5 = df5.append(DataFrame(table4,columns=['x','phi','theta']))
df5 = df5.append(DataFrame(table10,columns=['x','phi','theta']))
df5.index = range(len(df5['x']))
In [18]:
#Cria os conjuntos e as regras
conj5,reg5 = fuzzy_rule(df5,[5,5,5])
print_rule(reg5,'theta')
plot_conj2D(df5['x'],df5['phi'],conj5,['x','phi'])
O método abaixo implementa a cinética do veículo, calcundo $x_{t+1}$ e $\phi_{t+_1}$ a partir do estado atual $(x_t, \phi_t)$ e $\theta$:
$$ x_{t+1} = x_t + cos(\phi_t+\theta_t) + sen(\theta_t)sen(\phi_t) $$$$ \phi_{t+1} = \phi_t - sen^{-1}(\frac{2sen(\theta_t)}{b}) $$
In [19]:
def calcula_cinetica(xx,phix,thetax,bx):
pr = math.radians(phix)
tr = math.radians(thetax)
xtmp = xx + math.cos(pr + tr) + math.sin(tr)*math.sin(pr)
phitmp = phix - math.degrees(math.asin(2*math.sin(tr)/bx))
ret = [xtmp, phitmp]
return ret
In [20]:
#Testa o modelo criado
xt = 1.0
phit = 0.0
xf = 10.0
phif = 90.0
x = []
u = []
y = []
v = []
yy = 20
ax = plt.axes()
while yy < 50 and abs(xf - xt) > 0.1:
x.append(xt)
y.append(yy)
thetat = defuzzy_rule({'x':xt, 'phi': phit},['x','phi'],'theta',conj5,reg5)
xtn, phit = calcula_cinetica(xt,phit,thetat,4)
ax.arrow(xt, yy, xtn-xt-0.1, -0.8,length_includes_head=True,head_width=0.2)
yy -= 1
xt = xtn
scatter(x,y)
Out[20]:
A série temporal caótica de Mackey-Glass é gerada pela equação: $$ \frac{dx(t)}{dt} = \beta\frac{ x(t - \tau)}{1 + x(t - \tau)^k} - \gamma x(t) $$
In [21]:
def mackey_glass_ts(n,tau, beta, gamma,k):
xret = [rnd.random() for k in range(tau)]
for t in range(tau-1,n):
xret.append( beta*((xret[t - tau]) / (1 + math.pow(xret[t - tau],k))) - gamma*xret[t] )
return xret
In [22]:
#Cria o conjunto de dados Mackey-Glass para tau = 50, beta = 0.9, gamma = 0.5 e k = 10
ts = mackey_glass_ts(500,50,0.9,0.5,10)
fig = plt.figure(figsize=[20,6])
plot(ts)
Out[22]:
In [23]:
# Cria o conjunto de dados de treinamento
def cria_dataframe(timeserie,m):
labels = []
dados = []
for m1 in range(m,-1,-1):
labels.append('t'+str(m1))
for count in range(m + 1,len(timeserie)):
tmp = [0 for xx in range(0,m+1)]
for m1 in range(0,m+1):
tmp[m1] = timeserie[count - m-m1]
#print(tmp)
dados.append(tmp)
dfts = DataFrame(dados,columns=labels)
dfts.index = range(len(dfts['t0']))
return dfts
In [24]:
df_ts = cria_dataframe(ts,9)
df_ts[0:50]
Out[24]:
In [25]:
#Cria os conjuntos e as regras
particoes = [10 for xx in range(10)]
conj6,reg6 = fuzzy_rule(df_ts,particoes)
print_rule(reg6,'t0')
In [26]:
# Visualiza os conjuntos
plot_conj2D(df_ts['t7'],df_ts['t0'],conj6,['t7','t0'])
In [27]:
ii = {'t0':0}
for ic in range(9,0,-1):
ii['t'+str(ic)] = ts[200-ic]
tmp = defuzzy_rule(ii,['t'+str(xx) for xx in range(9,0,-1)],'t0',conj6,reg6)
In [28]:
#Testa a predição da série
zz = []
for t in range(200,400):
ii = {'t0':0}
for ic in range(9,0,-1):
ii['t'+str(ic)] = ts[t-ic]
tmp = defuzzy_rule(ii,['t'+str(xx) for xx in range(9,0,-1)],'t0',conj6,reg6)
ii['t0'] = tmp
zz.append(tmp)
#zz
In [29]:
#Visualiza a predição
fig = plt.figure(figsize=[20,6])
plt.ylabel('z(t)')
plt.xlabel('t')
po, = plot(ts[200:400],color='b',label='Original')
pp, = plot(zz,color='r',label='Prevista')
plt.legend(handles=[po,pp])
fig.show()
In [49]:
SSE(ts[200:400],zz)
Out[49]:
In [30]: