Fonctions Floues
unknown
plain_text
2 years ago
42 kB
9
Indexable
from fonctions_autres import *
import numpy as np
import matplotlib.pyplot as plt
import scipy
import skfuzzy as fuzz
from skfuzzy import control as ctrl
from tqdm import tqdm
def input_fuzz(variable, input_to_fuzz):
fuzz_input = []
for vari in variable.terms.keys():
fuzz_input.append(fuzz.interp_membership(variable.universe, variable.terms[vari].mf, input_to_fuzz))
return fuzz_input
def IRR_2var(SF_rules, input_fuzz1, input_fuzz2):
IRR = np.empty([len(SF_rules), 2])
k = -1
for j in range(0, len(SF_rules)):
if (j % len(input_fuzz2)) == 0:
k += 1
IRR[j, 0] = input_fuzz1[k]
IRR[j, 1] = input_fuzz2[j % len(input_fuzz2)]
return IRR
def IRR_3var(SF_rules, input_fuzz1, input_fuzz2, input_fuzz3):
IRR = np.empty([len(SF_rules), 3])
k = -1
m = -1
for j in range(0, len(SF_rules)):
if (j % len(input_fuzz3)) == 0:
k += 1
if (k % len(input_fuzz2)) == 0:
m += 1
IRR[j, 0] = input_fuzz1[m]
IRR[j, 1] = input_fuzz2[j % len(input_fuzz2)]
IRR[j, 2] = input_fuzz3[j % len(input_fuzz3)]
return IRR
def output_fuzz(csq, output_var):
first_rule = True
output = np.zeros(1)
for membi in output_var.terms:
csq_i = '[' + str(output_var.terms[membi]) + ']'
if csq[csq_i] != 0:
x_csq = fuzz.fuzzymath.interp_universe(output_var.universe, output_var.terms[membi].mf, 1)
if first_rule:
output = np.array([[np.max(x_csq), csq[csq_i]]])
first_rule = False
else:
output = np.append(output, [[(np.min(x_csq)), csq[csq_i]]], axis=0)
return output
def SF1_compute(input_temperature, input_d_speed, input_time, want_to_plot):
pas = 0.01
temperature = ctrl.Antecedent(np.arange(263, 323, pas), 'Temperature (K)')
d_speed = ctrl.Antecedent(np.arange(-100, 100, pas), 'Variation of wind speed (m/s)')
timeTC = ctrl.Antecedent(np.arange(0, 15, pas), 'Time since TC formation (days)')
lifespan = ctrl.Consequent(np.arange(0, 15, pas), 'Lifespan')
temperature['Low'] = fuzz.trapmf(temperature.universe, [263, 263, 293, 303])
temperature['Ok'] = fuzz.trapmf(temperature.universe, [293, 303, 323, 323])
d_speed['--'] = fuzz.trapmf(d_speed.universe, [-100, -100, -60, -20])
d_speed['-'] = fuzz.trimf(d_speed.universe, [-60, -30, 0])
d_speed['St'] = fuzz.trimf(d_speed.universe, [-30, 0, 20])
d_speed['+'] = fuzz.trimf(d_speed.universe, [0, 30, 60])
d_speed['++'] = fuzz.trapmf(d_speed.universe, [20, 60, 100, 100])
timeTC['Young'] = fuzz.zmf(timeTC.universe, 3, 6)
timeTC['Middle'] = fuzz.gaussmf(timeTC.universe, 6, 1)
timeTC['Old'] = fuzz.smf(timeTC.universe, 6, 9)
lifespan['Dying+'] = fuzz.trapmf(lifespan.universe, [0, 0, 3, 5])
lifespan['Dying'] = fuzz.trimf(lifespan.universe, [3, 5, 7])
lifespan['St'] = fuzz.trimf(lifespan.universe, [5, 7, 10])
lifespan['Living'] = fuzz.trimf(lifespan.universe, [7, 10, 13])
lifespan['Living+'] = fuzz.trapmf(lifespan.universe, [10, 13, 15, 15])
if want_to_plot:
lifespan.view()
temperature.view()
d_speed.view()
timeTC.view()
SF1_rules = [ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Young'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Middle'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['--'] & timeTC['Old'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Middle'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['-'] & timeTC['Old'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Middle'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Young'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Middle'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Old'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Middle'], lifespan['Living']),
ctrl.Rule(temperature['Low'] & d_speed['++'] & timeTC['Old'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Middle'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['--'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Young'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Middle'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['-'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Middle'], lifespan['St']),
ctrl.Rule(temperature['Ok'] & d_speed['St'] & timeTC['Old'], lifespan['Dying']),
ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Middle'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['+'] & timeTC['Old'], lifespan['Living']),
ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Young'], lifespan['Living+']),
ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Middle'], lifespan['Living+']),
ctrl.Rule(temperature['Ok'] & d_speed['++'] & timeTC['Old'], lifespan['Living+']),
]
SF1_ctrl = ctrl.ControlSystem(SF1_rules)
SF1 = ctrl.ControlSystemSimulation(SF1_ctrl)
input_temperature_fuzz = input_fuzz(temperature, input_temperature)
input_d_speed_fuzz = input_fuzz(d_speed, input_d_speed)
input_time_fuzz = input_fuzz(timeTC, input_time)
IRR = IRR_3var(SF1_rules, input_temperature_fuzz, input_d_speed_fuzz, input_time_fuzz)
declenchement = IRR.min(axis=1)
csq_SF1 = {} # variable de type dictionnaire
for rule_i in range(0, len(SF1_rules)): # on initialise à 0 pour chaque règle
csq_SF1[str(SF1_rules[rule_i].consequent)] = 0
for rule_i in range(0, len(SF1_rules)):
csq_i = str(SF1_rules[rule_i].consequent)
csq_SF1[csq_i] = max(csq_SF1[csq_i], declenchement[rule_i]) # Comme dans le TP6
return output_fuzz(csq_SF1, lifespan)
def SF2_compute(input_humidity, input_TC_size, want_to_plot):
pas = 0.01
humidity = ctrl.Antecedent(np.arange(0, 100, pas), 'humidity') # Humidity (%)'
TC_size = ctrl.Antecedent(np.arange(0, 350, pas), 'TC_size') # TC size (km)
TC_size_variation = ctrl.Consequent(np.arange(0, 60, pas), 'TC_var') # Variation of TC size for next day (km)
humidity['--'] = fuzz.trapmf(humidity.universe, [0, 0, 20, 40]) # 20%
humidity['-'] = fuzz.trimf(humidity.universe, [20, 40, 60]) # 40%
humidity['+'] = fuzz.trimf(humidity.universe, [40, 60, 80]) # 60%
humidity['++'] = fuzz.trapmf(humidity.universe, [60, 80, 100, 100]) # 80%
TC_size['Little+'] = fuzz.trapmf(TC_size.universe, [0, 0, 75, 120])
TC_size['Little'] = fuzz.trimf(TC_size.universe, [75, 150, 200])
# TC_size['Normal'] = fuzz.trimf(TC_size.universe, [120, 150, 200])
TC_size['Huge'] = fuzz.trimf(TC_size.universe, [150, 200, 250])
TC_size['Huge+'] = fuzz.trapmf(TC_size.universe, [200, 250, 350, 350])
# 10 13 15 20 25 50
TC_size_variation['Very Low'] = fuzz.trapmf(TC_size_variation.universe, [0, 0, 10, 13]) # 10
TC_size_variation['Low'] = fuzz.trimf(TC_size_variation.universe, [10, 13, 15]) # 13
TC_size_variation['Quite Low'] = fuzz.trimf(TC_size_variation.universe, [13, 15, 20]) # 15
TC_size_variation['Quite High'] = fuzz.trimf(TC_size_variation.universe, [15, 20, 25]) # 20
TC_size_variation['High'] = fuzz.trimf(TC_size_variation.universe, [20, 25, 50]) # 25
TC_size_variation['Very High'] = fuzz.trapmf(TC_size_variation.universe, [25, 50, 50, 50]) # 50
if want_to_plot:
humidity.view()
TC_size.view()
TC_size_variation.view()
SF2_rules = [ctrl.Rule(humidity['--'] & TC_size['Little+'], TC_size_variation['Very Low']),
ctrl.Rule(humidity['--'] & TC_size['Little'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['--'] & TC_size['Huge'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['--'] & TC_size['Huge+'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['-'] & TC_size['Little+'], TC_size_variation['Low']),
ctrl.Rule(humidity['-'] & TC_size['Little'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['-'] & TC_size['Huge'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['-'] & TC_size['Huge+'], TC_size_variation['Quite Low']),
ctrl.Rule(humidity['+'] & TC_size['Little+'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['+'] & TC_size['Little'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['+'] & TC_size['Huge'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['+'] & TC_size['Huge+'], TC_size_variation['Quite High']),
ctrl.Rule(humidity['++'] & TC_size['Little+'], TC_size_variation['High']),
ctrl.Rule(humidity['++'] & TC_size['Little'], TC_size_variation['Very High']),
ctrl.Rule(humidity['++'] & TC_size['Huge'], TC_size_variation['Very High']),
ctrl.Rule(humidity['++'] & TC_size['Huge+'], TC_size_variation['Very Low'])]
SF2_ctrl = ctrl.ControlSystem(SF2_rules)
SF2 = ctrl.ControlSystemSimulation(SF2_ctrl)
input_TC_size_fuzz = input_fuzz(TC_size, input_TC_size)
input_humidity_fuzz = input_fuzz(humidity, input_humidity)
IRR = IRR_2var(SF2_rules, input_humidity_fuzz, input_TC_size_fuzz)
declenchement = IRR.min(axis=1)
csq_SF2 = {} # variable de type dictionnaire
for rule_i in range(0, len(SF2_rules)): # on initialise à 0 pour chaque règle
csq_SF2[str(SF2_rules[rule_i].consequent)] = 0
for rule_i in range(0, len(SF2_rules)):
csq_i = str(SF2_rules[rule_i].consequent)
csq_SF2[csq_i] = max(csq_SF2[csq_i], declenchement[rule_i]) # Comme dans le TP6
return output_fuzz(csq_SF2, TC_size_variation) # Dictionnaire des conséquences, et variable de sortie
def fuzz_CAF1(uwind, vwind, x_TC, y_TC):
# Fuzzification des matrices, avec des NFT
# Soit x le nombre net, alors cela revient à l'ajout d'une dimension à la matrice, de valeur : [x+a,x+b,x+c,x+d]
a, b, c, d = [-2, -1, 1, 2] # Augmenter l'envergure augmente la zone des cases dans le graph final
uwind = uwind[:, :, np.newaxis] # Extension de la matrice en 3D
# Fuzzification uwind vwind
# Créer une dimension pour les coordonnées des nombres trapézoidaux, de forme [a,b,c,d]
uwind = np.dstack((uwind + a, uwind + b, uwind + c, uwind + d)) # La matrice est maintenant fuzzifiée
a, b, c, d = [-2, -1, 1,2]
vwind = vwind[:, :, np.newaxis]
vwind = np.dstack((vwind + a, vwind + b, vwind + c, vwind + d)) # La matrice est maintenant fuzzifiée
dNFT = [-10000, -5000, 5000, 10000] #Faut pas trop augmenter sinon ça deviens incohérent
x_TC = np.array(x_TC + dNFT)
y_TC = np.array(y_TC + dNFT)
return uwind, vwind, x_TC, y_TC
def CAF1_compute(uwind, vwind, x_TC, y_TC):
# Calcul PSI, PSI_TC
# ---------
# PSI : Ce dernier est calculé à partir d'une intégrale, donc d'une somme et d'une division par un nombre net
# U et V sont des IFT de la forme (a,b,c,d)
# Nous pouvons calculer l'intégrale pour a,b,c et d, 1 à 1
# Intégrales de U et V pour déterminer psi, à l'instant t
# Cf.Formule
intx = scipy.integrate.cumtrapz(vwind, lon, axis=1, initial=0)
inty = np.flip(scipy.integrate.cumtrapz(uwind, lat, axis=0, initial=0),axis=2)
#############
# POURQUOI FAUT FLIP INTY??????
#############
# La soustraction nécessite un np.flip
psi = (inty - np.flip(intx,axis=2)) / 2
# ----------
# PSI_TC : Comme au dessus, PSI_TC se base sur une moyenne, donc par le calcul d'additions et d'une division par un nombre net
# Nous pouvons calculer directement PSI_TC pour les 4 termes qui composent notre IFT
# Nous allons le déterminer avec les coordonnées polaires
# L'idée est de longer selon l'axe x, et de calculer PSI_TC sur le cercle de rayon x
res = 3000
psi_TC_xy = np.zeros(np.shape(psi))
psi_TC_r = np.zeros((res, 4))
K = np.zeros(4)
for i in range(0, 4):
univers = np.linspace(x_TC[i], max(lon) + 5000000, res)
# On va au delà de lon, nécessaire car on travaille avec des cercles, faut sortir du carré
# L'Univers est pour définir l'espace de travail de x
# On veut PSI_TC selon R (une liste), et PSI_TC selon (x,y) (un champ 2D)
psi_TC_r[:, i], psi_TC_xy[:, :, i] = psi_TC_compute(x_TC[i], y_TC[i], univers, psi[:, :, i], is_fuzzy=True)
# Soustraction par psi_inf.
psi_inf = np.mean(np.flip(psi, axis=2)[:, :, i])
psi_TC_r[:, i] = [x - psi_inf for x in psi_TC_r[:, i]]
psi_TC_xy[:, :, i] = psi_TC_xy[:, :, i] - psi_inf
K[i] = K_compute(psi_TC_r[:, i], univers)
# Le paramètre K nous sera aussi nécessaire
psi_large = (psi - np.flip(psi_TC_xy, axis=2))
return psi, psi_large, psi_TC_xy, K
def xy_pred(x0, y0, u_large, v_large, n_grad_x, n_grad_y, r_TC, t_1s, dts, K,SF3_ouput):
x_TC_out = x0 # On prépare la sortie
y_TC_out = y0
# A partir des champs et de la position du TC, nous pouvons prédire sa trajectoire pour dans un temps t_1s
# La vitesse varie dans le champ, nous ferons donc des "petits pas" de temps dts, pour arriver à l'instant final
# Retourne la position finale (x;y)
# Cependant il n'est pas intéressant de prendre les valeurs de vitesse en 1 point : on s'intéresse à une zone
# de r mètres autour du cyclone. On prendra la vitesse moyenne dans cette zone
# _____________
# Divers paramètres pour le calcul :
r_min = 20000
C = 0.7156
nm = 1.2830
a = 1.5
for k in range(0, int(np.round(t_1s)), int(np.round(dts))):
r_map, angle_map = cart_to_polar(x_TC_out, y_TC_out, lon, lat)
u_large_r1 = []
v_large_r1 = []
n_grad_y_r1 = []
n_grad_x_r1 = []
r_univers = np.arange(r_min, r_TC, 10000)
for r in r_univers:
u_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, u_large)))
v_large_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, v_large)))
n_grad_y_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_y)))
n_grad_x_r1.append(np.nanmean(array_r1_compute(r, r_map, angle_map, n_grad_x)))
# FLOU
# Ces lignes là seront à changer pour le flou
u = np.mean(u_large_r1)
v = np.mean(v_large_r1)
ny = np.mean(n_grad_y)
nx = np.mean(n_grad_x)
#
TC_lat = np.radians(met_to_deg(y_TC_out)) # Détermination de la latitude du cyclone, en rad
# Calcul du paramètre m. Celui-ci corrige 'l'erreur' lié au fait que nous travaillons sur une sphère, et donc
# que les distances sont courbées comparées à une étude sur un plan 2D
m = (C * (np.cos(TC_lat)) ** (-1)) * (np.tan(0.25 * (np.pi - 2 * TC_lat))) ** nm
# idx_x, idx_y = xy_to_indx_lonlat(x_TC_out, y_TC_out, lon, lat) # On récupère les index selon [lon,lat] du TC
Cx0 = a * u # - K * ny # Vélocité selon x au centre du cyclone, en m/s
Cy0 = a * v # + K * nx # Vélocité selon y au centre du cyclone, en m/s
# On actualise la position
x_TC_out = x_TC_out + Cx0 * dts * m # m > 0
y_TC_out = y_TC_out + Cy0 * dts * m
return x_TC_out, y_TC_out
def CAF2_compute(uwind, vwind, psi, psi_TC, psi_large, x_TC, y_TC, r_TC, temps_etude, K,SF3_output):
f = np.expand_dims(lat, axis=1) * np.ones(np.shape(psi[:, :, 0]))
f = parametre_coriolis(f)
f = np.dstack(([f, f, f, f]))
delta = lat[1] - lat[0] # Résolution de nos données
n = scipy.ndimage.laplace(psi_large) / (delta ** 2) + 1 * f
# u_TC = fuzzy_gradient(psi_TC, delta=lat[1] - lat[0], axis=0)
# v_TC = -np.flip(fuzzy_gradient(psi_TC, delta=lon[1] - lon[0], axis=1), axis=2)
u_TC = np.gradient(psi_TC, lat, axis=0)
v_TC = -np.flip(np.gradient(psi_TC, lon, axis=1),axis=2)
v_large = vwind[:, :, :] - np.flip(v_TC, axis=2) # La négation du nombre flou demande le np.flip()
u_large = uwind[:, :, :] - np.flip(u_TC, axis=2)
n_grad_y = np.array(np.gradient(n, lat, axis=0))
n_grad_x = np.array(np.gradient(n, lon, axis=1))
# Il vont être utilisés via une soustraction
n_grad_y = np.flip(n_grad_y, axis=2)
# On lisse les données
# smooth = 5
# n_grad_y = smooth_array(n_grad_y, lon, smooth, axis=1)
# n_grad_x = smooth_array(n_grad_x, lat, smooth, axis=0)
# On a tous les paramètres nécessaires au calcul de la vélocité, et donc de la position future du cyclone
# La vélocité dépend des champs n et PSI. Ces derniers n'étant pas uniformes, la vélocité n'est pas la même partout
# Il est donc nécessaire d'étudier la vélocité à petit pas, avec un pas de temps dt
# Soit t_1 (heures) l'instant futur du cyclone
# Nous travaillons en unités SI
t_1s = temps_etude * 60 * 60
dts = (temps_etude / 3) * 60 * 60
# Idée pour le flou : Ici r_TC c'est la taille max ; Mais ça sera un nb flou.
# On trace la moyenne de U,V,N en fonction du rayon
# Pour chaque valeur du rayon on multiplie par le degré de vérité de r_TC
# Pour chaque pas de temps dt, allant de l'instant t0 (maintenant) à l'instant t_1
# Désormais il sera obligatoire de différencier les x en fonction de leur hauteur
h = 1
x_TC = np.append(x_TC, h)
y_TC = np.append(y_TC, h)
# Constitué de plein de IFT. 1ère dim : chaque point.
x_field = np.zeros((1, 5))
y_field = np.zeros((1, 5))
x_field[0] = x_TC
y_field[0] = y_TC
# Au départ x_field n'est composé que d'un point
for t in range(0, t_1s, int(dts)):
# Pour chaque point i flou ; donc pour chaque ligne
x_field_1 = np.zeros((1, 5))
y_field_1 = np.zeros((1, 5))
print("Calcul...", t / dts)
for i in tqdm(range(0, np.size(x_field[:, 0]))):
# On récupère les coordonnées
xi_fuzz = np.transpose(np.expand_dims(x_field[i, :-1], axis=0)) # Vecteur colonne
yi_fuzz = np.transpose(np.expand_dims(y_field[i, :-1], axis=0))
# On met le nombre (singleton) sous l'écriture de Kaufman.
yi_fuzz = yi_fuzz * np.ones((4, 5))
xi_fuzz = xi_fuzz * np.ones((4, 5))
# Pour la hauteur :
xi_fuzz[:, 4] = x_field[i, 4]
yi_fuzz[:, 4] = y_field[i, 4]
# A l'extrémité des IFT, on a forcément h=0
xi_fuzz[3::4, 4] = 0
yi_fuzz[3::4, 4] = 0
xi_fuzz[::4, 4] = 0
yi_fuzz[::4, 4] = 0
# Ajout des alpha_coupes
# Alpha = 0.5 ; La hauteur est de moitié
to_insert = (xi_fuzz[0, :] + xi_fuzz[1, :]) / 2
to_insert1 = (xi_fuzz[2, :] + xi_fuzz[3, :]) / 2
xi_fuzz = np.insert(xi_fuzz, [1], to_insert, axis=0)
xi_fuzz = np.insert(xi_fuzz, [-1], to_insert1, axis=0)
to_insert = (yi_fuzz[0, :] + yi_fuzz[1, :]) / 2
to_insert1 = (yi_fuzz[2, :] + yi_fuzz[3, :]) / 2
yi_fuzz = np.insert(yi_fuzz, [1], to_insert, axis=0)
yi_fuzz = np.insert(yi_fuzz, [-1], to_insert1, axis=0)
# Nous avons xy_fuzz (net), en 1 point, nous calculons son déplacement (flou)
# Comme il s'agit plus précisément d'une addition, alors on procède élément par élément
for m in range(0, np.size(xi_fuzz[:, 0])):
for j in range(0, 4): # Il ne va pas à 4
xi_fuzz[m, j], yi_fuzz[m, j] = xy_pred(xi_fuzz[m, j], yi_fuzz[m, j], u_large[:, :, j],
v_large[:, :, j],
n_grad_x[:, :, j], n_grad_y[:, :, j],
r_TC, dts, dts/2, K[j],SF3_output)
# Nous avons le déplacement flou
x_field_1 = np.append(x_field_1, xi_fuzz, axis=0)
y_field_1 = np.append(y_field_1, yi_fuzz, axis=0)
# On enlève ceux avec une possibilité trop faible
x_field = x_field_1[x_field_1[:, 4] >= 0.1]
y_field = y_field_1[y_field_1[:, 4] >= 0.1]
print('u_large.mean = ', np.mean(u_large, axis=(0, 1)))
print('v_large.mean = ', np.mean(v_large, axis=(0, 1)))
print('n_grad_x.mean = ', np.mean(n_grad_x, axis=(0, 1)))
print('n_grad_y.mean = ', np.mean(np.flip(n_grad_y, axis=2), axis=(0, 1)))
return x_field, y_field
def IFT_appartenance(a, univers):
a1, a2, a3, a4, h = a
mu = float("nan") * univers
mu[(univers < a1)] = 0
mu[(univers > a1) & (univers < a2)] = h * (univers[(univers > a1) & (univers < a2)] - a1) / (a2 - a1)
mu[(univers > a2) & (univers < a3)] = h
mu[(univers > a3) & (univers < a4)] = h * (a4 - univers[(univers > a3) & (univers < a4)]) / (a4 - a3)
mu[(univers > a4)] = 0
return mu
def xy_fuzz_plot(x, y, lon_plot, lat_plot):
# Soit un point de coordonnées X et Y IFT notation de Kaufmann
# Ressort un champ /meshgrid avec les degrés de possibilités
# Unités SI (mètres)
grid = np.meshgrid(lon_plot, lat_plot)
plot_grid = 0 * grid[0]
print("___________AFFICHAGE_EN_COURS_________")
for i in tqdm(range(0, np.size(x[:, 0]))): # Pour chaque point flou
mu_x = IFT_appartenance(x[i, :], grid[0])
mu_y = IFT_appartenance(y[i, :], grid[1])
xy_grid = mu_x * mu_y # FAUT DEFINIR UNE T NORME
# là où il n'y a pas de 0 dans xy_grid on rajoute xy-grid, sinon on met plot_grid (= on change rien)
mask = np.invert(xy_grid == 0)
plot_grid = np.where(mask, np.maximum(plot_grid, xy_grid), plot_grid)
plot_array(plot_grid, lon_plot, lat_plot)
def SF3_compute(input_TC_velocity, input_surface_pressure, input_dsurface_pressure, want_to_plot):
velocity = ctrl.Antecedent(np.arange(0, 100, 0.01), 'Velocity') # Velocity (m/s)
surface_pressure = ctrl.Antecedent(np.arange(92000, 103000, 10), 'Surface_Pressure') # Surface pressure (Pa)
dsurface_pressure = ctrl.Antecedent(np.arange(-2000, 2000, 10),
'Surface_Pressure_Variation') # Surface pressure variation (Pa/h) Surface variation per hour or per day?
velocity_variation = ctrl.Consequent(np.arange(-25, 25, 0.1),
'Velocity_Variation') # Variation of the velocity next hour (m/s). Next hour / next day ?
velocity['Little'] = fuzz.trapmf(velocity.universe, [0, 0, 20, 40])
velocity['Intermediate'] = fuzz.trapmf(velocity.universe, [20, 30, 50, 60])
velocity['High'] = fuzz.trapmf(velocity.universe, [40, 60, 100, 100])
surface_pressure['Little'] = fuzz.trapmf(surface_pressure.universe, [92000, 92000, 99000, 101000])
surface_pressure['High'] = fuzz.trapmf(surface_pressure.universe, [99000, 101000, 103000, 103000])
dsurface_pressure['--'] = fuzz.trapmf(dsurface_pressure.universe, [-2000, -2000, -1500, -500])
dsurface_pressure['-'] = fuzz.trimf(dsurface_pressure.universe, [-1500, -500, 0])
dsurface_pressure['0'] = fuzz.trimf(dsurface_pressure.universe, [-500, 0, 500])
dsurface_pressure['+'] = fuzz.trimf(dsurface_pressure.universe, [0, 500, 1500])
dsurface_pressure['++'] = fuzz.trapmf(dsurface_pressure.universe, [500, 1500, 2000, 2000])
velocity_variation['--'] = fuzz.trimf(velocity_variation.universe, [-25, -25, -10])
velocity_variation['-'] = fuzz.trapmf(velocity_variation.universe, [-25, -10, -5, 0])
velocity_variation['St'] = fuzz.trimf(velocity_variation.universe, [-5, 0, 5])
velocity_variation['+'] = fuzz.trapmf(velocity_variation.universe, [0, 5, 10, 25])
velocity_variation['++'] = fuzz.trimf(velocity_variation.universe, [10, 25, 25])
if want_to_plot:
velocity.view()
surface_pressure.view()
dsurface_pressure.view()
velocity_variation.view()
SF3_rules = [
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['--'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['-'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['+'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['++'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['--'], velocity_variation['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['-'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['0'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['+'], velocity_variation['-']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['++'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['--'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['-'], velocity_variation['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['++'], velocity_variation['--']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['--'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['-'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['+'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['++'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['--'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['-'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['0'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['+'], velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['++'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['--'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['-'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['++'], velocity_variation['-'])]
SF3_ctrl = ctrl.ControlSystem(SF3_rules)
SF3 = ctrl.ControlSystemSimulation(SF3_ctrl)
input_velocity_fuzz = input_fuzz(velocity, input_TC_velocity)
input_surface_pressure_fuzz = input_fuzz(surface_pressure, input_surface_pressure)
input_dsurface_pressure_fuzz = input_fuzz(dsurface_pressure, input_dsurface_pressure)
IRR = IRR_3var(SF3_rules, input_surface_pressure_fuzz, input_velocity_fuzz, input_dsurface_pressure_fuzz)
declenchement = IRR.min(axis=1)
csq_SF3 = {} # variable de type dictionnaire
for i in range(0, len(SF3_rules)): # on initialise à 0 pour chaque règle
csq_SF3[str(SF3_rules[i].consequent)] = 0
for i in range(0, len(SF3_rules)):
csq_i = str(SF3_rules[i].consequent)
csq_SF3[csq_i] = max(csq_SF3[csq_i], declenchement[i]) # Comme dans le TP6
return output_fuzz(csq_SF3, velocity_variation)
def estimation2(t0, t1, dt, uwind, vwind, study_grid):
# Estimation flou de la trajectoire (sortie CAF2), de t0 à t1
print("_________________________SF2_______________________")
SF2_humidity = np.where(study_grid, humidity[t0, :, :], float('nan'))
SF2_humidity = np.nanmean(SF2_humidity)
SF2_taille_TC0 = taille_TC(t0, study_grid) # Rayon, en km
dtaille = SF2_compute(SF2_humidity, SF2_taille_TC0,
want_to_plot=False) # Prediction pour dans 24 heures, en pourcentage
dtaille[:, 0] = (t1 - t0) * (dtaille[:, 0] / 24) # On le ramène à une prévision au temps t1
dtaille[:, 0] = (dtaille[:, 0] * 0.01)
dtaille[:, 0] = (dt) * dtaille[:, 0] / 24 # Pour les calculs de vitesse nous voudrons l'évolution pour chaque dt
SF2_taille_TC1 = dtaille
SF2_taille_TC1[:, 0] = SF2_taille_TC0 * (1+dtaille[:, 0])
print("Taille Actuelle du cyclone : ", np.rint(SF2_taille_TC0), "Humidité :", np.rint(SF2_humidity))
print("_________________________SORTIE SF2 - Future Taille du Cyclone à dt___")
print("Prédiction :", SF2_taille_TC1, 'Réel :', np.rint(taille_TC(t0 + dt, study_grid)))
print("_____________________________CAF1____________")
uwind_fuzz, vwind_fuzz, x_TC_fuzz, y_TC_fuzz = fuzz_CAF1(uwind[t0, :, :], vwind[t0, :, :], x_TC[t0], y_TC[t0])
print('X =', x_TC_fuzz, 'Y = ', y_TC_fuzz)
print('uwind.mean = ', np.rint(np.mean(uwind_fuzz, axis=(0, 1))))
print('vwind.mean = ',np.rint(np.mean(vwind_fuzz, axis=(0, 1))))
print('____')
psi, psi_large, psi_TC, K = CAF1_compute(uwind_fuzz, vwind_fuzz, x_TC_fuzz, y_TC_fuzz)
print('Psi.mean =', np.rint(np.mean(psi, axis=(0, 1))))
print('Psi_large.mean = ', np.rint(np.mean(psi_large, axis=(0, 1))))
print('Psi_TC.mean = ', np.rint(np.mean(psi_TC, axis=(0, 1))))
print('K', np.rint(K))
print('_______________CAF1 - DONE___________________')
print("_______________SF3___________________")
#SF3_speed_TC = vitesse_TC(0,t0,calcul_pth())
#SF3_pressure_TC = pressure[t0,idx_y_TC[t0],idx_x_TC[t0]]
#SF3_dpressure_TC = pressure[t0,idx_y_TC[t0],idx_x_TC[t0]]
SF3_output = 0#SF3_compute(SF3_speedTC,SF3_pressure_TC,SF3_dpressure_TC)
print("_____________________SF3 - DONE______________")
x, y = CAF2_compute(uwind_fuzz, vwind_fuzz, psi, psi_TC, psi_large, x_TC_fuzz, y_TC_fuzz, 200000,
t1 - t0,
K,SF3_output)
print('Il y a ', np.size(x[:, 0]), "points !")
print('________________CAF2 - DONE________________________')
# plot_array(psi[:, :, 0], lon, lat)
# for i in range(0, np.size(x[:, 0])):
# if x[i, 4] == 1:
# plt.plot(met_to_deg(x[i, 0:4]), met_to_deg(y[i, 0:4]), ms=1, marker='x', markeredgecolor='red')
plt.plot(met_to_deg(x_TC[t0]), met_to_deg(y_TC[t0]), ms=10, marker='o', markeredgecolor='red')
plt.plot(met_to_deg(x_TC[t1]), met_to_deg(y_TC[t1]), ms=10, marker='x', markeredgecolor='red')
# plt.show()
# On augmente la résolution pour l'affichage des degrés de possibilités, sinon on n'est pas assez précis
lon_fuzz_plot = np.arange(np.min(lon), np.max(lon), ((lon[1]) - (lon[0])) / 6)
lat_fuzz_plot = np.arange(np.max(lat), np.min(lat), ((lat[1]) - (lat[0])) / 6)
xy_fuzz_plot(x, y, lon_fuzz_plot, lat_fuzz_plot)
plt.show()
def SF4_compute(input_actu, input_ant, want_to_plot):
tendance_actuelle = ctrl.Antecedent(np.arange(-90, 90, 0.1), 'Angle calcul')
tendance_antecedents = ctrl.Antecedent(np.arange(-90, 90, 0.1), 'Angle évalué par la forme')
angle = ctrl.Consequent(np.arange(-90, 90, 0.1), 'Prédiction courbature')
tendance_actuelle['Courbé Gauche Fort'] = fuzz.trapmf(tendance_actuelle.universe, [-90, -90, -60, -45])
tendance_actuelle['Courbé Gauche'] = fuzz.trapmf(tendance_actuelle.universe, [-60, -45, -10, 0])
tendance_actuelle['Linéaire'] = fuzz.trapmf(tendance_actuelle.universe, [-10, 0, 0, 10])
tendance_actuelle['Courbé Droite'] = fuzz.trapmf(tendance_actuelle.universe, [0, 10, 45, 60])
tendance_actuelle['Courbé Droite Fort'] = fuzz.trapmf(tendance_actuelle.universe, [45, 60, 90, 90])
tendance_antecedents['Courbé Gauche Fort'] = fuzz.trapmf(tendance_antecedents.universe, [-90, -90, -60, -45])
tendance_antecedents['Courbé Gauche'] = fuzz.trapmf(tendance_antecedents.universe, [-60, -45, -20, 0])
tendance_antecedents['Linéaire'] = fuzz.trapmf(tendance_antecedents.universe, [-20, 0, 5, 20])
tendance_antecedents['Courbé Droite'] = fuzz.trapmf(tendance_antecedents.universe, [5, 20, 45, 60])
tendance_antecedents['Courbé Droite Fort'] = fuzz.trapmf(tendance_antecedents.universe, [45, 60, 90, 90])
angle['Courbé Gauche Fort'] = fuzz.trapmf(angle.universe, [-90, -90, -60, -45])
angle['Courbé Gauche'] = fuzz.trapmf(angle.universe, [-60, -45, -20, 0])
angle['Linéaire'] = fuzz.trapmf(angle.universe, [-20, -5, 5, 20])
angle['Courbé Droite'] = fuzz.trapmf(angle.universe, [5, 20, 45, 60])
angle['Courbé Droite Fort'] = fuzz.trapmf(angle.universe, [45, 60, 90, 90])
if want_to_plot:
tendance_actuelle.view()
tendance_antecedents.view()
angle.view()
SF4_rules = [ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Linéaire'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_antecedents['Courbé Droite Fort'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Linéaire'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Droite'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Linéaire'], angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Gauche'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Linéaire'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Gauche Fort'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Gauche'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Linéaire'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Droite'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_antecedents['Courbé Droite Fort'],
angle['Courbé Droite Fort'])]
SF4_ctrl = ctrl.ControlSystem(SF4_rules)
SF4 = ctrl.ControlSystemSimulation(SF4_ctrl)
input_tendance_actuelle_fuzz = input_fuzz(tendance_actuelle, input_actu)
input_tendance_antecedents_fuzz = input_fuzz(tendance_antecedents, input_ant)
#print('Tendance calculée fuzz')
#print(input_tendance_actuelle_fuzz)
#print('Tendance évaluée fuzz')
#print(input_tendance_antecedents_fuzz)
IRR = IRR_2var(SF4_rules, input_tendance_actuelle_fuzz, input_tendance_antecedents_fuzz)
#print(IRR)
declenchement = IRR.min(axis=1)
#print(declenchement)
csq_SF4 = {} # variable de type dictionnaire
for i in range(0, len(SF4_rules)): # on initialise à 0 pour chaque règle
csq_SF4[str(SF4_rules[i].consequent)] = 0
for i in range(0, len(SF4_rules)):
csq_i = str(SF4_rules[i].consequent)
csq_SF4[csq_i] = max(csq_SF4[csq_i], declenchement[i]) # Comme dans le TP6
# print("Conséquences SF4 :", csq_SF4)
return (output_fuzz(csq_SF4, angle))
# Différence d'angle nécessaire au calcul de SF4, entre 2 points
# ________
def SF4_angle(t,dt):
return get_angle(x_TC[t],y_TC[t],x_TC[t+dt],y_TC[t+dt]) - get_angle(x_TC[t+dt],y_TC[t+dt],x_TC[t+2*dt],y_TC[t+2*dt])
def SF4_conclusion(angle_TC,angle_tendance,t):
# Prediction pour 12 heures
prediction = np.zeros(np.shape(uwind[0, :, :]))
SF4_conc = SF4_compute(angle_TC, angle_tendance,want_to_plot=False)
for lignes in range(np.shape(SF4_conc)[0]):
a = get_angle(x_TC[t],y_TC[t],x_TC[t-12],y_TC[t-12]) + SF4_conc[lignes][0]
for r in np.arange(0, 1000000, 20000):
xr, yr = polar_to_cart(r, a)
x = xr + x_TC[t]
y = - yr + y_TC[t]
idx_lon, idx_lat = xy_to_indx_lonlat(x,y, lon, lat)
prediction[idx_lat, idx_lon] = SF4_conc[lignes][1]
plot_array(prediction[:, :], lon, lat)
plt.plot(met_to_deg(x_TC[t]),met_to_deg(y_TC[t]),ms=5, marker='o', markeredgecolor='red')
plt.show()
return prediction,SF4_conc
Editor is loading...