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...