4. Compensazione (o Misura) dinamica (DARE)

Sommario

Esercizi sulla compensazione dinamica delle misure mediante analisi dello spettro e filtraggio offline.

1 Librerie

Quanto segue richiede questi pacchetti:

library(plotly)
library(tidyverse)
library(signal)
library(gsignal)
library(control)
library(latex2exp)

2 Funzioni di utilità

Per comodità definiamo alcune funzioni di utilità:

Bode plot mediante GGplot:

# FUNZIONE: Bodeplot
ggbodeplot <- function(tf, fmin=1, fmax=1e4, df=0.01) {
  # vector of points for each order of magnitude (OOM):
  pts <- 10^seq(0, 1, df) %>% tail(-1)
  # vector of OOMs:
  ooms <- 10^(floor(log10(fmin)):ceiling(log10(fmax)-1))
  # combine pts and ooms:
  freqs <- as.vector(pts %o% ooms)
  # warning: bode wants pulsation!
  bode(tf, freqs*2*pi) %>% {
    tibble(f=.$w/(2*pi), `magnitude (dB)`=.$mag, `phase (deg)`=.$phase)} %>%
    pivot_longer(-f) %>% 
    ggplot(aes(x=f, y=value)) +
    geom_line() +
    scale_x_log10(
      minor_breaks=scales::minor_breaks_n(10), 
      labels= ~ latex2exp::TeX(paste0("$10^{", log10(.), "}$"))) +
    facet_wrap(~name, nrow=2, scales="free") +
    labs(x="frequency (Hz)")
}

Generazione di un segnale sinusoidale con diverse armoniche, come definite in una tabella (pars) con le colonne w (ampiezza), f (frequenza) e phi (fase) per le varie frequenze (una riga per ogni armonica):

signal <- function(t, pars, rad = FALSE) { 
  stopifnot(is.data.frame(pars))
  with(pars, {
    if (!rad) {
      phi <- phi/180*pi
      f <- 2*pi*f
    }
    map_dbl(t, \(t) sum( map_vec(seq_along(w) , \(i) w[i]*sin(t*f[i] + phi[i] ))))
  })
}

3 Simuliamo l’uscita di un sistema di misura

# Parametri del misurando (input del sistema di misura)
pars <- tibble(
  w = c(1, 0.3, 0.3), 
  f = c(1/30, 1/10, 1/4), 
  phi = c(0, 0, 0)
)

# Generiamo il misurando u
# Specifica il tempo di campionamento:
Ts <- 0.01  
u <- tibble(
  t = 0:10000 * Ts,
  u = signal(t, pars),
  un = u + rnorm(length(t), 0, pars$w[1]/10)
)

# Definiamo la Funzione di Trasferimento dello strumento:
# Coefficienti del numeratore e denominatore
num <- c(5, 1)   # 5s + 1
den <- c(1, 2, 1) # s^2 + 2s + 1
H <- tf(num, den)
print(H)

y1:
    5 s^1 + 1 
  - - - - - - - - - -
     s^2 + 2 s + 1 


Transfer Function: Continuous time model 

# Diagramma di Bode
ggbodeplot(H, fmin=0.01, fmax=1e2, df=1/max(u$t))

(u %>% 
  mutate(
    out = lsim(H, un, t)$y[1,],
    ideal = lsim(H, u, t)$y[1,],
  ) %>% 
  pivot_longer(c(un, out, ideal), names_to = "series") %>% 
  ggplot(aes(x=t, y=value, color=series)) + 
  geom_line()) %>% 
  ggplotly()

# Simulazione dell'uscita con rumore
output <- lsim(H, u$un, u$t)

# Simulazione dell'uscita ideale
outputI <- lsim(H, u$u, u$t)
yI <- outputI$y[1,]

# Grafici
pp <- plot_ly() %>%
  add_lines(u$t, output$y, name = "output", line= list(color= "red")) %>%
  add_lines(u$t, u$un, name = "misurando", line= list(color= "blue")) %>%
  add_lines(u$t, yI, name = "output ideale", line= list(color= "green"))
# pp

Si nota come un rumore ‘bianco’, ovvero contenente tutte le armoniche, viene ridotto in uscita dal sistema di misura che si comporta in modo simile ad un passa-basso.

4 Stima del misurando a partire dal segnale in uscita

Se vogliamo invece stimare il misurando a partire dal segnale in uscita invertendo quindi la funzione di trasferimento, condizione che si verifica per le misure ‘dinamiche’ (per le statiche ricordate si inverte la caratteristica statica che si ottiene dopo taratura, statica). In questo caso, cosa succederà a parità di rumore bianco in uscita?

# Segnale in uscita con rumore
yn = yI + rnorm(length(u$t), 0, pars$w[1]/10)


# Ricaviamo l'inversa della Funzione di Trasferimento:
Hinv <- tf(den, num)

[1] "TFCHK: Transfer function may not be proper and may lead to errors. Num > Den"

print(Hinv)

y1:
     s^2 + 2 s + 1 
  - - - - - - - - - -
    5 s^1 + 1 

Transfer Function: Continuous time model 

# Diagramma di Bode
ggbodeplot(Hinv, fmin=0.01, fmax=1e2, df=1/max(u$t))

# SCOMMENTARE:
# # Simulazione dell'ingresso
# input <- lsim(Hinv, yn, y$t, x0 = rep(0, length(pole(Hinv))))
# 
# # Grafici dell'ingresso e dell'uscita
# pp <- plot_ly() %>%
#   add_lines(y$t, y$yn, name = "output", line= list(color= "blue")) %>%
#   add_lines(y$t, input$y, name = "Estimated input", line= list(color= "red"))
# pp

Otteniamo l’errore: “The order of the Numerator should be equal or lesser than the Denominator” Dunque R si rifiuta di impiegare la funzione di trasferimento inversa. Perchè? Trovare una soluzione e svilupparla continuando il codice fornito.

Esercitazione

Riconoscimento punti esame

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

Assegnazioni

Suggerimenti