L’idea alla base dell’algoritmo a differenze finite stocastiche è di utilizzare un campionamento stocastico per stimare la derivata di una funzione rispetto a una variabile di input, invece di utilizzare la sua esatta derivata analitica.Questo approccio consente di ottenere una stima più robusta della derivata, in presenza di incertezza e rumore nei dati.

Operativamente il primo passaggio è stato quello di definire una funzione NaK da poter derivare ciclicamente. La funzione richiede due argomenti di input, malto e luppolo, che vengono combinati in un vettore input. Successivamente,viene chiamata la funzione birra sull’input combinato, il cui risultato viene assegnato alla variabile output.

La variabile L viene quindi calcolata come la somma del primo e del secondo elemento del vettore output, che rappresentano i livelli di sodio e potassio nella birra.

source("birra.R")
NaK=function(malto,luppolo){
  input=c(malto,luppolo)
  output=birra(input)
  L=output[1]+output[2]
  return(L)
}
set.seed(123)

Estremi del range dei due parametri:

range1x<-50
range2x<-100
range1y<-10
range2y<-20

Numero massimo di iterazioni:

nmaxiter<-1000

Parametri del modello a differenze stocastiche finite:

alpha<-0.5
gamma<-0.2
a<-100
A<-80
c<-1

Guess iniziale dei parametri assegnato in modo casuale:

malto0<-runif(1,range1x,range2x) 
luppolo0<-runif(1,range1y,range2y)

Algoritmo

Il codice seguente utilizza un ciclo while per iterare fino a un numero massimo di iterazioni (1.000).

In ogni iterazione, l’algoritmo approssima le derivate parziali della funzione di perdita NaK rispetto ai parametri di input malto e luppolo utilizzando le differenze finite, i valori ottenuti vengono assegnati alle variabili g1 e g2.

In seguito, l’algoritmo calcola i nuovi valori dei parametri di input malto e luppolo utilizzando il parametro ak e i gradienti calcolati.

ak viene calcolato come una funzione dell’iterazione corrente niter e di due costanti a e A. I nuovi valori dei parametri di input malto e luppolo sono calcolati come una combinazione lineare dei valori correnti e dei gradienti calcolati, entro i rispettivi limiti definiti in precedenza.

Come ultimo passo, i nuovi valori di malto e luppolo vengono assegnati alle nuove variabili malto0 e luppolo0 per la prossima iterazione, e il valore corrente della funzione di perdita NaK viene salvato in un vettore L, che rappresenta la Loss Function.

niter<-0
L<-numeric()
df=data.frame()
while (niter<nmaxiter){
  niter<-niter+1
  ck<-c/(niter+1)^gamma
  g1<-(NaK(malto0+ck,luppolo0)-NaK(malto0-ck,luppolo0))/(2*ck)
  g2<-(NaK(malto0,luppolo0+ck)-NaK(malto0,luppolo0-ck))/(2*ck)
  ak<-a/(niter+1+A)^alpha
  malto1<-min(max(malto0-ak*g1,range1x),range2x)
  luppolo1<-min(max(luppolo0-ak*g2,range1y),range2y)
  malto0<-malto1
  luppolo0<-luppolo1
  L[niter]<-NaK(malto0,luppolo0)
  parametri=c(malto0,luppolo0,L[niter])
  df=rbind(df,parametri)
}

colnames(df)=c("malto","Luppolo","Loss")

DataFrame risultante

Statistiche di sintesi:

summary(df) 
##      malto           Luppolo           Loss      
##  Min.   : 50.00   Min.   :10.00   Min.   :29.77  
##  1st Qu.: 56.92   1st Qu.:10.00   1st Qu.:32.76  
##  Median : 77.39   Median :12.00   Median :34.25  
##  Mean   : 75.81   Mean   :14.31   Mean   :34.36  
##  3rd Qu.: 95.19   3rd Qu.:20.00   3rd Qu.:35.82  
##  Max.   :100.00   Max.   :20.00   Max.   :39.83

Valore di ottimo

In questo caso l’ottimo non si trova in corrispondenza dell’ultima iterazione a causa del rumore della funzione birra.

Nell’esperimento in esame il minimo pari a 29,77 mg/L è stato ottenuto alla 621esima iterazione in corrispondenza di 98,11 g/L di malto e 10,62 g/L di luppolo.

df[L==min(L),]
##        malto  Luppolo     Loss
## 621 98.10728 10.61726 29.77087

Grafico della funzione

Dopo aver fatto il plot dei valori della funzione di loss per tutte le 1.000 iterazioni viene evidenziato il punto di minimo con il comando points().

plot(df$Loss,type="p",main="Funzione di Loss",xlab="Iterazioni",ylab="Sodio+Potassio") 
minimo=df[L==min(L),]
points(y=minimo$Loss,x=rownames(minimo),col="red",pch=19)