Esercitazione 3

Taratura dinamica tramite regressione non lineare

Per la consegna dell’esercitazione si faccia riferimento a questo link. Mentre le librerie sono rimandate nell’icona Librerie.

1 Funzione di trasferimento e risposta al gradino

Funzione di trasferimento e risposta al gradino.

Il principio alla base per ricavare la funzione di trasferimento del sistema è il bilancio delle potenze termiche scambiate:

\[ m c \frac{d}{dt} T(t) = \alpha A(T_{amb}(t)-T(t)) \overset{\mathcal{F}}{\longrightarrow} mc \, j\omega T(j\omega) =\alpha A (T_{amb}(j\omega) - T(j\omega)). \]

Perciò, bilanciando la potenza scambiata tra sonda (uscita del sistema) e guaina, si ottiene:

\[\begin{aligned} m_s c_s \, j\omega T_s(j\omega) &= \alpha_{gs} A_{gs} (T_{g}(j\omega) - T_s(j\omega)) \\ \Leftrightarrow \frac{T_s(j\omega)}{T_g(j\omega)} &= \frac{1}{\tau_1 \, j\omega + 1} , \end{aligned}\]

dove \(\tau_2 = \frac{M_s c_s}{\alpha_{gs} A_{gs}}\).

Mediante gli stessi procedimenti, si può esprimere anche la potenza scambiata tra il misurando (input del sistema) e la guaina:

\[\begin{aligned} m_g c_g \, j\omega T_g(j\omega) &= \alpha_{mg} A_{gs} (T_{m}(j\omega) - T_g(j\omega)) \\ \Leftrightarrow \frac{T_g(j\omega)}{T_m(j\omega)} &= \frac{1}{\tau_2 \, j\omega + 1} , \end{aligned}\]

dove dove \(\tau_1 = \frac{M_g c_g}{\alpha_{mg} A_{mg}}\).

Pertanto, esprimendo \(T_g\) in funzione di \(T_m\), si ricava la funzione di risposta in frequenza tra misura (sonda) e misurando:

\[ H(j\omega) = \frac{1}{\tau_1 \tau_2 (j \omega)^2 + (\tau_1 + \tau_2)j\omega + 1}. \]

Riguardo alla risposta allo scalino da \(T_i\) a \(T_f\) nelle condizioni di \(\tau_1 > \tau_2\), si modella il gradino come segue, traslato di \(T_i\):

\[\begin{aligned} T_m(t) = T_i + (T_f - T_i)\,\text{sca}(t) &\overset{T_i = \text{cost}}{\Leftrightarrow} \vartheta_m(t) := T_m(t) - T_i = (T_f - T_i)\,\text{sca}(t) \\ &\overset{\mathcal{L}}{\longrightarrow} \Theta_m(s) = (T_f - T_i)\frac{1}s . \end{aligned}\]

Si definisce anche l’uscita traslata di \(T_i\):

\[ \vartheta_s(t) := T_s(t) - T_i \overset{\mathcal{L}}{\longrightarrow} \Theta_s(s). \]

In virtù della natura LTI del sistema, la risposta allo scalino nel dominio della frequenza è pari a

\[ \Theta_s(s) = H(s) \cdot (T_f - T_i)\frac{1}s. \]

Mediante la scomposizione per fratti semplici e poi anti-trasformando, si ottiene la risposta nel dominio nel tempo:

\[ T_s(t) = T_f + (T_i+T_f) \cdot \left( \frac{\tau_1}{\tau_1 - \tau_2} e^{-\frac{t}{\tau_1}} - \frac{\tau_2}{\tau_1 - \tau_2} e^{-\frac{t}{\tau_2}} \right) . \]

2 Segnale campionato

Il set di dati è disponibile a questo link.

df <- read.table("datiPT100_2005Taratura.txt", header=FALSE, sep="\t", dec=",", stringsAsFactors=FALSE)
colnames(df) <- c("t", "temp", "temp_f")

df <- df %>% dplyr::filter(t >= 150) %>% select(-temp_f) %>%
  mutate(t = t - 150)

df %>% plot_ly(x = ~t, y = ~temp) %>%
  add_lines() %>% 
  layout(xaxis = list(title="Tempo (s)"), yaxis = list(title="Temperatura (°C)"), title="Segnale di temperatura campionato")

# Verificato campionamento regolare con media e sd
fs <- 1/(df$t[2]-df$t[1]) # freq. di campionamento
Ta <- max(df$t) # tempo di acquisizione: primo istante a 0 sec

Parametri di campionamento (regolare):

• \(f_s =5.99\,\mathrm{Hz}\)

• \(T_a =244.67\,\mathrm{s}\)

3 Taratura costanti \(\tau_1\) e \(\tau_2\)

Regressione sulle sole costanti di tempo.

t0 <- 28                        # soglia temporale tra plateau e transizione
Ti <- mean(df$temp[df$t < t0])
Tf <- 82                        # attribuito visivamente

# Modello fisico: risposta allo scalino
temperature <- function(t, t0, Ti, Tf, tau1, tau2) {
  ifelse(
    t < t0,
    Ti,
    Tf + (Ti-Tf)*(
      tau1/(tau1-tau2)*exp( - (t-t0)/tau1 ) -
      tau2/(tau1-tau2)*exp( - (t-t0)/tau2 )
    )
  )
}

df.nls <- nls(temp ~ temperature(t, t0, Ti, Tf, tau1, tau2),
              data = df, start = list(tau1 = 30, tau2=25)
            )
summary(df.nls)

Formula: temp ~ temperature(t, t0, Ti, Tf, tau1, tau2)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
tau1 31.41414    0.06063  518.10   <2e-16 ***
tau2  0.86715    0.04374   19.83   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4801 on 1467 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 4.686e-06

df <- df %>% add_predictions(df.nls) %>% add_residuals(df.nls)

p1 <- df %>% ggplot(aes(x = t)) +
  geom_line(aes(y = temp)) +
  geom_line(aes( y = pred), color = "red") +
  labs(x="Tempo (s)", y="Temperatura (°C)")

p2 <- df %>% ggplot(aes(x = t, y = resid)) +
  geom_point() +
  geom_vline(xintercept = as.numeric(coef(df.nls)[1]), color = "orange", linetype = 2) +
  labs(x="Tempo (s)", y="Residui (°C)")

p1+p2

Osservando i residui, è evidente un pattern di questi durante la fase di transizione che suggerisce un underfitting dei dati. Questo potrebbe essere sintomo che il modello fisico ipotizzato non sia adeguato. Peraltro, i parametri regrediti suggeriscono un comportamento predominante dalla costante di tempo \(\tau_1\).

Tuttavia, graficamente è chiaro come la regressione nls segua bene le misure sperimentali e comunque si parla di una dispersione dei residui di circa \(\pm1 \,\mathrm{°C}\).

4 Taratura su tutti i parametri

Regressione su tutti i parametri.

params <- list(t0 = t0,
               Ti = Ti,
               Tf = Tf,
               tau1 = as.numeric(coef(df.nls)[1]),
               tau2 = as.numeric(coef(df.nls)[2])
            )

df.nls2 <- nls(temp ~ temperature(t, t0, Ti, Tf, tau1, tau2),
               data = df, start = params)
summary(df.nls2)

Formula: temp ~ temperature(t, t0, Ti, Tf, tau1, tau2)

Parameters:
     Estimate Std. Error t value Pr(>|t|)    
t0   24.89962    0.13417  185.58   <2e-16 ***
Ti   28.12231    0.03563  789.33   <2e-16 ***
Tf   81.80832    0.01983 4124.48   <2e-16 ***
tau1 30.08639    0.10480  287.10   <2e-16 ***
tau2  4.71385    0.16979   27.76   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4449 on 1464 degrees of freedom

Number of iterations to convergence: 12 
Achieved convergence tolerance: 5.334e-06

df <- df %>% add_predictions(df.nls2) %>% add_residuals(df.nls2)

p1 <- df %>% ggplot(aes(x = t)) +
  geom_line(aes(y = temp)) +
  geom_line(aes( y = pred), color = "red") +
  labs(x="Tempo (s)", y="Temperatura (°C)")

p2 <- df %>% ggplot(aes(x = t, y = resid)) +
  geom_point() +
  geom_vline(xintercept = as.numeric(coef(df.nls2)[1]), color = "orange", linetype = 2) +
  labs(x="Tempo (s)", y="Residui (°C)")

p1+p2

Anche in questo caso i residui continuano a mostrare un determinato andamento nella fase di transizione, ma visivamente la regressione si ben adatta alle misure sperimentali.

5 Intervallo di confidenza dei parametri

Intervallo di confidenza dei parametri mediante bootstrap.

L’intervallo di confidenza (IC) viene calcolato al \(95\%\): ossia è quell’intervallo al quale i parametri regrediti appartengono con tale probabilità.

# Funzione che calcola i parametri di nls ad ogni ripetizione 
# per questa ragione sono delle statistiche
stats <- function(temp, data) {
  fit <- nls(temp ~ temperature(t, t0, Ti, Tf, tau1, tau2), 
             data = data, start = params)
  return(fit$m$getPars())
}

# Bootstrap NON parametrico (con reinserimento di default)
df.b <- boot(df, \(x,i) stats(x[i, "temp"], x[i,]), R=10^4)

df.b

ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = df, statistic = function(x, i) stats(x[i, "temp"], 
    x[i, ]), R = 10^4)


Bootstrap Statistics :
     original        bias    std. error
t1* 24.899622 -0.0254079706  0.30031102
t2* 28.122310 -0.0001026519  0.03721231
t3* 81.808320 -0.0007831858  0.02444025
t4* 30.086388 -0.0176934271  0.16464437
t5*  4.713854  0.0372583421  0.38740033

È verificato che le stime regredite coincidono con la regressione precedente, per cui si accetta il calcolo senza ulteriori dimostrazioni (e.g. un grafico QQ di ciascun parametro per capire se il numero di ripetizioni \(R\) è adeguato).

# Intervallo di confidenza al 95% (metodo percentile)
ci <- list(
  t0 = boot.ci(df.b, type = "perc", index = 1)$percent[4:5],
  Ti = boot.ci(df.b, type = "perc", index = 2)$percent[4:5],
  Tf = boot.ci(df.b, type = "perc", index = 3)$percent[4:5],
  tau1 = boot.ci(df.b, type = "perc", index = 4)$percent[4:5],
  tau2 = boot.ci(df.b, type = "perc", index = 5)$percent[4:5]
)

Graficando l’IC, si mostra quanto sia incerta la stima della curva regredita.

Interessanti sono gli IC relativi a \(t_0, T_i\) e \(T_f\).

Mentre su tutta la curva, combinando i diversi IC, si ottiene:

# Funzione che calcola l'estremo alto o basso dell'IC
# per ogni possibile combinazione degli estremi di IC
# di ciascun parametro
f_conf <- function(t, temperature, ci, upper=TRUE) {
  df <- expand.grid(ci)       # ogni possibile combinazione di IC
  df$temp <- with(df, temperature(t, t0, Ti, Tf, tau1, tau2))
  ifelse(upper, max(df$temp), min(df$temp))
}

df %>% mutate(
  lower = map_dbl(t, ~f_conf(., temperature, ci, upper=FALSE)),
  upper = map_dbl(t, ~f_conf(., temperature, ci))
) %>% ggplot(aes(x=t)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill="blue", alpha=.5) +
  geom_line(aes(y=temp)) +
  geom_line(aes(y=pred), color="red") +
  labs(x="Tempo (s)", y="Temperatura (°C)")

Si osserva come l’incertezza nella fase di transizione risulti più ampia, verosimilmente riconducibile alla stima di \(\tau_1\) e \(\tau_2\). Tale comportamento può indicare una parziale inadeguatezza del modello fisico adottato oppure una difficoltà di nls nel raggiungere una soluzione ottimale.