Visto che puro ML non funziona per le serie tempo di cui mi sto occupando ed le regressioni basate su formule analitiche mostrano dei limiti in presenza del rumore degli strumenti proviamo con PINN (Physics informed Neural Networks) in cui si conosce a priori la legge che regola il fenomeno e si demanda alla rete neurale di approssimare al meglio le variabili sconosciute o di gestire il rumore strumentale
Data la mancanza di dati mi sono fatto (...ha fatto ChatGPT) un generatore di dati sintetici che puo' aggiungere un livello arbitrario di rumore
Sono stati creati due dataset
 |
| Dataset prova 1 |
 |
| Dataset prova 2 |
import csv
import random
import matplotlib.pyplot as plt
def generate_landslide_data(num_points=100, initial_displacement=10.0, velocity_decrease_rate=0.005, noise_level=0.01):
"""
Generates simulated landslide data with a *linear* decrease in inverse velocity,
and calculates inverse velocity. Correctly applies noise to the velocity.
"""
data = []
inverse_velocity = 1.0 # Start with an arbitrary inverse velocity
displacement = initial_displacement
previous_displacement = initial_displacement #set initial displacement
for day in range(1, num_points + 1):
# Linearly decrease the inverse velocity:
inverse_velocity -= velocity_decrease_rate
# Ensure inverse velocity doesn't go below a small positive value to avoid div by zero
inverse_velocity = max(inverse_velocity, 0.0001)
# Calculate velocity (V = 1 / (1/V)):
velocity = 1 / inverse_velocity
# Add noise to the velocity *directly*:
noise = random.uniform(-noise_level * velocity, noise_level * velocity)
velocity += noise
# Ensure velocity doesn't become negative after adding noise (important!)
velocity = max(velocity, 0.0001) #Very small positive value to avoid zero velocity.
inverse_velocity = 1/velocity #Recalculate inverse velocity based on the noisy velocity.
# Calculate displacement increase (ΔDisplacement = Velocity * ΔTime)
displacement_increase = velocity
# Calculate new displacement
displacement += displacement_increase
data.append((day, displacement, inverse_velocity))
return data
def write_to_csv(data, filename="landslide_data.csv"):
"""
Writes the simulated landslide data to a CSV file with semicolon as a separator,
including a column for inverse velocity.
"""
with open(filename, "w", newline="") as csvfile:
writer = csv.writer(csvfile, delimiter=";")
writer.writerow(["Day", "Cumulative_Displacement", "Inverse_Velocity"])
for row in data:
inverse_velocity_str = "{:.6f}".format(row[2]) if row[2] is not None else ""
writer.writerow([row[0], row[1], inverse_velocity_str])
print(f"Data written to {filename}")
def plot_inverse_velocity(filename="landslide_data.csv"):
"""
Reads the landslide data from the CSV file and plots inverse velocity vs. time using matplotlib.
"""
days = []
inverse_velocities = []
with open(filename, "r") as csvfile:
reader = csv.reader(csvfile, delimiter=";")
next(reader) # Skip the header row
for row in reader:
try:
day = int(row[0])
inverse_velocity = float(row[2]) if row[2] else None #Handle empty cells
if inverse_velocity is not None:
days.append(day)
inverse_velocities.append(inverse_velocity)
except (ValueError, IndexError) as e:
print(f"Error reading row: {row}. Skipping. Error: {e}")
# Plotting the Data:
if days and inverse_velocities: # Check to make sure data was read.
plt.figure(figsize=(10, 6)) # Set plot size
plt.plot(days, inverse_velocities, marker='o', linestyle='-', color='blue') #Plot data
plt.title("Fukuzono Plot: Inverse Velocity vs. Time") #Set title
plt.xlabel("Time (Days)") #set labels
plt.ylabel("Inverse Velocity (days/mm)")
plt.grid(True) #Add grid
plt.show() #Display plot
else:
print("No data to plot.")
# --- Main execution ---
if __name__ == "__main__":
landslide_data = generate_landslide_data(num_points=100, initial_displacement=10.0, velocity_decrease_rate=0.011, noise_level=0.05) #Tuned parameters
write_to_csv(landslide_data)
plot_inverse_velocity() #Call the plotting function
A questo punto ho provato a farci girare la PINN
Il metodo di Fukuzono dell'inverso della velocita' deve essere espresso in modo differenziale
La formula completa senza la semplificazione dell'esponente 2 (le lettere tra le due formule indicano grandezze non coincidenti)
con la rete che cerca di ottimizzare il parametro m
 |
| modello dataset 1 |
 |
| modello dataset 2 |
alla fine si vede che il modello ha sovrastimato il momento del collasso
come nota il modello ha bisogno di un mumero elevato (> 5000 epochs) per convergere in maniera minimamente soddisfacente
import tensorflow as tf
import numpy as np
import pandas as pd
# Carica i dati dal file CSV (Day, Cumulative_Displacement, Inverse_Velocity)
# Modifica qui con il nome del file
file_path = './fuku/landslide_data3.csv'
df = pd.read_csv(file_path, delimiter=';')
t = df['Day'].values.reshape(-1, 1) # Tempo (giorni)
inv_v = df['Inverse_Velocity'].values.reshape(-1, 1) # Inverso della velocità
# Parametro iniziale m (da ottimizzare)
m_init = tf.Variable(1.5, dtype=tf.float32, trainable=True)
# Modello PINN
class PINN(tf.keras.Model):
def __init__(self):
super(PINN, self).__init__()
self.hidden1 = tf.keras.layers.Dense(20, activation='tanh')
self.hidden2 = tf.keras.layers.Dense(20, activation='tanh')
self.hidden3 = tf.keras.layers.Dense(20, activation='tanh')
self.out = tf.keras.layers.Dense(1)
def call(self, t):
x = self.hidden1(t)
x = self.hidden2(x)
x = self.hidden3(x)
return self.out(x)
# Inizializza il modello
model = PINN()
def loss_fn(t, inv_v, m):
t = tf.convert_to_tensor(t, dtype=tf.float32) # Conversione a tensore
inv_v = tf.convert_to_tensor(inv_v, dtype=tf.float32)
with tf.GradientTape() as tape:
tape.watch(t)
u = model(t)
du_dt = tape.gradient(u, t)
if du_dt is None:
du_dt = tf.zeros_like(t) # Evita errori se il gradiente è None
eqn = du_dt - m * ((u*u) / tf.convert_to_tensor(t_f - t, dtype=tf.float32)) # Converte t_f - t in tensore
mse_data = tf.reduce_mean(tf.square(u - inv_v))
mse_phys = tf.reduce_mean(tf.square(eqn))
return mse_data + mse_phys
# Ottimizzatore
optimizer = tf.keras.optimizers.Adam(learning_rate=0.001)
# Tempo di collasso iniziale (da stimare)
t_f = tf.Variable(np.max(t) + 1.0, dtype=tf.float32, trainable=True)
# Training Loop
epochs = 10000
for epoch in range(epochs):
with tf.GradientTape() as tape:
loss = loss_fn(t, inv_v, m_init)
grads = tape.gradient(loss, model.trainable_variables + [t_f, m_init])
optimizer.apply_gradients(zip(grads, model.trainable_variables + [t_f, m_init]))
if epoch % 10 == 0:
print(f"Epoch {epoch}, Loss: {loss.numpy()}, t_f: {t_f.numpy()}, m: {m_init.numpy()}")
from scipy.interpolate import interp1d
from scipy.optimize import brentq
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
# Compute min and max time from dataset
t_min = np.min(t)
t_max = np.max(t)
# Generate predictions for a fine grid of time values
t_pred = np.linspace(t_min, t_max, 500).reshape(-1, 1).astype(np.float32)
inv_v_pred = model(t_pred).numpy().flatten()
# Ensure `t_pred_original` has the correct shape
t_pred_original = t_pred.flatten()
# Fit a linear model to estimate collapse time
lin_reg = LinearRegression()
lin_reg.fit(t_pred_original.reshape(-1, 1), inv_v_pred)
slope = lin_reg.coef_[0]
intercept = lin_reg.intercept_
# Compute collapse time as x-intercept (1/V = 0)
if slope != 0:
collapse_time = -intercept / slope
print(f"Estimated collapse time: {collapse_time:.2f} days")
else:
collapse_time = None
print("Could not estimate collapse time (slope is zero).")
# Visualization
plt.figure(figsize=(10, 6))
# Scatter plot of real data
plt.scatter(t, inv_v, label='Data', color='red', marker='o')
# PINN-predicted curve
plt.plot(t_pred_original, inv_v_pred, label='PINN Prediction', color='blue')
# Linear regression line (dashed black)
t_reg_line = np.linspace(t_min, t_max, 500)
inv_v_reg_line = lin_reg.predict(t_reg_line.reshape(-1, 1))
plt.plot(t_reg_line, inv_v_reg_line, linestyle='dashed', color='black', label="Linear Fit")
# Vertical line at estimated collapse time
if collapse_time is not None and t_min <= collapse_time <= t_max:
plt.axvline(collapse_time, color='green', linestyle='dashed', label=f'Collapse Time: {collapse_time:.2f} days')
plt.xlabel('Time (days)')
plt.ylabel('Inverse Velocity (1/V)')
plt.legend()
plt.grid()
plt.title(f"Forecasted Collapse Time: {collapse_time:.2f}" if collapse_time else "Forecasted Collapse Time: Not Found")
plt.show()
