FONCTION_FLOU
unknown
plain_text
3 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[k % 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,
np.max(output_var.terms[membi].mf))
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, [-20, -20, -5, -1])
d_speed['-'] = fuzz.trimf(d_speed.universe, [-5, -1, 0])
d_speed['St'] = fuzz.trimf(d_speed.universe, [-1, 0, 1])
d_speed['+'] = fuzz.trimf(d_speed.universe, [0, 1, 5])
d_speed['++'] = fuzz.trapmf(d_speed.universe, [1, 5, 20, 20])
timeTC['Young'] = fuzz.trapmf(timeTC.universe, [0, 0, 1, 2])
timeTC['Middle'] = fuzz.trimf(timeTC.universe, [1, 2, 4])
timeTC['Old'] = fuzz.trapmf(timeTC.universe, [2.5, 4, 15, 15])
lifespan['Dying+'] = fuzz.trapmf(lifespan.universe, [0, 0, 0.5, 2])
lifespan['Dying'] = fuzz.trimf(lifespan.universe, [0.5, 2, 3])
lifespan['St'] = fuzz.trimf(lifespan.universe, [2, 3, 4])
lifespan['Living'] = fuzz.trimf(lifespan.universe, [3, 4, 8])
lifespan['Living+'] = fuzz.trapmf(lifespan.universe, [4, 8, 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['St']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Middle'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['St'] & timeTC['Old'], lifespan['Dying+']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Young'], lifespan['Living']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Middle'], lifespan['St']),
ctrl.Rule(temperature['Low'] & d_speed['+'] & timeTC['Old'], lifespan['Dying']),
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['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['-'] & 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, 400, pas), 'TC_size') # TC size (km)
TC_size_variation = ctrl.Consequent(np.arange(0, 70, 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, 400, 400])
# 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, 70]) # 25
TC_size_variation['Very High'] = fuzz.trapmf(TC_size_variation.universe, [25, 70, 70, 70]) # 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 IFT
# 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 + 0.1 * SF3_ouput
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
# EN VRAI Y A UNE SOUSTRACTION DANS LE CALCUL DE M, IL FAUDRAIT NP.FLIP LA POSITION, CEPENDANT LA DIFFERENCE
# EST EXTREMEMENT MINIME (moins de 1%), par manque de temps nous laisserons l'approximation suivante
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 = np.zeros(np.shape(psi))
for i in range(0, 4):
n[:, :, i] = scipy.ndimage.laplace(psi_large[:, :, i]) / (delta ** 2) + 1 * f[:, :, i]
# 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 / 4) * 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, np.flip(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(n_grad_y, 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)
return plot_grid
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, 0.01), 'Surface_Pressure') # Surface pressure (Pa)
dsurface_pressure = ctrl.Antecedent(np.arange(-2000, 2000, 0.01),
'Surface_Pressure_Variation') # Surface pressure variation (Pa/h) Surface
# variation per hour or per day?
velocity_variation = ctrl.Consequent(np.arange(-30, 30, 0.01),
'Velocity_Variation')
velocity['Little'] = fuzz.trapmf(velocity.universe, [0, 0, 20, 30])
velocity['Intermediate'] = fuzz.trapmf(velocity.universe, [20, 30, 30, 60])
velocity['High'] = fuzz.trapmf(velocity.universe, [30, 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, -500, -50])
dsurface_pressure['-'] = fuzz.trimf(dsurface_pressure.universe, [-500, -50, 0])
dsurface_pressure['0'] = fuzz.trimf(dsurface_pressure.universe, [-50, 0, 50])
dsurface_pressure['+'] = fuzz.trimf(dsurface_pressure.universe, [0, 50, 150])
dsurface_pressure['++'] = fuzz.trapmf(dsurface_pressure.universe, [50, 150, 2000, 2000])
velocity_variation['--'] = fuzz.trapmf(velocity_variation.universe, [-30, -30, -10, -5])
velocity_variation['-'] = fuzz.trapmf(velocity_variation.universe, [-10, -5, -1, 0])
velocity_variation['St'] = fuzz.trimf(velocity_variation.universe, [-1, 0, 1])
velocity_variation['+'] = fuzz.trapmf(velocity_variation.universe, [0, 1, 5, 15])
velocity_variation['++'] = fuzz.trapmf(velocity_variation.universe, [5, 15, 30, 30])
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['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['-'], velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['0'], velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['+'], velocity_variation['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Little'] & dsurface_pressure['++'], velocity_variation['++']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['++'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['0'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['+']),
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['High'] & dsurface_pressure['-'], velocity_variation['+']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['+'], velocity_variation['St']),
ctrl.Rule(surface_pressure['Little'] & velocity['High'] & dsurface_pressure['++'], velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['--'], velocity_variation['++']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['-'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['0'], velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['+'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Little'] & dsurface_pressure['++'], velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['--'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['-'],
velocity_variation['+']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['0'],
velocity_variation['St']),
ctrl.Rule(surface_pressure['High'] & velocity['Intermediate'] & dsurface_pressure['+'],
velocity_variation['+']),
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['High'] & dsurface_pressure['-'], velocity_variation['-']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & dsurface_pressure['0'], velocity_variation['--']),
ctrl.Rule(surface_pressure['High'] & velocity['High'] & 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
print(csq_SF3)
return output_fuzz(csq_SF3, velocity_variation)
def estimation1(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])
SF2_defuzz = (np.sum(SF2_taille_TC1[:, 0] * SF2_taille_TC1[:, 1]) / np.sum(SF2_taille_TC1[:, 1]))
print("Taille Actuelle du cyclone, km: ", np.rint(SF2_taille_TC0), "Humidité, %:", np.rint(SF2_humidity))
print("_________________________SORTIE SF2 - Future Taille du Cyclone à dt__________")
print("Prédiction :", np.round(SF2_taille_TC1, 1), 'Defuzz, km', np.rint(SF2_defuzz), 'Réel, km :',
np.rint(taille_TC(t1, 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 (mètres) =', x_TC_fuzz, 'Y (mètres) = ', y_TC_fuzz)
print('uwind.mean , m/s = ', np.rint(np.mean(uwind_fuzz, axis=(0, 1))))
print('vwind.mean , m/s = ', 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)
# UNITESSSS ??????????????
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(t0, study_grid)
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]] - pressure[
t0 - 12, idx_y_TC[t0 - 12], idx_x_TC[t0 - 12]]
print('Vitesse, m/s=', np.round(SF3_speed_TC, 1))
print('Pression, Pa =', np.round(SF3_pressure_TC, 1))
print('Variation de la pression, Pa =', np.round(SF3_dpressure_TC, 1))
SF3_output = SF3_compute(SF3_speed_TC, SF3_pressure_TC, SF3_dpressure_TC, want_to_plot=False)
print("_____________________SORTIE SF3 - Variation de vitesse______________")
SF3_defuzz = (np.sum(SF3_output[:, 0] * SF3_output[:, 1]) / np.sum(SF3_output[:, 1]))
print(np.round(SF3_output, 2), 'Defuzz :', SF3_defuzz)
x, y = CAF2_compute(uwind_fuzz, vwind_fuzz, psi, psi_TC, psi_large, x_TC_fuzz, y_TC_fuzz, SF2_defuzz * 1000,
t1 - t0,
K, SF3_defuzz)
print('Il y a ', np.size(x[:, 0]), "points !")
print('________________CAF2 - DONE________________________')
plt.figure(0)
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')
CAF2_plot = xy_fuzz_plot(x, y, lon_fuzz_plot, lat_fuzz_plot)
return CAF2_plot
def SF4_compute(input_actu, input_reg, want_to_plot):
tendance_actuelle = ctrl.Antecedent(np.arange(-90, 90, 0.1), 'Angle calcul')
tendance_regression = 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_regression['Courbé Gauche Fort'] = fuzz.trapmf(tendance_regression.universe, [-90, -90, -60, -45])
tendance_regression['Courbé Gauche'] = fuzz.trapmf(tendance_regression.universe, [-60, -45, -20, -5])
tendance_regression['Linéaire'] = fuzz.trapmf(tendance_regression.universe, [-20, -5, 5, 20])
tendance_regression['Courbé Droite'] = fuzz.trapmf(tendance_regression.universe, [5, 20, 45, 60])
tendance_regression['Courbé Droite Fort'] = fuzz.trapmf(tendance_regression.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, 25, 25, 60])
angle['Courbé Droite Fort'] = fuzz.trapmf(angle.universe, [25, 60, 90, 90])
if want_to_plot:
tendance_actuelle.view()
tendance_regression.view()
angle.view()
SF4_rules = [ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_regression['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_regression['Courbé Gauche'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_regression['Linéaire'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_regression['Courbé Droite'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche Fort'] & tendance_regression['Courbé Droite Fort'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_regression['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_regression['Courbé Gauche'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_regression['Linéaire'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_regression['Courbé Droite'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Gauche'] & tendance_regression['Courbé Droite Fort'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_regression['Courbé Gauche Fort'],
angle['Courbé Gauche Fort']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_regression['Courbé Gauche'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_regression['Linéaire'], angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_regression['Courbé Droite'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Linéaire'] & tendance_regression['Courbé Droite Fort'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_regression['Courbé Gauche Fort'],
angle['Courbé Gauche']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_regression['Courbé Gauche'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_regression['Linéaire'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_regression['Courbé Droite'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite'] & tendance_regression['Courbé Droite Fort'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_regression['Courbé Gauche Fort'],
angle['Linéaire']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_regression['Courbé Gauche'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_regression['Linéaire'],
angle['Courbé Droite']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_regression['Courbé Droite'],
angle['Courbé Droite Fort']),
ctrl.Rule(tendance_actuelle['Courbé Droite Fort'] & tendance_regression['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_regression_fuzz = input_fuzz(tendance_regression, input_reg)
IRR = IRR_2var(SF4_rules, input_tendance_actuelle_fuzz, input_tendance_regression_fuzz)
declenchement = IRR.min(axis=1)
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
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, 2000000, 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]
return prediction, SF4_conc
Editor is loading...