3. Taratura Dinamica tramite Regressione non lineare

taratura
regressione
Autore/Autrice

Mariolino De Cecco e Paolo Bosetti

Data di Pubblicazione

15 maggio 2025

Sommario

Esempio per la taratura dinamica di un sensore PT100 mediante regressione lneare e non-lineare ai minimi quadrati.

Attenzione: 1. Scriviamo dplyr::filter(i>80) anziché semplicemente filter(i>80) in quanto la funzione si trova in due pacchetti con finalità diverse

1 Taratura dinamica

Consideriamo il caso di una sonda PT100 a temperatura \(T_i\) immersa repentinamente in un mezzo termostatato a temperatura \(T_f\). Questa azione la possiamo modellare come un ingresso ideale a gradino in quanto la dinamica termica, essenso lenta, è praticamente insensibile al lasso di tempo in cui si inserisce la sonda nel mezzo. Sotto queste ipotesi la temperatura della sonda raggiunge quella del mezzo secondo la legge:

\[ T(t) = (T_i - T_f)e^{-\frac{t}{\tau}} + T_f \tag{1}\]

Leggiamo un file di taratura acquisito in laboratorio consistente nei dati ottenuti inserendo una PT100 in un fornelletto termostatato, simulando quindi l’ingresso a gradino appena modellato matematicamente.

# Leggiamo il file assumendo che sia separato da tabulazioni (\t)
dati <- read.table("0_Dati/gio_1_3_n1bis.txt", header=FALSE, sep="\t", dec=",", stringsAsFactors=FALSE)
colnames(dati) <- c("Tempo", "Temp", "Temp_f")
dati[] <- lapply(dati, function(x) as.numeric(gsub(",", ".", x)))

s <- tibble(
  t = dati[,1],
  Temp = dati[,2]
)

# Grafici
pp <- plot_ly(s) %>%  # 's' è il dataframe
  add_lines(x = ~t, y = ~Temp, name = "Temperatura", line = list(color = "red")) %>%
  layout(
    title = "Grafico Temperatura vs Tempo",
    xaxis = list(title = "Tempo (s)"),
    yaxis = list(title = "Temperatura (°C)")
  )
pp
Con GGPlot2 
(s %>% 
  ggplot(aes(x=t, y=Temp)) +
  geom_line() +
   labs(x="Tempo (s)", y="Temperatura (°C)")) %>% 
  ggplotly()

1.1 Secondo metodo: linearizzazione

La Equazione 1 può essere resa lineare nei coefficienti così:

\[ \begin{align} \frac{T(t) -T_f}{T_i - T_f} &= e^{-\frac{t}{\tau}} \\ \ln\left(\frac{T(t) -T_f}{T_i - T_f}\right) &= \ln(e^{-\frac{t}{\tau}}) \\ \ln\left(\frac{T(t) -T_f}{T_i - T_f}\right) &= -\frac{t}{\tau} \end{align} \]

Per ottenere tale andamento lineare possiamo quindi riorganizzare i dati come segue ed eseguire un fitting lineare con la funzione lm():

# Definizione dei valori iniziali e finali
Ti <- 33
Tf <- 95
ti <- 46

# Selezione e trasformazione dei dati
si <- s %>% 
  select(t, Temp) %>% 
  mutate(t = t - ti) %>%
  dplyr::filter(t > 0 ) %>%
  mutate(y = log((Temp-Tf)/(Ti - Tf))) %>% 
  dplyr::filter(!is.nan(y))
Warning: There was 1 warning in `mutate()`.
ℹ In argument: `y = log((Temp - Tf)/(Ti - Tf))`.
Caused by warning in `log()`:
! NaNs produced
# Fitting Lineare
si.lm <- si %>% 
  dplyr::filter(t< 70) %>%
  lm(y~t, data=.) 
si.lm %>% summary()

Call:
lm(formula = y ~ t, data = .)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.11338 -0.01097 -0.00022  0.01401  0.07064 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)  5.262e-02  3.071e-03   17.14   <2e-16 ***
t           -2.928e-02  7.605e-05 -384.98   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02558 on 277 degrees of freedom
Multiple R-squared:  0.9981,    Adjusted R-squared:  0.9981 
F-statistic: 1.482e+05 on 1 and 277 DF,  p-value: < 2.2e-16
gain = si.lm$coefficients[2]
tau <- round(-1/gain, 1)
cat(paste("tau:", tau))
tau: 34.2
# Grafici
pp <- plot_ly(si) %>% 
  add_lines(x = ~t, y = ~y, name = "Andamento stimato", line = list(color = "red")) %>%
  add_lines(x = ~t, y = ~gain*t, name = "Log Temperatura", line = list(color = "blue")) %>%
  layout(
    title = "Grafico Linearizzato Temperatura vs Tempo",
    xaxis = list(title = "Tempo [s]"),
    yaxis = list(title = "Temperatura Linearizzata [°C]")
  )
pp
Con GGPlot2 
library(modelr)

Attaching package: 'modelr'
The following object is masked from 'package:gsignal':

    resample
pp <- si %>% 
  add_predictions(si.lm) %>% 
  select(t, y, pred) %>% 
  pivot_longer(-t) %>% 
  ggplot(aes(x=t, y=value, color=name)) +
  geom_line() +
  labs(x="Tempo (s)", y="Temperatura (°C)")
ggplotly(pp)

1.2 Terzo metodo: regressione non-lineare

Con il terzo metodo si usa la regressione non-lineare ai minimi quadrati per ottenere \(\tau\) o direttamente tutti i possibili parametri incogniti \(T_i, T_f, \tau, t_i\). Solitamente però la sonda PT100 costituia da filo di Platino avvolto a spirale è protetta da una guaina metallica che assorbe e scambia calore con il sensore e con l’ambiente di cui ha il compito di misurare la temperatura. Il modello più realistico è leggermente più complicato.

Figura: Modello 2° ord di una PT100

Con \(\alpha\) ed A sono indicati il coefficiente di conduzione e l’area di contatto tra guaina e sonda e tra misurando e guaina, la massa m ed il calore specifico c della sonda e della guaina sono indicati rispettivamente con i pedici S e G.

In questo caso il modello è del secondo ordine e la funzione di trasferimento tra la trasformata del misurando \(T_M(\omega)\) e la trasformata dell’uscita \(T_{S}(\omega\) diviene:

\[ \frac{T_{S}(\omega)}{T_M(\omega)} = \frac{1}{\tau_1 \tau_2 i\omega^{2} + (\tau_1 + \tau_2)i\omega + 1} \]

Un sistema termico è generalmente sovrasmorzato per cui le radici dell’equazione caratteristica sono reali e distinte. Definendo le costanti di tempo come \(\tau_1\) > \(\tau_2\), la risposta per una variazione a gradino da \(T_i\) a \(T_f\) seguirà l’andamento temporale:

\[ T(t) = T_f + (T_i - T_f) * (\frac{\tau_1}{\tau_1 - \tau_2} e^{-t/\tau_1} - \frac{\tau_2}{\tau_1 - \tau_2} e^{-t/\tau_2}) \]

Esercitazione

Riconoscimento punti esame

Il completamento di questo esercizio comporta il riconoscmiento di ±0.5 punti.

Assegnazioni

Elenco assegnazioni

  • A partire dal modello del sistema con guaina mostrato in figura,
    mostrare i passaggi necessari per ottenere la funzione di trasferimento
    e la risposta al gradino

  • Impiegando i dati forniti, eseguire la regressione non lineare in cui i
    parametri incogniti sono solo \(\tau_1\) e \(\tau_2\)

  • Impiegando i dati forniti, eseguire la regressione non lineare in cui i
    parametri incogniti sono \(T_i, T_f, \tau_1, \tau_2, t_i\)

  • Mediante il metodo bootstrap, ottenere gli intervalli di confidenza sui
    parametri

Suggerimenti

2 OMETTERE E FAR FARE COME ESERCIZIO

# FUNZIONE: Generazione risposta esponenziale del 1° ordine
ModelloPrimoOrdine <- function(t, Ti, Tf, tau){ 
    map_dbl(t, \(t)  (Ti - Tf)*exp(-t/tau) + Tf)
}

# FUNZIONE: Generazione risposta esponenziale del 2° ordine
ModelloSecondoOrdine <- function(t, Ti, Tf, tau1, tau2){
    map_dbl(t, \(t)  Tf + (Ti - Tf)*( (tau1 / (tau1-tau2) ) * exp(-t/tau1) - (tau2 / (tau1-tau2) ) * exp(-t/tau2)) )
}


# Selezione e trasformazione dei dati
si <- s %>% 
  select(t, Temp) %>% 
  mutate(t = t - 46) %>%
  filter(t > 0)

# Definizione dei valori iniziali
Ti <- 33
Tf <- 95

# Fit del modello di secondo ordine
fit <- nls(Temp ~ ModelloSecondoOrdine(t, Ti, Tf, tau1, tau2), 
    data = si, 
    start = list(
      tau1 = 2,    # Prima costante di tempo
      tau2 = 3      # Seconda costante di tempo
    ),
    algorithm = "port",  # Usa un algoritmo più robusto
    lower = list(tau1 = 0.1, tau2 = 0.0001),  # Limiti inferiori
    upper = list(tau1 = 100, tau2 = 10)  # Limiti superiori 
)

# Estrazione dei coefficienti stimati
tau1 <- coef(fit)["tau1"]
tau2 <- coef(fit)["tau2"]

# Visualizzazione dei risultati
summary(fit)

Formula: Temp ~ ModelloSecondoOrdine(t, Ti, Tf, tau1, tau2)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
tau1  32.0515     0.1901  168.62   <2e-16 ***
tau2   2.8109     0.1426   19.71   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.252 on 1117 degrees of freedom

Algorithm "port", convergence message: relative convergence (4)
# 1° ordine
# Ti <- 33
# Tf <- 95
# fit <- nls(Temp ~ ModelloPT100(t, Ti, Tf, tau), 
#     data = si, 
#     start = list(
#       tau=10
#     ))
# 
# fit
# tau <- coef(fit)["tau"]
# tau

# Inseriamo l'andamento stimato da modello nella struttura dati:
si <- si %>% mutate(Temp_fit = ModelloSecondoOrdine(t, Ti, Tf, tau1, tau2))
print(si)
# A tibble: 1,119 × 3
       t  Temp Temp_fit
   <dbl> <dbl>    <dbl>
 1  0.25  33.2     33.0
 2  0.5   33.7     33.1
 3  0.75  32.5     33.2
 4  1     33.1     33.3
 5  1.25  32.9     33.5
 6  1.5   34.0     33.6
 7  1.75  33.3     33.8
 8  2     33.5     34.1
 9  2.25  34.3     34.3
10  2.5   34.3     34.6
# ℹ 1,109 more rows
# Grafici
pp <- plot_ly(si) %>%  # 'si' è il dataframe
  add_lines(x = ~t, y = ~Temp, name = "Temperatura", line = list(color = "blue")) %>%
  add_lines(x = ~t, y = ~Temp_fit, name = "Temperatura \n prevista \n da modello", line = list(color = "red")) %>%
  layout(
    title = "Grafico Temperatura vs Tempo",
    xaxis = list(title = "Tempo (s)"),
    yaxis = list(title = "Temperatura (°C)")
  )
pp

3 Junks

#| include: false
tau <- 33.5
Ti <- 33.4
Tf <- 95
t0 <- 50
temp <- tibble(
  t = seq(0, 350, 0.5),
  Tn = ifelse(t<t0, Ti, (Ti - Tf)*exp(-(t-t0)/tau) + Tf),
  T = Tn + rnorm(length(t), 0, 0.5)
)

Sappiamo che temperatura iniziale \(T_i = 35\) °C, \(T_f = 95\) °C e stimiamo che l’immersione inizi a \(t_0 = 50\) s.

ggplot(s, aes(x=t, y=Temp)) + 
  geom_line() +
  labs(x="tempo (s)", y="temperatura (°C)")

Acquisizione gradino con termocoppia PT100