POLITECNICO DI MILANO
CORSO DI LAUREA IN INGEGNERIA MATEMATICA
ORIENTAMENTO CALCOLO SCIENTIFICO
REGOLAZIONE DELLA CONVEZIONE
FORZATA PER RAFFREDAMENTO
AD ARIA DI UN CHIP
PROGETTO DI
ANALISI NUMERICA DELLE EQUAZIONI A DERIVATE PARZIALI 2
PROGRAMMAZIONE AVANZATA PER IL CALCOLO SCIENTIFICO
AZZURRA PORPORA Matr.
FRANCESCA STORTI Matr.
Anno Accademico 2008-2009
708835
708131
Indice
1 Introduzione
3
2 Problema in esame
2.1 Equazioni del modello . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.2 Parametri del modello . . . . . . . . . . . . . . . . . . . . . . . . . . . .
2.3 Analisi di buona posizione . . . . . . . . . . . . . . . . . . . . . . . . . .
4
4
6
8
3 Problema di controllo ottimo
3.1 Analisi di esistenza e unicità . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Risoluzione numerica: metodo iterativo . . . . . . . . . . . . . . . . . .
11
15
21
4 Implementazione del codice
4.1 Pre-processing e mesh computazionale . . . .
4.2 Integrazione numerica . . . . . . . . . . . . .
4.3 Assemblaggio e risoluzione dei sistemi lineari
4.4 Ciclo iterativo . . . . . . . . . . . . . . . . . .
4.5 Post-processing . . . . . . . . . . . . . . . . .
25
25
28
31
36
37
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
5 Risultati numerici
38
6 Conclusioni
46
Appendice - Codici Freefem
47
Riferimenti Bibliografici
59
2
1
Introduzione
In questo lavoro ci proponiamo di calcolare il valore ottimale di regolazione di una
ventola, che ha lo scopo di mantenere sotto una certa soglia la temperatura su di un chip.
Quest’ultimo dispositivo funziona infatti soltanto all’interno di un certo range di valori
di temperatura: il superamento di tali valori può causare un malfunzionamento del chip
stesso fino a portarlo a rottura. Lo strumento matematico utilizzato per raggiungere tale
scopo è la teoria del controllo ottimo per equazioni alle derivate parziali. L’andamento
della temperatura viene modellato tramite un’equazione di diffusione con opportune
condizioni al contorno, mentre il parametro di cui si vuole trovare il valore ottimale
viene modellato tramite un funzionale costo ad esso associato.
In particolare, il sistema fisico è costituito da tre componenti: il chip, affiancato da
un dissipatore passivo di calore e dalla ventola. Il dissipatore, posto a diretto contatto
con il chip, permette di disperdere nell’ambiente circostante il calore ricevuto dal chip
tramite flusso conduttivo. Il processo dissipativo avviene però sia per convezione libera
che forzata: il lavoro del dissipatore è infatti coadiuvato dalla ventola, che accelera il
raffredddamento dell’intero sistema.
Trovare il punto di funzionamento ottimale della ventola significa trovare la sua
potenza e la sua regolazione ideali, in base alla temperatura di lavoro del chip e alla capacità di dispersione termica del dissipatore. L’EDP considerata modella la temperatura
all’interno del dissipatore passivo; il funzionale costo da minimizzare riguarda invece la
temperatura sul chip, che non deve discostarsi troppo dal suo valore di funzionamento
ideale; infine, il parametro di cui vogliamo trovare il valore ottimale è il coefficiente
di convezione forzata, direttamente collegato al punto di funzionamento ottimale della
ventola.
Nella sezione 2 presentiamo il modello fisico considerato con tutti i parametri ad esso
associati, insieme all’equazione di diffusione che descrive l’andamento della temperatura
ed alla sua analisi di buona posizione. Nella sezione 3 presentiamo il vero e proprio
problema di controllo ottimo: descriviamo l’espressione del funzionale costo e ricaviamo
le equazioni che legano tale funzionale all’EDP per la temperatura; presentiamo inoltre
l’analisi di esistenza e unicità della variabile di controllo ottimo. Nella sezione 4 presentiamo diversi metodi numerici necessari per la risoluzione del problema di controllo
ottimo, in particolare descriviamo il ciclo iterativo che permette il raggiungimento del
valore ottimale della variabile di controllo. Nella sezione 5 presentiamo il codice da noi
sviluppato tramite il linguaggio di programmazione C++ per la risoluzione del problema
di controllo. Inizialmente abbiamo però implementato il problema tramite FreeFem++,
in modo da poter effettuare un confronto tra i risultati ottenuti con i due diversi programmi; i programmi FreeFem++ sono riportati in Appendice. Infine nella sezione 6
riportiamo e commentiamo i principali risultati ottenuti.
3
2
2.1
Problema in esame
Equazioni del modello
Per studiare la propagazione del calore e il campo di temperatura all’interno del dissipatore, consideriamo i diversi scambi termici che avvengono nell’intero sistema. All’interno
del dissipatore il calore si propaga tramite conduzione termica. La distribuzione di temperatura nel dissipatore, in regime stazionario, è quindi descritta dalla seguente legge di
diffusione:
−k∆T = 0 ,
dove k è il coefficiente di conducibilità termica e T è la temperatura.
Il calore generato dal chip si propaga nel dissipatore attraverso conduzione termica: la
quantità di calore scambiata dipende dalla geometria e dalle caratteristiche del corpo e
dalla differenza di temperatura tra i due diversi mezzi. In regime stazionario il flusso termico, ovvero la potenza termica per unità di superficie, trasmesso dal chip al dissipatore
per conduzione è definito dalla legge di Fourier per la conduzione:
q = −k∇T ,
dove q indica il flusso termico.
Infine il dissipatore cede calore all’ambiente esterno tramite convezione termica. Si
può distinguere tra due diversi tipi di convezione: convezione naturale se lo scambio di
calore avviene per una semplice differenza di temperatura e convezione forzata se invece
lo scambio di calore è favorito dall’azione di dispositivi esterni che generano trasporto
di massa, come ad esempio la ventola nel caso in esame. Il flusso termico legato alla
convezione è direttamente proporzionale alla differenza di temperatura, come espresso
dalla legge di Newton per la convezione:
q = h(T − T∞ ) ,
dove h è il coefficiente di convezione termica e T∞ è la temperatura del fluido (nel
nostro caso aria) a distanza sufficientemente grande dal dissipatore. Il coefficiente h
non è una proprietà del fluido, ma è un parametro determinato sperimentalmente, il cui
valore dipende da tutte le variabili che influenzano il processo di convezione, quali la
geometria della superficie di scambio, la natura del moto, le proprietà e la velocità del
fluido circostante. Tale parametro può essere legato alla semplice convezione naturale
oppure può includere il contributo di un’eventuale convezione forzata (come nel nostro
caso).
Consideriamo ora uno schema semplificato e bidimensionale del dissipatore in esame
(si veda la figura 1), dove con Ω indichiamo il dominio spaziale, ovvero l’intero volume
del dissipatore, con Γ1 la superficie a contatto con l’aria circostante, interessata alla
convezione, con Γ2 la superficie a contatto con il supporto solido del dispositivo, considerata adiabatica, cioè che non permette la trasmissione del calore verso l’ambiente
esterno (attraverso la quale il flusso di calore q è nullo) e con Γ3 la superficie a contatto
con il chip, interessata alla conduzione. Si osservi la presenza di superfici alettate, al
4
fine di aumentare la superficie di contatto tra il dissipatore e l’aria circostante, cosı̀ da
aumentare la quantità di calore scambiato per convezione. Il chip è fissato alla superficie
orizzontale inferiore del dissipatore e da esso si propaga la potenza termica.
Figura 1: Schema semplificato del dissipatore.
Il problema trattato viene quindi modellato tramite il seguente sistema di equazioni:


−∇ · (k∇T ) = 0
in Ω




 k∇T · n = −(h + u)(T − T∞ )
su Γ1
,
(1)

k∇T · n = 0
su Γ2




 k∇T · n = q
su Γ3
dove T = T (x, y) è la temperatura, funzione a valori scalari delle due variabili spaziali x
e y, k è la conducibilità termica dell’alluminio, materiale di cui è costituito il dissipatore
in esame, q è il flusso di calore tra chip e dissipatore, n è il vettore normale esterno alle
diverse parti del dominio, T∞ è la temperatura dell’ambiente circostante, (h + u) è il
coefficiente di convezione forzata dell’aria (infatti, essendo h, coefficiente di convezione
naturale dell’aria, molto piccolo rispetto ad u, coefficiente di convezione forzata, possiamo accorpare i due coefficienti in modo da definirne uno singolo globale).
La prima equazione del sistema (1) è un’equazione di semplice diffusione e descrive la
distribuzione di temperatura all’interno del dissipatore; la seconda è una condizione di
Robin e definisce il flusso di calore scambiato per convezione tra il dissipatore e l’aria
circostante; la terza e la quarta sono condizioni di Neumann e definiscono rispettivamente il flusso di calore nullo attraverso la superficie adiabatica del dissipatore e il flusso
di calore scambiato per conduzione tra il dissipatore e il chip. Si tratta quindi di un
5
problema ellittico con condizioni di bordo miste, di cui studieremo le proprietà di buona
posizione (esistenza e unicità della soluzione).
2.2
Parametri del modello
Vediamo ora in dettaglio il calcolo dei parametri fisici presenti nel modello espresso dal
sistema (1).
Il chip considerato è il tipo T-5000 della Intel, mentre il dissipatore è una versione
semplificata del tipo HS 6115 della Vemaline Products, Inc. ([Ç98, pag. 360]): esso è un
dispositivo anodizzato nero sagomato di alluminio, costituito da 10 alette di lunghezza
25 mm, larghezza 76 mm e spessore 3 mm. La lunghezza totale del dispositivo è di
102 mm.
Le principali proprietà termodinamiche del sistema chip-dissipatore sono evidenziate
dal coefficiente di conducibilità termica del dissipatore k e dal coefficiente di convezione
naturale dell’aria h.
La conducibilità termica k è quella propria dell’alluminio, cioè k = 237 mW◦ C .
La convezione naturale h riferita all’aria dipende dalla particolare geometria del dispositivo considerato:
karia
h = 1.31
,
(2)
Sopt
dove karia = 0.026 mW◦ C è la conducibilità termica dell’aria e Sopt indica la spaziatura
ottimale tra le alette del dissipatore. Tale spaziatura può essere calcolata tramite la
seguente formula:
L
Sopt = 2.714
,
(3)
Ra
dove L è la lunghezza caratteristica del dissipatore, corrispondente nel nostro caso
alla larghezza delle alette (76 mm), e Ra è il numero di Rayleigh (adimensionale).
Quest’ultimo viene a sua volta calcolato per definizione con la seguente formula:
Ra =
gγ(Ts − T∞ )δ 3
Pr ,
ν2
(4)
dove:
• g = 9.81
m
s2
è l’accelerazione di gravità,
• γ è il coefficiente di dilatazione cubica misurato in K −1 ,
• Ts e T∞ , misurate entrambe in ◦ C, sono rispettivamente la temperatura sulla
superficie del dispositivo e la temperatura dell’aria sufficientemente lontano dalla
superficie,
• δ è una lunghezza caratteristica, coincidente in questo caso con L,
2
• ν = 1.57 × 10−5 ms è la viscosità cinematica dell’aria,
• P r = 0.712 è il numero di Prandtl (adimensionale).
6
Consideriamo T∞ coincidente con la temperatura ambiente, pari cioè a 20◦ C, e Ts pari
alla temperatura massima di funzionamento del chip, cioè 90◦ C; γ può essere calcolato
come l’inverso della temperatura media tra T∞ e Ts misurate in K, cioè nel seguente
modo:
2
γ=
.
(5)
(Ts + T∞ )
Il numero di Rayleigh calcolato con la formula (4) risulta quindi pari a: Ra = 0.849×106 .
Inserendo questo valore nella formula (3) per calcolare Sopt , si ottiene: Sopt = 6.8 mm.
A questo punto possiamo calcolare tramite la formula (2) la convezione naturale dell’aria
h, che risulta pari a:
W
h = 5.7 2◦ .
m C
Il coefficiente di convezione forzata dell’aria u, che è direttamente collegato alle caratteristiche termodinamiche della ventola, è il parametro che possiamo variare in modo
che il chip non si riscaldi troppo; il suo valore è solitamente compreso tra 25 mW
2◦ C e
250 mW
.
2◦ C
Il flusso termico q che passa dal chip al dissipatore può essere calcolato tramite una
semplice equazione di bilancio, che eguaglia il flusso di calore prodotto dal chip al flusso
di calore disperso dal dissipatore Integrando la prima equazione del sistema (1) su tutto
il dominio Ω e sfruttando il teorema della divergenza otteniamo:
Z
Z
k∇T · ndΓ = 0 ,
− ∇ · (k∇T )dΩ = −
∂Ω
Ω
che può essere riscritto separatamente sulle tre diverse porzioni di bordo del dominio
come:
Z
Z
Z
−
k∇T · ndΓ −
k∇T · ndΓ −
k∇T · ndΓ = 0 ,
Γ1
Γ2
Γ3
e sostituendo le condizioni al contorno del sistema (1):
Z
Z
(h + u)(T − T∞ )dΓ −
qdΓ = 0 .
Γ1
Γ3
Quest’ultima equazione può essere riscritta nel seguente modo:
Z
Z
qdΓ =
(h + u)(T − T∞ )dΓ ,
Γ3
Γ1
e ricordando che le quantità q e (h + u) sono costanti, possiamo ricavare l’espressione
per il calcolo del flusso di calore q:
Z
(h + u)
q=
(T − T∞ )dΓ ,
|Γ3 |
Γ1
dove |Γ3 | rappresenta la lunghezza del bordo Γ3 . Scegliamo a questo punto un valore
ragionevole per il coefficiente di convezione forzata pari a (h̄ + ū) = 25 mW
2◦ C e un valore
7
medio per la temperatura pari a T̄ = 90◦ C: il valore T̄ − T∞ diventa cosı̀ una costante,
per cui l’integrale può essere riscritto nel seguente modo:
q=
(h̄ + ū)
(T − T∞ ) |Γ1 | ,
|Γ3 |
dove |Γ1 | rappresenta la lunghezza del bordo Γ1 . Sostituendo i valori dei parametri
presenti nell’ultima equazione otteniamo infine il valore del flusso termico q:
q = 15080
2.3
W
.
m2
Analisi di buona posizione
Ci occupiamo ora dell’analisi di buona posizione del problema ellittico di partenza
espresso dal sistema (1), prima di affrontare il problema di controllo ottimo vero e
proprio.
Gli spazi funzionali che utilizzeremo sono: H 1 (Ω) e R, in particolare T ∈ H 1 (Ω) e
u ∈ R. Prima di affrontare il problema di controllo vero e proprio, verifichiamo la buona
posizione del problema ellittico di partenza.
Ricaviamo innanzitutto la formulazione debole del problema (1). Consideriamo quindi una funzione test v ∈ H 1 (Ω) e moltiplichiamo l’equazione differenziale del sistema
(1) per tale funzione:
Z
−
∇ · (k∇T )vdΩ = 0 .
Ω
Applichiamo la regola di integrazione per parti, ottenendo:
Z
Z
−
k∇T · nvdΓ +
k∇T · ∇vdΩ = 0 .
∂Ω
Ω
Sostituendo le condizioni al contorno del sistema (1) e ricordando che h, u e q sono
grandezze costanti, otteniamo:
Z
Z
Z
(h + u)
(T − T∞ )vdΓ − q
vdΓ +
k∇T · ∇vdΩ = 0 ,
Γ1
Γ3
Ω
da cui, portando al secondo membro gli integrali che non contengono l’incognita T :
Z
Z
Z
Z
k∇T · ∇vdΩ + (h + u)
T vdΓ = (h + u)T∞
vdΓ + q
vdΓ .
Ω
Γ1
Γ1
A questo punto definiamo la forma bilineare a : H 1 (Ω) × H 1 (Ω) → R:
Z
Z
a(T, v) =
k∇T · ∇vdΩ + (h + u)
T vdΓ ,
Ω
Γ3
(6)
Γ1
e il funzionale lineare F : H 1 (Ω) → R:
Z
F (v) = (h + u)T∞
Z
vdΓ + q
Γ1
8
vdΓ .
Γ3
(7)
La formulazione debole del problema (1) risulta quindi essere:
trovare T ∈ H 1 (Ω) : a(T, v) = F (v)
∀v ∈ H 1 (Ω) .
(8)
A questo punto per dimostrare l’esistenza e l’unicità della soluzione del problema scritto
nella forma debole (8), utilizziamo il lemma di Lax-Milgram:
Lemma 2.1. Sia V uno spazio di Hilbert, a(·, ·) : V × V → R una forma bilineare
continua e coerciva. Sia F (·) : V → R un funzionale lineare e continuo. Allora esiste
unica la soluzione del problema:
trovare u ∈ V : a(u, v) = F (v)
∀v ∈ V .
Verifichiamo ora le ipotesi di questo lemma.
Per dimostrare la continuità della forma bilineare a(·, ·), prendiamo in considerazione i due membri della forma bilineare separatamente e applichiamo ad entrambi la
disuguaglianza di Schwarz :
Z
k∇T · ∇vdΩ ≤ k k∇T k 2 k∇vk 2
L (Ω)
L (Ω) ,
Ω
Z
(h + u)
T vdΓ ≤ (h + u) kT kL2 (Γ1 ) kvkL2 (Γ1 ) ,
Γ1
ricordando che i parametri k, h ed u sono costanti e positivi. Utilizziamo poi la
definizione di norma nello spazio H 1 (Ω) per il primo membro e la disuguaglianza di
traccia (con costante C) per il secondo:
Z
k∇T · ∇vdΩ ≤ k kT k 1 kvk 1
H (Ω)
H (Ω) ,
Ω
Z
(h + u)
T vdΓ ≤ C 2 (h + u) kT kH 1 (Ω) kvkH 1 (Ω) .
Γ1
Dalle ultime due disequazioni si ricava quindi:
|a(T, v)| ≤ max k, C 2 (h + u) kT kH 1 (Ω) kvkH 1 (Ω) .
(9)
La forma bilineare è dunque continua.
Per dimostrare la coercività della forma bilineare a(·, ·), utilizziamo il lemma di
Petree-Tartar (si veda in proposito [EG04]):
Lemma 2.2. Siano X, Y , Z tre spazi di Banach. Siano inoltre A : X → Y un operatore
lineare, continuo e iniettivo e B : X → Z un operatore lineare, continuo e compatto. Se
esiste una costante c > 0 tale che c kxkX ≤ kAxkY + kBxkZ , allora esiste una costante
c̄ tale che:
c̄ kxkX ≤ kAxkY
∀x ∈ X .
Per i nostri scopi possiamo riformulare il lemma nella forma seguente:
9
Lemma 2.3. Sia Ω un dominio sufficientemente regolare (ad esempio tale che ∂Ω sia
lipschitziano) e sia f : H 1 (Ω) → R un funzionale lineare tale che f (v) 6= 0
∀v 6= 0,
v ∈ R. Allora esiste una costante positiva cR tale che:
cR kvkH 1 (Ω) ≤ k∇vkL2 (Ω) + |f (v)| .
(10)
Dalla (10), detta anche disuguaglianza di Poincarè-Friedrichs, elevando al quadrato e
sfruttando la disuguaglianza di Young si ottiene:
c2R kvk2H 1 (Ω) ≤ (k∇vkL2 (Ω) + |f (v)|)2
= k∇vk2L2 (Ω) + |f (v)|2 + 2 k∇vkL2 (Ω) |f (v)|
≤ 2 k∇vk2L2 (Ω) + |f (v)|2 .
(11)
Per trattare il nostro problema di partenza, ovvero il caso di equazione di diffusione con
condizione di Robin e Neumann, scegliamo opportunamente come funzionale f (v) =
(h+u)1/2
k1/2
kvk2L2 (Γ1 ) . Riscriviamo quindi la (11) utilizzando l’espressione del funzionale:
h+u
2
2
2
2
kvkL2 (Γ1 ) ,
cR kvkH 1 (Ω) ≤ 2 k∇vkL2 (Ω) +
k
che equivale alla seguente:
k k∇vk2L2 (Ω) + (h + u) kvk2L2 (Γ1 ) ≥
k 2
c kvk2H 1 (Ω) .
2 R
(12)
Dalla definizione della forma bilineare (6) ricaviamo che:
a(v, v) = k k∇vk2L2 (Ω) + (h + u) kvk2L2 (Γ1 ) ,
che combinata con la (12) ci fornisce la seguente:
a(v, v) ≥
k 2
c kvk2H 1 (Ω) .
2 R
(13)
La forma bilineare è dunque coerciva.
Per dimostrare la continuità del funzionale F (·), utilizziamo nell’ordine la disuguaglianza triangolare e la disuguaglianza di Cauchy-Schwarz :
Z
Z
|F (v)| ≤ (h + u)T∞
vdΓ + q
vdΓ
Γ1
Γ3
≤ (h + u)T∞ kvkL2 (Γ1 ) + q kvkL2 (Γ3 ) ,
ricordando che i parametri h, u, T∞ e q sono costanti e positivi. Applichiamo infine la
disuguaglianza di traccia (con costanti C e C̃) e otteniamo:
h
i
|F (v)| ≤ (h + u)T∞ C + q C̃ kvkH 1 (Ω) .
(14)
Il funzionale è dunque continuo.
Grazie alle (9), (13) e (14), sono quindi soddisfatte le ipotesi del lemma di LaxMilgram, da cui segue l’esistenza e l’unicità della soluzione del problema (1) di partenza.
10
3
Problema di controllo ottimo
Lo studio della distribuzione di temperatura all’interno del dissipatore tramite il problema (1), è finalizzato al controllo della temperatura raggiunta dal chip, cioè quella
relativa al bordo Γ3 (superficie di contatto tra il chip e il dissipatore). Questa temperatura non può infatti superare un certo valore di soglia stabilito a priori, oltre cui non
è più garantito il corretto funzionamento del dispositivo. Per tenere conto di questa
condizione, introduciamo il seguente funzionale costo, che rappresenta la differenza tra
la temperatura effettiva T a cui si trova il chip e la sua temperatura di lavoro ideale
media Td :
Z
1
1
J(u) =
(T − Td )2 dΓ + βu2 ,
(15)
2 Γ3
2
dove il primo termine valuta appunto il quadrato della differenza tra le due temperature
T e Td sulla superficie Γ3 di contatto tra il dissipatore e il chip, mentre il secondo
termine pesato con un opportuno coefficiente β (positivo e sufficientemente piccolo)
e contenente la variabile di controllo u, è un termine opportunamente aggiunto per
garantire la coercività del problema, condizione necessaria per la buona posizione dello
stesso.
Unendo il sistema di partenza (1) all’espressione del funzionale costo (15) si ottiene un
problema di controllo ottimo. In termini generali per controllo si intende la possibilità
di intervenire sullo stato di un sistema con lo scopo di influenzarne il comportamento.
Da un punto di vista matematico un problema di controllo è caratterizzato da:
• la variabile di controllo u, scelta in un opportuno spazio funzionale, in generale può
essere definita sul dominio Ω o solo su una parte di bordo Γ; nel nostro caso u ∈ R,
è definito solo sulla porzione di bordo Γ1 e dal punto di vista fisico rappresenta il
parametro di convezione forzata;
• le variabili che descrivono lo stato del sistema, dette funzioni di stato, definite in
un opportuno spazio funzionale e dipendenti dal valore assunto dal controllo u; nel
nostro caso si tratta della temperatura T ∈ H 1 (Ω);
• l’insieme di equazioni e relative condizioni al bordo che costituiscono il cosiddetto
problema di stato, il quale descrive il sistema fisico che si intende controllare e
che deve essere soddisfatto dalle variabili di stato e di controllo; in forma generale
l’equazione di stato può essere scritta come:
AT (u) = f + Bu ,
(16)
dove A è un operatore differenziale che comprende anche le condizioni al bordo e
che agisce sullo stato T , f è il termine noto e B è un altro operatore che introduce
il controllo u nell’equazione di stato; nel nostro caso si fa riferimento al sistema
(1);
• un funzionale obiettivo (o costo) che si vuole minimizzare, definito nel nostro caso
dalla (15);
11
• la funzione di osservazione, dipendente dalla variabile di stato e quindi a sua volta
dal controllo, definita anch’essa in un opportuno spazio funzionale, va confrontata
con un certo valore desiderato stabilito a priori che rappresenta l’obiettivo da
raggiungere attraverso la minimizzazione del funzionale costo; nel nostro caso si
tratta della stessa temperatura T valutata sulla porzione di bordo Γ3 .
Il nostro è quindi un problema di controllo ottimo con controllo e osservazione al bordo.
Esso si traduce in un problema di minimizzazione vincolata: trovare un controllo u e
uno stato T che minimizzino il funzionale J(u) sotto il vincolo dato dalle equazioni
del problema di stato (1). Per la risoluzione di tale problema utilizziamo il metodo
della lagrangiana: attraverso l’uso dei moltiplicatori di Lagrange ci si riconduce ad un
problema di ottimizzazione libera, si cercano cioè i punti critici liberi del funzionale
lagrangiano definito in generale nel modo seguente:
L(T, p, u) = J(u) + hp, f + Bu − AT (u)i ,
(17)
dove p ∈ H 1 (Ω) rappresenta il moltiplicatore di Lagrange (o variabile aggiunta), il
simbolo h·, ·i denota la dualità tra H 1 (Ω) e H −1 (Ω) e si è fatto riferimento alla forma generale dell’equazione di stato (16). Il funzionale lagrangiano equivale quindi al
funzionale costo con l’aggiunta di un termine, formalmente nullo, che rappresenta il
prodotto tra la variabile aggiunta e il vincolo dato dall’equazione di stato espressa in
forma implicita. Nel nostro caso otteniamo:
L(T, p, u) = J(u) + hp, ∇ · (k∇T )i .
(18)
Il problema di ottimizzazione libera risulta quindi essere:
trovare (T, p, u) : ∇L(T, p, u) = 0 ,
(19)

L =0


 T
Lp = 0 ,



Lu = 0
(20)
esprimibile anche come:
dove si è usata la notazione abbreviata LT per indicare la derivata di Gateaux del funzionale L rispetto alla variabile T ; analogo significato hanno le notazioni Lp e Lu . Le tre
equazioni del sistema (20) forniranno rispettivamente l’equazione aggiunta, l’equazione
di stato e la condizione di ottimalità, mediante le quali determineremo lo stato, la variabile aggiunta e il controllo ottimi.
Vediamo ora, secondo due diversi approcci, come ottenere queste tre equazioni che
costituiscono il vero e proprio problema di controllo.
Consideriamo dapprima l’approccio utilizzato in [Q08].
Il funzionale lagrangiano in forma debole, considerando cioè la versione debole della
(17), diventa:
L(T, p, u) = J(u) − c(T, p) ,
12
dove c(T, p) indica la forma bilineare associata all’operatore differenziale A e dove i
termini legati al termine noto f e all’operatore B non sono presenti, in quanto nulli
nel nostro problema. Nel nostro caso il funzionale lagrangiano in forma debole risulta
quindi essere:
Z
1
1
L(T, p, u) =
(T − Td )2 dΓ + βu2
2 Γ3
2
Z
Z
Z
pdΓ .
(21)
(T − T∞ )pdΓ + q
k∇T · ∇pdΩ − (h + u)
−
Γ3
Γ1
Ω
A questo punto il problema di controllo diventa:
trovare (T, p, u) ∈ (H 1 (Ω) × H 1 (Ω) × R) :
∇L(T, p,u)[ϕ, φ, ψ] = 0
∀(ϕ, φ, ψ) ∈ (H 1 (Ω) × H 1 (Ω) × R) ,
traducibile, tramite le derivate di Gateaux,

L =0


 T
Lp = 0



Lu = 0
nel sistema di tre equazioni:
∀ϕ ∈ H 1 (Ω)
∀φ ∈ H 1 (Ω)
.
(22)
∀ψ ∈ R
Calcolando le tre equazioni del sistema (22), annullando cioè le derivate di Gateaux
del funzionale lagrangiano espresso dalla (21) rispetto alle tre variabili (considerate in
questo caso indipendenti tra loro), otteniamo la formulazione cercata per il problema di
controllo. L’equazione di stato in forma debole risulta:
Z
Z
Z
Lp [φ] = − k∇T · ∇φdΩ − (h + u)
(T − T∞ )φdΓ + q
φdΓ = 0
∀φ ∈ H 1 (Ω) .
Ω
Γ1
Γ3
(23)
L’equazione aggiunta in forma debole risulta:
Z
Z
Z
LT [ϕ] = − k∇ϕ · ∇pdΩ − (h + u)
ϕpdΓ +
(T − Td )ϕdΓ = 0
Ω
Γ1
∀ϕ ∈ H 1 (Ω) ,
Γ3
che, cambiando i segni e integrando per parti il primo termine, diventa:
Z
Z
Z
LT [ϕ] = − ϕ∇ · (k∇p)dΩ +
kϕ∇p · ndΓ + (h + u)
ϕpdΓ
Ω
∂Ω
Γ1
Z
−
(T − Td )ϕdΓ = 0
∀ϕ ∈ H 1 (Ω) . (24)
Γ3
Dalla (24) si ricava quindi il problema aggiunto:

−∇ · (k∇p) = 0





 k∇p · n + (h + u)p = 0
k∇p · n = 0





 k∇p · n − (T − Td ) = 0
13
in Ω
su Γ1
su Γ2
su Γ3
.
(25)
Infine l’equazione che esprime il principio di minimo (o condizione di ottimalità) risulta:
Z
(T − T∞ ) pdΓ
Lu [ψ] =βuψ − ψ
Γ1
Z
= βu −
(T − T∞ )pdΓ ψ ≥ 0
∀ψ ∈ R .
(26)
Γ1
Le equazioni (23), (24) e (26) sono quelle che andranno risolte per risolvere il problema
di ottimo.
L’approccio basato sul funzionale lagrangiano può essere sviluppato anche in un
altro modo, partendo direttamente dalla formulazione forte del problema (si veda [S07]).
Consideriamo l’equazione di stato come un vincolo su T e u:
G(T, u) = −∇ · (k∇T ) = 0 ,
dove G : H 1 (Ω) × R(Γ1 ) → H −1 (Ω). Per tenere conto di tale vincolo, scriviamo il
funzionale costo nella seguente forma:
Z
Z
1
1
J(T ) =
(T − Td )2 dΓ + βu2 +
p∇ · (k∇T )dΩ ,
2 Γ3
2
Ω
dove abbiamo aggiunto al funzionale costo originale (15) un termine praticamente nullo.
La variabile p appare come un moltiplicatore a nostra disposizione, che scegliamo dipendente solo dalla variabile di controllo u. Considerando questo approccio, la condizione
di ottimalità J 0 (u)(v − u) ≥ 0 diventa:
Z
Z
(T (u) − Td ) (T (v) − T (u)) dΓ+βu(v−u)+ p∇·(k∇ (T (v) − T (u))) dΩ ≥ 0 . (27)
Γ3
Ω
Integriamo due volte per parti l’ultimo termine della (27):
Z
p∇ · (k∇(T (v) − T (u))) dΩ
Ω
Z
Z
=
kp∇(T (v) − T (u)) · n dΓ −
k∇p · ∇(T (v) − T (u))dΩ
Z∂Ω
ZΩ
=
kp∇(T (v) − T (u)) · n dΓ −
k(T (v) − T (u))∇p · n dΓ
∂Ω
∂Ω
Z
+
k(T (v) − T (u))∆p dΩ .
(28)
Ω
Consideriamo ora il primo termine della (28) e sostituiamo in esso le condizioni al
contorno del problema di partenza (1):
Z
Z
Z
kp∇(T (v) − T (u)) · n dΓ =
h(T (u) − T (v))dΓ +
T∞ (v − u)dΓ
∂Ω
Γ1
Γ1
Z
+
(uT (u) − vT (v))dΓ .
(29)
Γ1
14
Riscriviamo il termine non lineare (l’ultimo della (29)) nel seguente modo:
Z
Z
Z
(uT (u) − vT (v))dΓ =
(u − v)T (u)dΓ +
u(T (u) − T (v))dΓ
Γ1
Γ1
Γ1
Z
+
(v − u)(T (u) − T (v))dΓ ,
(30)
Γ1
in modo da separare le variabili u e v da T (u) e T (v), rispettivamente. Raggruppando
tutti i termini trovati, la condizione di ottimalità (27) diventa:
Z
(T (u) − Td )(T (v) − T (u))dΓ + βu(v − u)
Γ3
Z
Z
T∞ (v − u)dΓ
h(T (u) − T (v))dΓ +
+
Γ1
Γ1
Z
Z
Z
+
(u − v)T (u)dΓ +
u(T (u) − T (v))dΓ +
(v − u)(T (u) − T (v))dΓ
Γ1
Γ1
Γ1
Z
Z
−
k(T (v) − T (u))∇p · ndΓ +
k∆p(T (v) − T (u))dΩ ≥ 0 .
(31)
∂Ω
Ω
Da quest’ultima disequazione ricaviamo il problema aggiunto:

−k∆p = 0
in Ω





 k∇p · n + (h + u)p = 0
su Γ1
k∇p · n = 0





 k∇p · n − (T (u) − Td ) = 0
su Γ2
,
su Γ3
e il principio di minimo:
Z
βu −
(T − T∞ )p dΓ (v − u) ≥ 0
∀v ∈ R .
Γ1
Questi coincidono con i risultati trovati con l’approccio lagrangiano in forma debole.
Notiamo che è possibile esprimere il principio di minimo anche in un’altra forma:
Z
J 0 (u)(v − u)dΓ ≥ 0
∀v ∈ R ,
Γ1
J 0 (u)
dove
indica la derivata del funzionale costo; tale forma è quella che viene effettivamente utilizzata nella risoluzione numerica del problema di controllo.
3.1
Analisi di esistenza e unicità
Il nostro problema di controllo ottimo in forma forte può essere considerato come un
problema di minimo vincolato (si veda [S07]), esprimibile nel caso generale come:

 min F (x, u)
con il vincolo:
(32)

G(x, u) = 0 ,
15
dove F : X × U → R è un funzionale della forma F (x, u) = P (x) + N (u) e G è un altro
funzionale tale che G : X × U → Y . X, Y e U sono spazi di Banach riflessivi e U è
un covesso chiuso (quindi anche debolmente chiuso). Diciamo che (x̂, û) è soluzione del
problema (57) se è un punto di minimo vincolato per F . In particolare, vale il seguente
teorema di esistenza:
Teorema 3.1. Assumiamo che:
1. inf x∈X P = p0 ≥ −∞ ed esiste (x,u) tale che G(x,u)=0.
2. N (u) ≥ β||u||γU , con β, γ > 0, per ogni u ∈ U .
3. Se uk * u in U e xk * x in X, dove (xk , uk ) ∈ X ×U , allora G(xk , uk ) * G(x, u)
in Y .
4. F è debolmente semicontinuo inferiormente in X × U .
5. Se (xk , uk ) ∈ X × U è una successione minimizzante ed è tale che |P (xk )| ≤ M e
G(xk , uk ) = 0, allora ||xk ||X ≤ M0 .
Allora esiste una soluzione (x̂, û) ∈ X × U del problema (57).
Nel caso specifico che abbiamo preso in considerazione, gli spazi funzionali con cui
lavoriamo sono: X = H 1 (Ω), U = R(Γ1 ) e Y = H −1 (Ω). La variabile x corrisponde alla
temperatura T e la variabile u rappresenta la variabile di controllo. Il funzionale F (T, u)
coincide con il funzionale
da minimizzare J(T (u), u), che può essere scomposto nelle due
R
parti P (T ) = 21 Γ3 (T − Td )2 dΓ e N (u) = 12 βu2 . Il vincolo G(T, u) è l’equazione ellittica
da cui partiamo: G(T, u) = −∇ · (k∇T ).
Per dimostrare l’esistenza di una soluzione (T̂ , û) applichiamo il teorema (4.1),
verificando la validità di tutte le sue ipotesi.
1. Dal
R momento 2che la funzione integranda è elevata al quadrato, il funzionale P (T ) =
1
2 Γ3 (T − Td ) dΓ è sempre non negativo, quindi il suo valore minimo è raggiunto
quando la funzione T −Td si annulla, cioè quando T = Td . Di conseguenza abbiamo
che:
inf
T ∈H 1 (Ω)
P = 0 ≥ −∞.
Inoltre, dall’analisi di buona posizione del problema (10), sappiamo che l’equazione
G(T, u) = 0 ammette soluzione unica.
2. La seconda ipotesi richiede la coercività del funzionale N , facilmente dimostrabile
grazie alla presenza del parametro β > 0:
R
N (u) = 12 β Γ1 u2 dΓ = β|u|2 ,
avendo inglobato il fattore 21 all’interno del parametro β e avendo indicato con | · |
la norma in R(Γ1 ), cioè il valore assoluto. Dall’ultima equazione si vede quindi
che vale l’uguaglianza, con γ = 2.
16
3. Consideriamo le due successioni debolmente convergenti:
uk * u in R(Γ1 )
Tk * T in H 1 (Ω),
dove (Tk , uk ) ∈ H 1 (Ω) × R(Γ1 ). Vogliamo dimostrare che:
G(Tk , uk ) * G(T, u) in H −1 (Ω),
cioè che:
R
R
− Ω ∇ · (k∇Tk )vdΩ → − Ω ∇ · (k∇T )vdΩ
∀v ∈ H 1 (Ω).
Integriamo per parti due volte:
R
− Ω ∇ · (k∇Tk )vdΩ
e sostituiamo le condizioni al contorno:
Z
Z
Z
∂Tk
k
k∇Tk ∇vdΩ −
− ∇ · (k∇Tk )vdΩ =
vdΓ
∂n
∂Ω
Ω
Ω
Z
Z
Z
∂Tk
∂v
kTk ∆vdΩ −
k
vdΓ
=
kTk dΓ −
∂n
∂n
ZΩ
Z∂Ω
Z∂Ω
∂v
kTk ∆vdΩ +
(h + uk )(Tk − T∞ )vdΓ+
=
kTk dΓ −
∂n
Ω
Γ1
∂Ω
Z
−
qvdΓ .
Γ3
Consideriamo adesso gli ultimi quattro integrali separatamente:
R
R
∂v
∂v
• ∂Ω kTk ∂n
dΓ → ∂Ω kT ∂n
dΓ
per la convergenza debole della successione {Tk }.
R
R
• − Ω kTk ∆vdΩ → − Ω kT ∆vdΩ
per la convergenza debole della successione {Tk }.
R
• Γ1 (h + uk )(Tk − T∞ )vdΓ può essere scomposto in altri quattro integrali:
R
R
(a) Γ1 hTk vdΓ → Γ1 hT vdΓ
per la convergenza debole della successione {Tk }.
R
(b) − Γ1 hT∞ vdΓ è una costante e rimane inalterata nel passaggio al limite
per k → ∞.
R
(c) Γ1 uk Tk vdΓ è una forma trilineare nelle variabili u, T e v.
R
(d) − Γ1 uk T∞ vdΓ è una costante e rimane inalterata nel passaggio al limite
per k → ∞.
Analizziamo più in dettaglio la forma trilineare (c); essa può essere riscritta
nel seguente modo:
Z
Z
(c)
uk Tk vdΓ =
(uk Tk ± uTk ± uT )vdΓ
Γ1
Γ1
Z
=
[(uk − u)Tk + (Tk − T )u + T u)]vdΓ
Γ1
17
e, dal momento che uk * u in R(Γ1 ) e Tk * T in H 1 (Ω):
R
(uk − u)Tk vdΓ → 0
RΓ1
Γ1 (Tk − T )uvdΓ → 0,
quindi:
R
R
T
u
vdΓ
→
k
k
Γ1
Γ1 T uvdΓ.
R
• − Γ3 qvdΓ è una costante e rimane inalterata nel passaggio al limite per
k → ∞.
In conclusione:
Z
Z
Z
Z
∂v
kT
hT vdΓ+
kT ∆vdΩ +
− ∇ · (k∇Tk )vdΩ →
dΓ −
∂n
∂Ω
Γ1
Ω
Ω
Z
Z
Z
Z
−
hT∞ vdΓ +
hT∞ vdΓ −
uT vdΓ −
Γ1
Γ1
qvdΓ ,
Γ3
Γ1
che, controintegrando per parti, significa:
R
R
− Ω ∇ · (k∇Tk )vdΩ → − Ω ∇ · (k∇T )vdΩ,
cioè:
G(Tk , uk ) * G(T, u) in H −1 (Ω).
4. In generale, F (x, u) si dice debolmente semicontinuo inferiormente in un punto
(x0 , u0 ) ∈ X × U se:
F (x0 , u0 ) = P (x0 ) + N (u0 ) ≤ lim inf P (x) + lim inf N (u)
x→x0
u→u0
e F (x, u) si dice debolmente semicontinuo inferiormente in X × U se è debolmente
semicontinuo inferiormente in ogni suo punto (x0 , u0 ) ∈ X × U .
Nel nostro caso, dato (T0 , u0 ) ∈ H 1 (Ω) × R(Γ1 ), possiamo dire che:
F (T0 , u0 ) = P (T0 ) + N (u0 )
Z
1
1
=
(T0 − Td )2 dΓ + βu20
2 Γ3
2
Z
1
1
lim (T − Td )2 dΓ + β lim u2
=
2 Γ3 T →T0
2 u→u0
Z
1
1
=
lim inf (T − Td )2 dΓ + β lim inf u2 .
2 Γ3 T →T0
2 u→u0
A questo punto possiamo applicare il lemma di Fatou al funzionale P (T ) valutato
in T = T0 , dal momento che (T − Td )2 è una successione di funzioni misurabili a
valori non negativi:
Z
1
1
lim inf (T − Td )2 dΓ + β lim inf u2 ≤
2 Γ3 T →T0
2 u→u0
18
1
≤ lim inf
2 T →T0
Z
1
(T − Td )2 dΓ + β lim inf u2 ,
2 u→u0
Γ3
quindi:
F (T0 , u0 ) = P (T0 ) + N (u0 ) ≤ lim inf P (T ) + lim inf N (u).
u→u0
T →T0
Abbiamo perciò dimostrato che il nostro funzionale F (T, u) è debolmente semicontinuo inferiormente in un punto (T0 , u0 ) ∈ H 1 (Ω) × R(Γ1 ). Tuttavia, (T0 , u0 )
rappresenta un punto qualsiasi dell’insieme H 1 (Ω)×R(Γ1 ), quindi F è debolmente
semicontinuo inferiormente in tutto H 1 (Ω) × R(Γ1 ).
5. Sia (Tk , uk ) ∈ H 1 (Ω) × R(Γ1 ) una successione minimizzante tale che |P (Tk )| ≤ M
e G(Tk , uk ) = 0, cioè:
R
| Γ3 (Tk − Td )2 dΓ| ≤ M
−∇ · (k∇Tk ) = 0.
Allora possiamo dire che:
R
− Ω ∇ · (k∇Tk )vdΩ = 0
∀v ∈ H 1 (Ω),
cioè che:
Z
−
Z
∇ · (k∇Tk )vdΩ =
Ω
Z
k∇Tk ∇vdΩ +
(h + uk )Tk vdΓ +
Γ1
Z
Z
−
(h + uk )T∞ vdΓ −
qvdΓ = 0 .
Ω
Γ1
Γ3
v è una qualsiasi funzione appartenente allo spazio H 1 (Ω), quindi possiamo scegliere
v = Tk e sostituirla nell’ultima equazione:
Z
Z
k∇Tk ∇Tk dΩ +
(h + uk )Tk Tk dΓ +
Ω
Γ1
Z
Z
−
(h + uk )T∞ Tk dΓ −
qTk dΓ = 0 ,
Γ1
Γ3
che, portando a destra dell’uguale gli ultimi tre termini, risulta:
Z
k(∇Tk )2 dΩ
Ω
Z
Z
Z
2
=−
(h + uk )(Tk ) dΓ +
(h + uk )T∞ Tk dΓ +
Γ1
Γ1
Analizziamo gli integrali dell’ultima equazione separatamente:
R
• Ω k(∇Tk )2 dΩ = k||∇Tk ||2L2 (Ω) .
19
Γ3
qTk dΓ .
•
+ uk )(Tk )2 dΓ = (h + uk )||Tk ||2L2 (Γ1 ) ,
che, applicando la disuguaglianza di traccia con costante C1 > 0, diventa:
R
2
2
2
2
Γ1 (h + uk )(Tk ) dΓ = (h + uk )||Tk ||L2 (Γ1 ) ≤ C1 (h + uk )||Tk ||H 1 (Ω) .
R
Γ1 (h
• Applicando la disuguaglianza di Cauchy-Schwarz al terzo integrale, otteniamo:
Z
(h + uk )T∞ Tk dΓ ≤ ||(h + uk )T∞ ||L2 (Γ1 ) ||Tk ||L2 (Γ1 )
Γ1
= (h + uk )T∞ |Γ1 | ||Tk ||L2 (Γ1 ) ,
e, applicando la disuguaglianza di traccia con costante C2 > 0:
R
Γ1 (h + uk )T∞ Tk dΓ ≤ C2 (h + uk )T∞ |Γ1 | ||Tk ||H 1 (Ω) .
Poichè C2 è una costante arbitraria ma positiva, possiamo sceglierla in modo
tale da rendere tutto il secondo membro della disequazione maggiore di 1. In
tal caso, possiamo scrivere l’ulteriore disuguaglinza:
R
2
2
2 2
2
Γ1 (h + uk )T∞ Tk dΓ ≤ C2 (h + uk ) T∞ |Γ1 | ||Tk ||H 1 (Ω) .
• Sappiamo per ipotesi che:
R
| Γ3 (Tk − Td )2 dΓ| ≤ M
quindi, la funzione integranda (Tk − Td )2 è limitata e, in particolare, abbiamo che ||Tk − Td ||2L2 (Γ3 ) è limitata, quindi anche la norma di ||Tk ||2L2 (Γ3 ) è
limitata da una costante che chiamiamo M 0 . Considerando l’ultimo integrale
e applicando la disuguaglianza di Cauchy-Schwarz, possiamo dire che:
R
2
2
2
0
Γ3 qTk dΓ ≤ ||q||L2 (Γ3 ) ||Tk ||L2 (Γ3 ) = q |Γ3 | ||Tk ||L2 (Γ3 ) ≤ q |Γ3 | M .
Raggruppando tutti i termini trovati, otteniamo:
2 |Γ |2 ||T ||2
0
k||∇Tk ||2L2 (Ω) ≤ −C12 (h+uk )||Tk ||2H 1 (Ω) +C22 (h+uk )2 T∞
1
k H 1 (Ω) +q |Γ3 | M ,
che, scomponendo la norma H 1 (Ω), diventa:
k||∇Tk ||2L2 (Ω) + C12 (h + uk )||Tk ||2L2 (Ω) + C12 (h + uk )||∇Tk ||2L2 (Ω) +
2
2
−C22 (h + uk )2 T∞
|Γ1 |2 ||Tk ||2L2 (Ω) − C22 (h + uk )2 T∞
|Γ1 |2 ||∇Tk ||2L2 (Ω) ≤ q |Γ3 | M 0 .
Raccogliendo le costanti al primo membro per tornare alla norma in H 1 (Ω):
2
max{k + C12 (h + uk ) − C22 (h + uk )2 T∞
|Γ1 |2 , C12 (h + uk )+
2
− C22 (h + uk )2 T∞
|Γ1 |2 }||Tk ||2H 1 (Ω) ≤
≤ q |Γ3 | M 0 ,
e chiamando:
2 |Γ |2 , C 2 (h+u )−C 2 (h+u )2 T 2 |Γ |2 },
C̃ = max{k +C 2 (h+uk )−C 2 (h+uk )2 T∞
1
1
k
k
∞
otteniamo:
C̃||Tk ||2H 1 (Ω) ≤ q |Γ3 | M 0 ,
20
cioè:
0
||Tk ||2H 1 (Ω) ≤ q |Γ3 | M
.
C̃
0
, arriviamo alla dimostrazione dell’ultima ipotesi:
Ribattezzando M0 = q |Γ3 | M
C̃
||Tk ||2H 1 (Ω) ≤ M0 .
Nel problema di controllo da studiare, esiste quindi una soluzione (T, u) ∈ H 1 (Ω) ×
R(Γ1 ). In generale non ci si aspetta unicità del controllo ottimo.
3.2
Risoluzione numerica: metodo iterativo
La risoluzione numerica di un problema di controllo ottimo richiede un approccio di
tipo iterativo (si vedano [Q08], [D08]). Una volta definito l’algoritmo iterativo, sarà poi
necessario operare una scelta sul tipo di discretizzazione numerica da utilizzare per la
risoluzione delle edp che caratterizzano il problema di controllo in esame.
In corrispondenza dell’ottimo devono essere soddisfatte le tre equazioni (23), (24) e (26),
che rappresentano rispettivamente l’equazione di stato, l’equazione aggiunta e il principio
di minimo, le quali garantiscono che la soluzione trovata sia quella che minimizza il
funzionale lagrangiano espresso dalla (21). Una volta definito quindi il problema di
controllo, il metodo iterativo per la ricerca dell’ottimo può essere schematizzato nel
seguente algoritmo:
• si sceglie un valore iniziale u0 per la variabile di controllo u;
• dato u si risolve l’equazione di stato ricavando cosı̀ il valore della variabile di stato
T;
• noto T si calcola il valore del funzionale costo J (o del funzionale lagrangiano
L: è possibile utilizzare indifferentemente J o L, dato che dalla definizione di
quest’ultimo abbiamo L = J a patto di risolvere esattamente il problema di stato);
• sempre noto T , si risolve l’equazione aggiunta e si ricava il valore della variabile
aggiunta p;
• dato p, si calcola il valore della derivata del funzionale lagrangiano Lu (da qui in
poi indicata con J 0 ) tramite il principio di minimo;
• a questo punto, se il valore del funzionale J e della sua derivata J 0 (che indicano
la distanza dall’ottimo, rispetto ad una tolleranza scelta) rispettano opportuni
criteri di arresto, il ciclo termina e il valore di u all’iterazione corrente coincide
con l’ottimo; in caso contrario si aggiorna opportunamente il valore del controllo
u e si esegue un’altra iterazione, fino a quando non saranno soddisfatti i criteri di
arresto del ciclo.
In particolare abbiamo scelto un criterio di arresto (che indica il raggiungimento dell’ottimo) basato sulla valutazione della norma della derivata del funzionale costo rispetto
21
ad una tolleranza fissata. Il valore di tale tolleranza viene scelto sufficientemente piccolo
in relazione al valore iniziale della quantità monitorata e al livello di vicinanza all’ottimo
reale che si vuole ottenere. Il criterio di arresto per ogni iterazione k si traduce quindi
nella seguente:
0 k J (u )
≤ 10−6 .
(33)
kJ 0 (u0 )k
Per quanto riguarda l’aggiornamento della variabile di controllo u è possibile utilizzare
un metodo di tipo gradiente, cioè un algoritmo iterativo che data la soluzione al passo
k calcola quella al passo successivo k + 1 nel modo seguente (si utilizza la notazione
vettoriale in quanto più generale):
uk+1 = uk + τ k dk
k = 0, 1, 2, ... ,
(34)
dove τ k indica un parametro reale e strettamente positivo detto parametro di accelerazione (o di rilassamento), mentre dk rappresenta la direzione di discesa, cioè tale per
cui vale:
dk T · ∇J(uk ) < 0 se ∇J(uk ) 6= 0.
A seconda delle possibili scelte per la direzione dk , si ottengono diversi metodi iterativi
di tipo gradiente, quali ad esempio il metodo di Newton, il metodo di quasi-Newton,
il metodo di discesa più ripida e il metodo del gradiente coniugato. In particolare per
il problema di controllo in esame abbiamo scelto di utilizzare gli ultimi due metodi sopraelencati.
Il metodo di discesa più ripida (SD, steepest descent) corrisponde alla scelta dk =
−∇J(uk ) e quindi per il caso in esame (dove u ∈ R) l’algoritmo iterativo risulta essere:
uk+1 = uk − τ k J 0 (uk )
k = 0, 1, 2, ... .
(35)
La valutazione della derivata del funzionale J 0 in un dato punto del dominio di controllo fornisce un’indicazione della sensitività del funzionale J in quel punto rispetto alle
variazioni della variabile di controllo u: questo fa sı̀ che il metodo di discesa più ripida
sia, tra gli altri possibili, il metodo più utile alla comprensione del ruolo giocato da J 0
e quindi dalla variabile aggiunta p.
Il metodo del gradiente coniugato (CG, conjugate gradient) invece corrisponde alla
seguente scelta per la direzione di discesa dk :
(
−∇J(uk )
per k = 0,
dk =
k
k
k−1
−∇J(u ) + γ d
per k ≥ 1,
dove γ k è uno scalare da scegliere in modo tale che dk T · dk−1 = 0. Questo parametro
può essere calcolato in uno dei possibili modi seguenti:
∇J(uk )T ∇J(uk ) − ∇J(uk−1 )
k
γP RP =
(Polak-Ribiere-Poylak),
∇J(uk−1 )T ∇J(uk−1 )
n
o
γPk RP + = max γPk RP , 0
(Polak-Ribiere-Poylak modificato),
k
γCD
=−
∇J(uk )T ∇J(uk )
dk−1 T ∇J(uk−1 )
(conjugate descent).
22
Dunque per il problema di controllo in esame il metodo del gradiente coniugato si traduce
in:
uk+1 = uk + τ k dk
k = 0, 1, 2, ... ,
(
−J 0 (uk )
per k = 0,
con dk =
0
k
k
k−1
−J (u ) + γ d
per k ≥ 1,
(36)
e γ k può essere calcolato in uno dei seguenti modi:
γPk RP
=
J 0 (uk )
J 0 (uk ) − J 0 (uk−1 )
kJ 0 (uk−1 )k
n
o
γPk RP + = max γPk RP , 0 ,
0 k 2
J (u )
k
γCD
= − k−1 0 k−1 .
d
J (u )
2
,
(37)
(38)
(39)
Entrambi i metodi presentati comportano la valutazione di un appropriato valore del
parametro di accelerazione τ k . La scelta di τ k deve essere in grado di garantire la
convergenza del funzionale costo al minimo, deve dunque valere ad ogni passo del ciclo
iterativo la seguente condizione di discesa:
J(uk+1 ) < J(uk ) .
(40)
Vi sono diversi metodi approssimati basati su criteri euristici per il calcolo di τ k ;
presentiamo qui alcune possibili scelte:
1. si pone τ k = τ̄ ∀k = 0, 1, 2, ..., dove il valore fissato τ̄ viene scelto in modo da
soddisfare ad ogni iterazione la condizione di discesa (40); per alcuni semplici
problemi di ottimizzazione è infatti possibile determinare un valore ottimale τopt
a priori dall’algoritmo iterativo sulla base delle proprietà del problema stesso;
2. fissato un valore costante τ̄ di partenza, possibilmente grande, la successione τ k
è scelta in modo che τ k = τ̄ /2j dove j = minq∈N J(uk + τ̄ /2q dk ) < J(uk ); in
altre parole ad ogni iterazione τ k viene posto pari al valore fissato τ̄ , si calcolano
il controllo uk+1 e il corrispondente valore del funzionale J(uk+1 ) e si verifica la
condizione di discesa: se la (40) è soddisfatta si procede con il ciclo iterativo, in
caso contrario, ovvero se uk+1 fornisce una soluzione peggiore della precedente, si
ritorna al valore uk e si ricalcola uk+1 con un passo ridotto, ponendo cioè τ k = τ̄ /2,
e cosı̀ via finchè la condizione di discesa non risulta verificata; questo metodo è
una versione particolare e semplificata del metodo di Armijo, che rientra nella
categoria dei metodi line search;
3. nel caso in cui risulti noto il valore finito del funzionale costo in corrispondenza
dell’ottimo, Jopt = inf u∈U J(u) ≥ 0, si può scegliere il parametro di accelerazione
2
come τ k = J(uk ) − Jopt / ∇J(uk ) ; in questo caso però non è più garantita
23
la monotonia della successione J(uk ) , il che significa che può accadere che la
condizione (40) non sia rispettata ad ogni iterazione; per questo motivo il valore
attribuito a Jopt deve essere il più preciso possibile, in quanto valori anche poco
distanti da quello esatto possono rendere molto lenta la convergenza dell’algoritmo
iterativo e in alcuni casi anche portare a divergenza .
La scelta del parametro τ k , oltre a soddisfare la (40), la quale assicura la convergenza
del funzionale al minimo, deve anche garantire che la convergenza venga raggiunta il più
rapidamente possibile. È possibile quindi introdurre criteri più severi della (40) nella
scelta dei possibili valori di τ k , con il fine di evitare due diversi tipi di problemi: una
velocità di decrescita troppo bassa e l’uso di passi troppo piccoli. La prima difficoltà si
può ovviare richiedendo che:
J(uk ) − J(uk+1 ) ≥ −στ k dk T · ∇J(uk ) ,
(41)
dove σ ∈ (0, 1/2). Questa condizione equivale a richiedere che la velocità media di
decrescita lungo dk di J in uk+1 sia almeno pari ad una frazione assegnata della velocità
di decrescita iniziale in uk . La seconda difficoltà invece può essere superata richiedendo
che la velocità di decrescita in uk+1 non sia inferiore ad una frazione assegnata di quella
in uk , attraverso la seguente:
kT
kT
k+1 k (42)
d · ∇J(u ) ≤ η d · ∇J(u ) ,
dove η ∈ (σ, 1), cosı̀ da garantire che valga contemporaneamente anche la (41). In
particolare per il nostro problema di controllo le condizioni precedenti si scrivono come:
J(uk ) − J(uk+1 ) ≥ −στ k dk J 0 (uk )
k 0 k+1 d J (u ) ≤ η dk J 0 (uk ) .
Nella pratica si scelgono σ ∈ [10−5 , 10−1 ] e η ∈ [0.1, 0.5]. Dunque per il calcolo del
parametro τ k si sceglierà uno dei metodi sopraelencati compatibilmente con il rispetto
delle condizioni (40), (41) e (42).
24
4
Implementazione del codice
In questa sezione diamo una descrizione generale del codice da noi implementato, senza
soffermarci troppo sui dettagli, per i quali rimandiamo alla documentazione Doxygen.
Sviluppiamo il codice per le simulazioni numeriche tramite il linguaggio di programmazione C++. La risoluzione del problema in esame richiede l’implementazione delle
seguenti operazioni:
1. lettura dei parametri fisici e numerici legati al problema;
2. creazione della mesh computazionale;
3. creazione e risoluzione dei sistemi lineari associati all’equazione di stato e all’equazione aggiunta;
4. ciclo iterativo legato alla ricerca del parametro di controllo ottimale;
5. post-processing dei risultati.
Ognuna delle precedenti operazioni è costituita da una o più unità di programmazione,
ciascuna suddivisa in due file: il primo, con estensione .hpp, contiene gli header file da
includere (librerie o altri file .hpp) e le dichiarazioni delle classi e/o delle funzioni, che
vengono poi definite nel secondo file corrispondente con estensione .cpp. Questa suddivisione non vale ovviamente per i file che utilizzano al loro interno strutture template:
in questo caso l’intero codice è contenuto nel file .hpp.
Verranno utilizzate alcune delle strutture dati della Standard Template Library, tra
cui vector, list, set e map. Inoltre per la risoluzione dei sistemi lineari si è scelto di
utilizzare la libreria UMFPACK, che implementa un metodo di risoluzione diretto molto
efficace nel caso di sistemi lineari con matrici sparse.
Spieghiamo adesso brevemente il programma, passo per passo.
4.1
Pre-processing e mesh computazionale
Per la lettura dei parametri fisici e numerici che caratterizzano il problema, si è utilizzato
GetPot: un parser che permette di processare un file di testo (chiamato datafile) che
contiene tutte le informazioni utili per l’esecuzione del programma. L’enorme vantaggio
fornito dall’impiego di un file di questo tipo per il passaggio dei parametri consiste nella
possibilità di operare modifiche su di essi (e quindi anche di risolvere problemi differenti)
senza dover necessariamente ricompilare il codice. All’interno del main si instanzia un
oggetto di classe GetPot, che verrà poi passato ai costruttori delle classi che necessitano
di leggere parametri dall’esterno. Scegliamo di creare la mesh computazionale tramite
FreeFem++, salvandola in un file di estensione .msh. Il nostro programma leggerà poi
tale file (tramite la libreria fstream), memorizzando la mesh in una opportuna struttura
dati. Tale struttura è costituita da quattro classi:
class Point;
25
class Triangle;
class Edge;
class MeshFF.
Un oggetto di classe Point contiene tutte le informazioni relative ai nodi che costituiscono la griglia.
class Point{
private:
double _x;
double _y;
int _id_node;
int _id_boundary;
...
}
Ogni nodo viene memorizzato tramite: le coordinate x e y, il numero identificativo
id node ad esso associato e il numero id boundary che indica il bordo a cui appartiene
(pari a 0 se si tratta di un nodo interno). I metodi di questa classe permettono di creare
un determinato nodo e di ricavare o stampare a video i parametri relativi a tale nodo.
Ogni elemento (triangolare) della griglia viene definito tramite i suoi tre vertici,
ognuno dei quali è un oggetto di classe Point. L’oggetto di classe Triangle definisce
quindi il singolo elemento della mesh e contiene: un vettore di tre puntatori ad oggetti
di classe Punto, element (si è preferito creare un vettore di puntatori e non di oggetti
di classe Punto perchè, dal momento che nella struttura dati della mesh ci sarà già un
vettore contenente tutti i nodi della griglia, sarebbe inutile memorizzarli una seconda
volta) ; il numero identificativo di ogni elemento, id element. I metodi di questa
classe permettono di creare un certo triangolo e di ricavare o stampare a video tutte le
informazioni relative all’elemento stesso o ai nodi ad esso associati.
class Triangle{
private:
int _id_element;
std::vector <Point*> _element;
std::vector <double> _cof_jacob;
...
}
Dal momento che il calcolo degli integrali bidimensionali, necessari per la costruzione
dei sistemi lineari associati al problema di stato e al problema aggiunto, verrà effettuato
passando all’elemento di riferimento (il triangolo di vertici (0, 0), (0, 1) e (1, 0)), è necessario conoscere i coefficienti della matrice Jacobiana associata alla trasformazione lineare
affine tra elemento generico della griglia ed elemento di riferimento. Tali coefficienti si
possono ricavare direttamente dalle coordinate dei tre vertici del triangolo corrente e
26
sono memorizzati nel vettore cof jacob; la classe contiene inoltre un metodo per il calcolo dell’inverso del determinante della matrice Jacobiana associata ad ogni elemento,
anch’esso necessario per il calcolo degli integrali bidimensionali.
Il calcolo del funzionale costo e l’aggiunta delle condizioni al contorno all’interno dei
sistemi lineari richiedono il calcolo di integrali monodimensionali sui lati di bordo. Per
questo motivo è necessario creare una classe che contenga tutte le informazioni relative
ai lati di bordo.
class Edge{
private:
int _id_sx;
int _id_dx;
int _id_bound;
double _length;
std::vector <Point*> _vertices;
...
}
L’oggetto di classe Edge contiene, per ogni lato di bordo: il numero identificativo del
bordo a cui appartiene ( id bound), i numeri identificativi dei nodi che rappresentano
l’estremo sinistro ( id sx) e l’estremo destro ( id dx), tutte le informazioni sui due
estremi all’interno di vertices, che è un vettore di puntatori ad oggetti di classe Punto.
Dal momento che anche gli integrali monodimensionali vengono calcolati sull’elemento
di riferimento (il segmento di lunghezza unitaria), è necessario memorizzare la lunghezza
di ogni lato di bordo ( length).
Le tre classi appena descritte, Point, Triangle e Edge, vengono utilizzate dalla
classe MeshFF per memorizzare la griglia computazionale, leggendo tutte le informazioni
necessarie dal file di estensione .msh.
class MeshFF{
private:
int _np;
int _ne;
int _nb;
std::vector <Point> _nodes;
std::vector <Triangle> _grid;
std::vector <Edge> _edges;
std::list <int>* _near;
double _length1;
double _length2;
double _length3;
...
}
27
Tutti gli elementi della griglia sono memorizzati nel vettore grid, che contiene elementi
di classe Triangle; i nodi della griglia sono salvati invece nel vettore nodes, che contiene
elementi di classe Point; infine tutti i lati di bordo sono memorizzati nel vettore edges,
che contiene elementi di classe Edge. I parametri np, ne e nb indicano, rispettivamente,
il numero totale di nodi, di triangoli e di lati di bordo della mesh, necessari per allocare lo
spazio per i tre vettori corrispondenti nodes, grid e edges. Inoltre la classe contiene
le lunghezze dei bordi Γ1 , Γ2 e Γ3 nei tre parametri length1, length2 e length3.
Infine, per eseguire il calcolo di tutti gli integrali, sia quelli monodimensionali che quelli
bidimensionali, è necessario conoscere, per ogni nodo della griglia, tutti i nodi a lui vicini,
cioè collegati ad esso tramite un lato. L’elenco di tutti i nodi vicini viene memorizzato
attraverso near, un puntatore ad una lista di interi, che verrà istanziato in modo tale
che per ogni nodo contenga la lista di tutti i suoi vicini, se stesso compreso. I vicini sono
memorizzati in ordine crescente. I metodi della classe Mesh permettono di costruire la
struttura dati contenente la griglia; consentono inoltre di trovare e salvare la lista dei
vicini per ogni nodo e di stampare a video tutte le altre informazioni relative alla mesh
computazionale.
Tutta la struttura dati per la mesh appena illustrata è pensata per essere costruita
a partire da un file di tipo .msh fornito da FreeFem++, che noi abbiamo scelto come
generatore di griglia per il nostro problema. Il programma è comunque in grado di
supportare l’utilizzo di una generica classe Mesh, diversa da MeshFF, in quanto tutto il
codice è template su Mesh: è dunque possibile costruire e utilizzare una struttura dati
diversa da quella presentata, che possa ad esempio leggere le informazioni a partire da
un diverso generatore di mesh, mantenendo il resto del codice invariato.
4.2
Integrazione numerica
Come già accennato, la risoluzione del problema di controllo richiede l’implementazione
di due integratori, uno bidimensionale e uno monodimensionale. In entrambi i casi, il
calcolo viene riportato sull’elemento di riferimento.
Gli integrali bidimensionali da calcolare sono associati ai termini diffusivi dell’equazione di stato e dell’equazione aggiunta e assumono entrambi la seguente forma:
Z
k
∇ω · ∇vdΩ ,
Ω
dove ω indica l’incognita T o la variabile aggiunta p e v indica la funzione test associata.
Si devono perciò valutare espressioni del tipo:
Z
k
∇φi (x) · ∇φj (x)dΩ ,
(43)
Ω
dove φi (x) e φj (x) indicano le funzioni di base associate al nodo i e al nodo j definite
localmente su ogni singolo elemento. Per riportare il calcolo di ogni singolo integrale
sull’elemento di riferimento, è necessario eseguire il seguente cambio di variabili:
x = J x̂ + c ,
28
dove x̂ indica il sistema di coordinate cartesiane sull’elemento di riferimento (il triangolo
di vertici (0, 0), (0, 1) e (1, 0)), x indica il sistema di coordinate sull’elemento generico
corrente e la matrice J è lo Jacobiano della trasformazione affine. In particolare, la
matrice J e il vettore c sono definiti nel seguente modo:
x2 − x1 x3 − x1
x1
J=
, c=
.
y 2 − y 1 y 3 − y1
y1
Sull’elemento di riferimento le funzioni di base degli elementi finiti lineari sono definite,
per ogni vertice, nel seguente modo:
φ1 (x̂) = 1 − x̂ − ŷ ,
φ2 (x̂) = ŷ ,
φ3 (x̂) = x̂ ,
dove φ1 è associata al nodo (0, 0), φ2 è associata al nodo (0, 1) e φ3 è associata al nodo
(1, 0). Il calcolo dell’integrale bidimensionale richiede però la valutazione dei gradienti
delle funzioni di base che, riportati sull’elemento di riferimento, assumono per il vertice
i la seguente forma:
∇x φi (x) = ∇x φi (J x̂ + c) = (J −1 )T ∇x̂ φi (x̂) ,
cioè, per ogni singolo nodo:
∇x̂ φ1 (x̂) = (−1, −1)T ,
∇x̂ φ2 (x̂) = (0, 1)T ,
∇x̂ φ3 (x̂) = (1, 0)T .
Applicando il cambio di variabili da x a x̂, cioè passando dall’elemento corrente all’elemento di riferimento, l’integrale da calcolare (43) diventa:
Z
k
∇Tx̂ φi (x̂)J −1 (J −1 )T ∇x̂ φi (x̂) det(J) dx̂ ,
(44)
Ω̂
dove Ω̂ indica l’insieme degli elementi correnti riportati sugli elementi di riferimento.
Tale formula può essere riscritta esplicitando i prodotti matrice-vettore.
Per implementare l’integratore bidimensionale è quindi necessario definire le funzioni
di base sull’elemento di riferimento e creare una classe che implementi una formula
di integrazione numerica di tipo Gaussiano. I file base.hpp e base.cpp contengono
rispettivamente le dichiarazioni e le definizioni dei prodotti scalari tra i gradienti delle
funzioni di base, necessari per il calcolo dell’espressione (44). La classe Quadratura
contiene invece la formula di quadratura Gaussiana a tre nodi:
29
class Quadratura {
private:
std::vector<double> nodes_x;
std::vector<double> nodes_y;
double weight;
pfn integrand;
public:
Quadratura();
double integrate(pfn f);
}
I vettori nodes_x e nodes_y contengono le coordinate dei nodi di quadratura in direzione
x e y rispettivamente; weight contiene il valore del peso di quadratura; pfn indica un
puntatore ad una funzione di due variabili ed è definito nel file typedef.hpp: integrand
contiene la funzione da integrare, cioè una opportuna funzione definita nel file base.cpp.
Il metodo integrate implementa la formula di quadratura gaussiana sommando il valore
della funzione integrand nei nodi di quadratura e moltiplicando il valore finale per il
peso di quadratura.
Gli integrali monodimensionali da calcolare servono sia per la valutazione del funzionale costo e della sua derivata, sia per includere le condizioni al contorno dei sistemi
associati all’equazione di stato e all’equazione aggiunta. Le condizioni al bordo vanno in
parte aggiunte alla matrice del sistema e in parte al termine noto. Tutti questi termini
possono essere ridotti alla valutazione di due diversi integrali:
Z
ψi (x)ψj (x)dΓ ,
ZΓ
ψi (x)dΓ ,
Γ
dove le funzioni di base ψi sono la restrizione delle funzioni di base φi sui bordi del
dominio di integrazione. Per passare all’elemento di riferimento (segmento di lughezza
unitaria: estremo sinistro x̂1 = 0 ed estremo destro x̂2 = 1) è necessario il cambio di
variabili:
x = (x2 − x1 )x̂ + a ,
dove x1 e x2 sono riferiti all’elemento corrente, x indica la coordinata dell’elemento
corrente e x̂ indica la coordinata dell’elemento di riferimento. Le funzioni di base ψ
sull’elemento di riferimento sono:
ψ1 (x̂) = 1 − x̂ ,
ψ2 (x̂) = x̂ ,
30
dove ψ1 è la funzione di base associata all’estremo sinistro e ψ2 è la funzione di base
associata all’estremo destro. Il calcolo degli integrali monodimensionali richiede quindi
la valutazione dei seguenti termini:
Z
(x2 − x1 )
ψ1 (x̂)ψ2 (x̂)dΓ̂ ,
ZΓ̂
(x2 − x1 )
ψ1 (x̂)ψ1 (x̂)dΓ̂ ,
ZΓ̂
(x2 − x1 )
ψ2 (x̂)ψ2 (x̂)dΓ̂ ,
ZΓ̂
(x2 − x1 )
ψ1 (x̂)dΓ̂ ,
ZΓ̂
(x2 − x1 )
ψ2 (x̂)dΓ̂ ,
Γ̂
dove Γ̂ indica l’elemento di riferimento monodimensionale. Grazie alla simmetria delle
funzioni di base e a quella dell’elemento di riferimento, il valore del secondo integrale
coincide con quello del terzo, mentre il valore del quarto coincide con quello dell’ultimo.
Pertanto, basta calcolare tre soli integrali. Le funzioni di base necessarie (diversa,
uguale e bordo) sono definite nel file base.cpp e la formula di quadratura è definita
nella classe Quadratura1d:
class Quadratura1d {
private:
std::vector<double> nodes_x;
std::vector<double> weight;
pfn1d integrand;
public:
Quadratura1d();
double integrate(pfn1d f);
}
Il vettore nodes_x indica la posizione dei nodi di quadratura; il vettore weight contiene
il valore dei pesi di quadratura per ogni nodo; pfn1d indica un puntatore a una funzione
di una variabile ed è definito nel file typedef.hpp: integrand contiene la funzione da
integrare, cioè una opportuna funzione definita nel file base.cpp. Il metodo integrate
implementa la formula di quadratura sommando il valore della funzione integrand nei
nodi di quadratura e moltiplicando il valore finale per il peso di quadratura: la formula
di quadratura scelta è quella di Cavalieri-Simpson.
4.3
Assemblaggio e risoluzione dei sistemi lineari
L’assemblaggio e la risoluzione dei sistemi lineari associati al problema di stato e al problema aggiunto, il calcolo del funzionale costo e il calcolo della sua derivata, richiedono l’u31
tilizzo di tutte le strutture definite precedentemente. E’ inoltre necessaria una apposita
struttura dati dove poter memorizzare le matrici associate ai sistemi lineari.
I sistemi lineari vengono risolti tramite la libreria UMFPACK, che richiede in ingresso
una matrice nel formato CSC (Compressed Sparse Column). Pertanto, tutte le matrici
andranno salvate in tale formato. La struttura dati per le matrici in formato CSC consiste
in tre puntatori: matrix contiene gli elementi non nulli della matrice per colonna, rows
contiene gli indici di riga degli elementi di matrix e columns contiene nella i-esima
componente la posizione in rows del primo elemento della i-esima colonna della matrice.
Il pattern della matrice viene creato dalla classe Matrice:
class Matrice {
private:
int* columns;
int* rows;
double* matrix;
int ncol;
int nz;
int nrow;
public:
void fill_Matrice(int i, int j, double value);
void modify_Matrice(int i, int j, double value);
...
}
ncol indica il numero totale di colonne della matrice, nrow indica il numero totale di
righe (uguale al numero di colonne dal momento che la matrice è quadrata) e nz indica il
numero totale di elementi non nulli. Una volta creata la matrice, è possibile accedere ai
singoli elementi e attribuirgli un valore diverso da zero tramite il metodo fill_Matrice
oppure modificarli aggiungendo un certo valore a quello già presente tramite il metodo
modify_Matrice.
La classe Problem permette di calcolare tutte le incognite del problema tramite
la risoluzione del sistema di stato e il sistema aggiunto e tramite la valutazione del
funzionale costo. Si tratta di una classe template su Mesh. Si noti inoltre che in questo
caso tutti gli attributi sono dichiarati protected per implementare successivamente
l’ereditarietà tra classi.
template <typename Mesh>
class Problem{
protected:
Mesh& _grid;
GetPot& _datafile;
double _J;
double _dJ;
int _pbdim;
32
Matrice _Astiffness;
double* _sol_stato;
double* _sol_aggiunto;
...
public:
Problem(Mesh& grid, GetPot& datafile);
void stato(const double u);
void aggiunto(const double u);
double* solve(Matrice A, double* b);
void funzionale(const double u);
void der_funzionale(const double u);
...
}
La parte legata al termine di diffusione del sistema di stato è identica a quella del
sistema aggiunto e rimane inalterata ad ogni passo del ciclo iterativo dal momento che
non contiene la variabile di controllo u. La matrice _Astiffness viene perciò creata
una sola volta per entrambi i problemi dal costruttore della classe Problem, chiamando
la funzione diffusione definita nel file assembla.hpp.
template <typename Mesh>
void diffusione(const Mesh& grid, Matrice& A, const double kappa){
...
double inv_jacob(0);
std::vector<double> integral(9,0);
Quadratura regola;
for (int i=0;i<ne;++i){
...
for(int j=0;j<3;++j){
col=grid.get_element(i).get_id_vertex(j);
for(int k=0;k<3;++k){
row=grid.get_element(i).get_id_vertex(k);
A.fill_Matrice(row,col,integral[3*j+k]*inv_jacob);}
}
}
}
La funzione diffusione itera su tutti gli elementi della mesh e per ognuno di essi
calcola i nove integrali associati alle combinazioni dei tre vertici di ogni triangolo: gli
integrali bidimensionali devono essere valutati per tutti i possibili prodotti dei gradienti
delle funzioni di base associate ad ogni nodo, come indicato dalla formula (44). Gli
integrali calcolati per ogni elemento vengono poi inseriti all’interno della matrice tramite
il metodo fill_matrice che riceve in ingresso gli indici di riga e colonna dell’elemento
da riempire e il valore dell’integrale riportato sull’elemento corrente.
33
Tornando alla classe Problem, il metodo stato costruisce e risolve il sistema associato
all’equazione di stato. Tramite la funzione bc definita nel file assembla.hpp si aggiunge
alla matrice A_stiffness il termine incognito che deriva dalle condizioni al contorno:
si tratta di calcolare degli integrali monodimensionali sul bordo Γ1 (tramite un ciclo
sui lati di bordo) e inserirli nella matrice tramite il metodo modify_Matrice. Tramite
la funzione termine_noto_stato, sempre definita nel file assembla.hpp, si costruisce
invece il termine noto b: si tratta di calcolare degli integrali monodimensionali sul bordo
Γ1 e sul bordo Γ3 (tramite un ciclo sui lati di bordo).
template <typename Mesh>
void bc(const Mesh& grid, Matrice& A, const double h, const double u){
double lunghezza(0.0);
Quadratura1d regola1d;
std::vector<double> integral1d(2,0);
integral1d[0]=(h+u)*regola1d.integrate(uguale);
integral1d[1]=(h+u)*regola1d.integrate(diversa);
for (int i=0;i<grid.get_nb();++i){
...
}
}
template <typename Mesh>
double* termine_noto_stato(const Mesh& grid, const double h,
const double u, const double Tinf, const double q){
double lunghezza(0.0);
Quadratura1d regola1d;
std::vector<double> integral1d(2,0);
integral1d[0]=Tinf*(h+u)*regola1d.integrate(bordo);
integral1d[1]=q*regola1d.integrate(bordo);
double* b;
b= new double[grid.get_num_points()];
...
for (int i=0;i<grid.get_nb();++i){
lunghezza=grid.get_edges(i).get_length();
if (grid.get_edges(i).get_id_bound()==1){
...
}
if (grid.get_edges(i).get_id_bound()==3){
...
}
}
}
Il metodo aggiunto costruisce e risolve il sistema associato all’equazione aggiunta in
modo totalmente analogo al metodo stato. Tramite la funzione bc e la funzione
34
termine_noto_aggiunto, entrambe definite nel file assembla.hpp, il termine incognito derivante dalle condizioni al contorno viene aggiunto alla matrice A_stiffness e
viene costruito il termine noto associato al sistema aggiunto.
Sia il metodo stato che il metodo aggiunto utilizzano il metodo solve, che permette
di risolvere i sistemi lineri associati al problema di stato e al problema aggiunto tramite
l’utilizzo della libreria UMFPACK e calcolare quindi le incognite T e p rispettivamente.
template <typename Mesh>
double* Problem::solve(Matrice _A, double* _F){
double* _solution;
_solution=new double[_pbdim];
for(int i=0;i<_pbdim;++i){
_solution[i]=0;
}
double*null = (double *) NULL ;
void *Symbolic, *Numeric;
(void) umfpack_di_symbolic (_pbdim, _pbdim, _A.get_columns(),
_A.get_rows(), _A.get_matrix(), &Symbolic, null, null);
(void) umfpack_di_numeric (_A.get_columns(), _A.get_rows(),
_A.get_matrix(), Symbolic, &Numeric, null, null);
umfpack_di_free_symbolic (&Symbolic);
(void) umfpack_di_solve (UMFPACK_A, _A.get_columns(),
_A.get_rows(), _A.get_matrix(), _solution, _F, Numeric, null,null);
umfpack_di_free_numeric (&Numeric) ;
return _solution;
...
}
Il funzionale costo J viene calcolato tramite il metodo funzionale, che utilizza l’integratore monodimensionale e la soluzione T del problema di stato _sol_stato. La parte
contenente la variabile di controllo u non ha bisogno di integrazione dal momento che
tale variabile è reale e non appartiene quindi allo spazio degli elementi finiti.
template <typename Mesh>
void Problem::funzionale(const double u){
_J=0.0;
Quadratura1d regola1d;
std::vector<double> integral1d(3,0);
integral1d[0]=regola1d.integrate(uguale);
integral1d[1]=regola1d.integrate(bordo);
integral1d[2]=regola1d.integrate(diversa);
...
}
Infine, la derivata del funzionale costo dJ viene calcolata tramite il metodo der_funzionale,
che funziona in modo totalmente analogo al metodo funzionale.
35
4.4
Ciclo iterativo
Come illustrato nella sezione 3.2, per risolvere il problema di controllo ottimo è necessario
effettuare un ciclo iterativo. Tale ciclo è contenuto nel main del nostro programma.
int main(int argc, char** argv)
{
...
Ciclo<MeshFF> controllo(mesh,data_file);
controllo.inizializza();
while ((controllo.get_dJtol()>controllo.get_toll()) &&
(controllo.get_iter()<=controllo.get_nmax())){
controllo.itera();
...
}
}
La struttura che gestisce i passi iterativi del ciclo while presente nel main è costituita
dalla classe Ciclo. Si tratta anche in questo caso di una classe template su Mesh.
Essa eredita pubblicamente dalla classe Problem e può quindi utilizzare tutti i metodi
e tutti gli attributi di tale classe. In aggiunta a questi, la classe Ciclo contiene anche
tutti i parametri necessari per l’implementazione dei metodi numerici di ottimizzazione:
discesa più ripida e gradiente coniugato.
template <typename Mesh>
class Ciclo:
public Problem<Mesh>{
private:
flag _flag_metodo;
flag _flag_tau;
flag _flag_gamma;
flag _flag_cond;
public:
Ciclo(Mesh& grid, GetPot& datafile);
void inizializza();
void itera();
void itera_SD();
void itera_CG();
...
}
La classe contiene in particolare quattro flag, che vengono lette nel costruttore dal
datafile, e che indicano il metodo iterativo che si intende utilizzare: _flag_metodo
permette di scegliere tra il metodo di discesa più ripida (35) e metodo del gradiente
36
coniugato (35); _flag_tau consente la scelta tra due diversi metodi per il calcolo del
parametro di accelerazione τ (si vedano le scelte 2 e 3 nella sezione 3.2); _flag_gamma
permette di scegliere tra i tre possibili modi di calcolare il parametro γ del metodo del
gradiente coniugato ((37), (38), (39)); _flag_cond consente di scegliere quale delle tre
condizioni (40), (41), (42) usare nel ciclo di aggiornamento di τ , nel caso si sia scelta
l’opzione 2 per il calcolo di questo parametro. Il metodo inizializza a partire da
un valore iniziale del controllo u0 , calcola il valore iniziale della temperatura T , della
variabile aggiunta p e il valore iniziale del funzionale costo e della sua derivata. Il metodo itera compie l’iterazione: tramite _flag_metodo si sceglie se chiamare il metodo
itera_SD (che implementa il metodo steepest descent) o itera_CG (che implementa
il metodo conjugate gradient). Il ciclo si arresta quando risulta essere soddisfatta la
condizione (33) sulla derivata del funzionale costo, che indica la vicinanza all’ottimo,
oppure quando viene raggiunto il numero massimo di iterazioni consentite.
4.5
Post-processing
Il post-processing dei risultati è affidato al metodo postprocess della classe Ciclo, che
permette di salvare i risultati ottenuti in opportuni file per la visualizzazione.
Tramite la libreria fstream, la distribuzione finale della temperatura T e della variabile
aggiunta p vengono scritti su due file; viene inoltre creato il file postproc.edp, che
contiene le istruzioni necessarie per leggere i due file contenenti le variabili T e p, e che
permette la lora visualizzazione tramite FreeFem++; il metodo inoltre lancia direttamente
l’esecuzione di FreeFem++ tramite il comando system.
Il metodo postprocess consente inoltre di scrivere su file .dat i valori di alcune variabili
del problema, salvate in opportuni vettori durante l’esecuzione del ciclo iterativo. In
particolare vengono salvati in questo modo i valori del funzionale costo J, della sua
derivata J 0 e della variabile di controllo u. Vengono salvati inoltre su un altro file .dat
i valori della temperatura T corrispondenti ai nodi sul bordo Γ3 , per poi visualizzarne
l’andamento. Tali file verranno successivamente letti tramite il programma GnuPlot,
che renderà possibile la visualizzazione dei grafici di interesse.
37
5
Risultati numerici
In questa sezione riportiamo alcuni dei risultati numerici ottenuti con le simulazioni
effettuate con il codice C++ illustrato nella sezione 4.
Consideriamo inizialmente il metodo di discesa più ripida, con un valore iniziale della
variabile di controllo u0 = 0. Scegliamo di fissare il valore di partenza del parametro τ
pari a 500 e di aggiornare di volta in volta il parametro τ dimezzandolo (opzione 2 per τ
nella sezione 3.2), in base al sottociclo iterativo che ha come condizione di uscita J > Jold ,
dove J e Jold indicano il valore del funzionale costo rispettivamente all’iterazione corrente
e a quella precedente. Il valore della tolleranza che permette l’uscita dal ciclo while
principale è fissato pari a 10−6 , mentre il numero massimo di iterazioni è posto pari
a 100. Decidiamo inoltre di fissare il valore del parametro di coercività β (si veda la
definizione del funzionale costo (15)) pari a 10−3 . Il valore ottimale del parametro di
controllo u e del funzionale costo J, raggiunti dopo n iterazioni, sono i seguenti:
numero totale di iterazioni = 11
valore ottimale di u = 38.2395
valore finale del funzionale J = 0.773907
Il miglior funzionamento del dissipatore di calore per il chip è quindi garantito con un
coefficiente di convezione forzata dell’aria (h + u) pari a 43.94 mW
2◦ C , ricordando che il
coefficiente di convezione naturale h nel nostro caso è pari a 5.7 mW
2◦ C .
Il valore ottimo trovato per u è quello che minimizza il valore del funzionale costo J(u):
plottando infatti il valore del funzionale J al variare delle iterazioni, si può vedere come
il valore di J decresca ad ogni iterazione e che il minimo è raggiunto in corrispondenza
dell’ultima iterazione, associata appunto al valore ottimale di u (figura 2). La derivata
del funzionale costo J 0 assume un valore dell’ordine di −102 alla prima iterazione, corrispondente al caso con controllo nullo, ma dalla seconda iterazione in poi si mantiene
molto piccolo fino a convergere a 0 (figura 3).
La distribuzione della temperatura sul dissipatore (figura 4 in alto) presenta valori più
elevati in corrispondenza della fonte di calore (il chip) che si trova a contatto con la superficie in basso e valori che diminuiscono gradualmente allontanandosi da esso, grazie
al lavoro di dispersione termica effettuato passivamente dal dissipatore e alla presenza
della ventola. La distribuzione della variabile aggiunta ha lo stesso andamento della temperatura, ma i suoi valori sono molto vicini allo zero (figura 4 in basso). La temperatura
desiderata di funzionamento del chip è pari a 65◦ C: i valori di temperatura sul bordo
Γ3 , ovvero la superficie inferiore del dissipatore, che si trova a contatto con il chip, sono
di poco superiori a tale valore (figura 5). Il dissipatore e la ventola permettono quindi
di mantenere la temperatura sul chip intorno al valore ideale desiderato, garantendo il
corretto funzionamento del chip stesso. L’andamento della temperatura sul bordo Γ3
è influenzato in modo evidente dalla geometria del dissipatore: si nota infatti che la
temperatura risulta più bassa nella zona centrale del bordo, al di sopra della quale si
trova il tronco centrale del dissipatore dove avviene la convezione, mentre rimane più
elevata nelle regioni laterali, più vicine alla parte di bordo adiabatica e dunque meno
influenzata dal raffredamento della zona superiore.
38
Figura 2: Metodo SD: funzionale costo J in funzione delle iterazioni.
Figura 3: Metodo SD: derivata del funzionale costo J 0 in funzione delle iterazioni.
39
Figura 4: Distribuzione della temperatura T e della variabile aggiunta p sul dissipatore in
corrispondenza dell’ottimo (metodo SD).
40
Figura 5: Grafico della temperatura T su Γ3 in corrispondenza dell’ottimo (metodo SD).
Risultati confrontabili a quelli finora illustrati si ottengono tramite l’utilizzo del
metodo del gradiente coniugato. Per esempio, scegliendo il metodo PRP (37) per il
calcolo del parametro γ e utilizzando l’opzione 2 per l’aggiornamento di τ abbinato
alla condizione (40) per il ciclo, si ottengono gli andamenti del funzionale J e della sua
derivata J 0 mostrati in figura 6. In questo caso il numero di iterazioni per giungere a
convergenza è 15.
Nella Tabella 1 riassumiamo i risultati finali ottenuti per tutti i metodi disponibili:
il metodo steepest descent con calcolo di τ tramite ciclo iterativo (abbinato alle tre
diverse condizioni (40), (41) e (42) per il ciclo), oppure con calcolo di τ tramite Jott ; il
metodo conjugate gradient con le tre diverse scelte (PRP, PRP+ e CD) per il calcolo
del parametro γ, abbinato alle diverse opzioni per il calcolo di τ . Come si può notare
dalla tabella, i valori di u e J in corrispondenza dell’ottimo che si ottengono sono tutti
molto simili. La principale differenza tra i diversi metodi consiste nel numero di iterazioni
necessarie per raggiungere l’ottimo: tale valore rimane comunque molto limitato, essendo
compreso tra un minimo di 9 ed un massimo di 25, dunque si può affermare che tutti i
metodi utilizzati hanno una velocità di convergenza confrontabile. Per quanto riguarda il
confronto diretto tra il metodo steepest descent e il gradiente coniugato, come si osserva
dai risultati riportati nella tabella, non vi sono differenze significative, o perlomeno esse
non si presentano nel particolare caso qui studiato. I risultati relativi alla scelta della
condizione (42) per l’aggiornamento di τ non vengono riportati, dal momento che con
questa condizione non si raggiunge il valore ottimo: questo potrebbe essere dovuto al
41
Figura 6: Metodo CG: funzionale costo J e derivata J 0 in funzione delle iterazioni.
42
fatto che tale condizione risulta più restrittiva delle altre due condizioni possibili e di
fatto ostacola l’avanzamento del ciclo. Occorre far notare inoltre, come già accennato
nella sezione 3.2, che il metodo per il calcolo di τ che utilizza il valore ottimale del
funzionale, Jott , funziona correttamente solo nel caso in cui viene fornito il valore esatto
preciso di Jott : per valori che si discostano da quello esatto, anche solo dell’ordine del
10%, il metodo impiega molte più iterazioni per giungere a convergenza e in certi casi
può anche accadere che non converga affatto.
Metodo
SD
SD
SD
CG
CG
CG
CG
CG
CG
CG
CG
CG
τ
τ /2
τ /2
Jott
τ /2
τ /2
τ /2
τ /2
τ /2
τ /2
Jott
Jott
2Jott
γ
/
/
/
PRP
PRP
PRP+
PRP+
CD
CD
PRP
PRP+
CD
condizione
1
2
/
1
2
1
2
1
2
/
/
/
n
11
10
17
15
9
17
16
11
9
25
17
12
u
38.2395
38.2397
38.2229
38.229
38.2408
38.2395
38.2397
38.2368
38.2345
38.2206
38.2229
38.2216
J
0.773907
0.773907
0.773908
0.773907
0.773907
0.773907
0.773907
0.773907
0.773907
0.773909
0.773908
0.773908
Tabella 1: Tabella con i risultati: per ogni metodo è riportato il numero di iterazioni n necessarie per
la convergenza, il controllo ottimo u trovato e il corrispondente valore del funzionale costo J. Le scelte
1 e 2 per la condizione del ciclo per τ si riferiscono rispettivamente alla (40) e alla (41).
43
Avendo eseguito tutte le simulazioni numeriche anche utilizzando dei codici FreeFem++, è
possibile effettuare un confronto tra i risultati ottenuti con i due diversi programmi. Per
ogni metodo, il numero di iterazioni totali è pari al numero di iterazioni ottenute con C++,
mentre i valori del parametro ottimale u e del funzionale costo J differiscono leggermente
da quelli già calcolati. Tali differenze sono però dell’ordine di 10−5 e possono quindi, in
relazione ai valori ottenuti, essere considerate trascurabili; anche la distribuzione della
temperatura e della variabile aggiunta presentano solo differenze trascurabili.
Per i risultati illustrati finora, si è utilizzato un parametro di coercività β pari a
10−3 : in questo modo si ottiene la coercività del problema introducendo solo una piccola
perturbazione. Abbiamo però voluto confrontare diversi valori di β: 1, 10−1 , 10−2 , 10−3
e 10−5 . In questo caso si è scelto di utilizzare il metodo SD con dimezzamento di τ
secondo la condizione 40. Man mano che il valore di β diminuisce è sempre più difficile
trovare il minimo: per esempio nel caso già considerato con β = 10−3 bastano 11
iterazioni per raggiungere il valore ottimo di u, mentre per β = 10−5 sono necessarie
ben 598 iterazioni. Questo comportamento è dovuto al fatto che al diminuire del valore
di β il funzionale costo J assume una forma sempre più appiattita ed è quindi più difficile
trovarne il minimo. In figura 7 è possibile notare il diverso andamento del funzionale
costo rispetto ad u, nel caso β = 10−2 e β = 10−3 . Inoltre si osserva che utilizzando
valori di β sempre più piccoli, il valore ottimale di u aumenta, mentre diminuiscono il
valore del funzionale costo J e i valori di temperatura T , che si avvicinano sempre di
più al valore desiderato.
44
Figura 7: Funzionale costo J in funzione del controllo u nel caso β = 10−2 e β = 10−3 .
45
6
Conclusioni
L’obiettivo di questo lavoro è stato quello di applicare la teoria del controllo ottimo
ad un caso reale: il mantenimento della temperatura su di un chip ad un certo livello
desiderato, in modo da evitare il malfunzionamento o la rottura del chip stesso. Per
raggiungere tale scopo il chip è posto a contatto con un dissipatore passivo di calore,
coadiuvato da una ventola nel suo lavoro di dispersione termica.
L’andamento della temperatura sul dissipatore è stato modellato tramite un’equazione
alle derivate parziali (equazione di diffusione con condizioni al contorno di Neumann e
Robin) a cui è stato associato un funzionale costo per tenere conto sia della temperatura
desiderata sul chip, sia del parametro di controllo. Dopo un’analisi matematica teorica del modello formulato, la sua discretizzazione e successiva implementazione tramite
FreeFem++ e C++, hanno permesso di trovare il valore ottimale del parametro di controllo
che permette il corretto funzionamento del chip.
Per la risoluzione numerica del problema sono stati utilizzati diversi metodi di ottimizzazione: diverse versioni del metodo di steepest descent e del metodo del gradiente coniugato. Gli schemi implementati hanno fornito risultati confrontabili e in buon
accordo con la fisica del problema.
Inoltre, questo lavoro prevedeva l’apprendimento del linguaggio di programmazione
C++, conseguito tramite lo sviluppo del codice a elementi finiti necessario per risolvere il
problema fornito dalla teoria del controllo ottimo. L’analisi approfondita di un problema
completo, partendo dall’analisi teorica fino ad arrivare all’implementazione del codice
per le simulazioni, ci ha permesso di comprendere tutte le fasi in cui viene affrontato un
problema a elementi finiti.
46
Appendice
chip sd taumezzi.edp
//**** CHIP - STEEPEST DESCENT - TAU/2 ****
//-è possibile scegliere tra 3 diverse condizioni per il ciclo di aggiornamento di tau
//mesh, spazio FE, incognite, funzione test
mesh Th=readmesh("chip.msh");
plot(Th,wait=1);
fespace Vh(Th,P1);
Vh w,p,wtest,ptest;
real u,u0,uold,
J,J0,Jold,
dJ,dJold,
dJnorm,dJ0norm,dJtol;
//dati problema
real k=237.0;
real h=5.7;
real q=15.08e3;
real winf=25.0;
real wd=65.0;
real b=0.001;
func f=0.0;
//conducibilità termica alluminio [W/m◦ C]
//convezione naturale [W/m^2◦ C]
//flusso di calore [W/m^2]
//temperatura ambiente [◦ C]
//temperatura desiderata [◦ C]
//coefficiente termine coercivo [◦ C^4m^4/W^2]
//forzante
//problema di stato
problem stato(w,wtest) = int2d(Th)(k*(dx(w)*dx(wtest)+dy(w)*dy(wtest)))
-int2d(Th)(f*wtest)
+int1d(Th,1)((h+u)*w*wtest)
-int1d(Th,1)((h+u)*winf*wtest)
-int1d(Th,3)(q*wtest);
//problema aggiunto
problem aggiunto(p,ptest) = int2d(Th)(k*(dx(p)*dx(ptest)+dy(p)*dy(ptest)))
-int2d(Th)(f*ptest)
+int1d(Th,1)((h+u)*p*ptest)
-int1d(Th,3)((w-wd)*ptest);
//inizializzazione algoritmo iterativo
u=0.0;
u0=u;
real taufix=500.0;
real tau=taufix;
real sigma=1.e-3;
real eta=0.5;
real tol=1.e-6;
int n=0;
int nmax=100;
real[int] Jplot(nmax+1);
//controllo iniziale [W/m^2◦ C]
//parametro di accelerazione iniziale
//parametro per seconda condizione
//parametro per terza condizione
//tolleranza
//contatore algoritmo iterativo
//numero massimo iterazioni
//vettore per plot funzionale J
stato;
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T iniziale");
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2; //J(u0)
J0=J;
Jplot[0]=J0;
47
aggiunto;
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p iniziale");
dJ=b*u-int1d(Th,1)((w-winf)*p); //J’(u0)
dJnorm=abs(dJ); //|J’(u0)|
dJ0norm=dJnorm;
dJtol=1;
//per far effettuare il while almeno una volta
//informazioni sull’iterazione 0
cout<<"********iterazione n="<<n<<endl;
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//algoritmo iterativo di ottimizzazione: steepest descent
while(dJtol>tol){
n=n+1;
cout<<"********iterazione n="<<n<<endl;
//memorizzo u, J, dJ al passo precedente
uold=u;
Jold=J;
dJold=dJ;
tau=taufix;
//aggiornamento del controllo u
u=uold-tau*dJold;
//risolvo problema di stato
stato;
//calcolo funzionale costo
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2;
//risolvo problema aggiunto
aggiunto;
//calcolo derivata del funzionale costo J
dJ=b*u-int1d(Th,1)((w-winf)*p);
//aggiornamento parametro di accelerazione tau
//3 diverse possibili condizioni
while (J>Jold){
//while (Jold-J<sigma*tau*dJold^2){
//while (abs(dJold*dJ)>eta*dJold^2){
tau=tau/2;
u=uold-tau*dJold;
cout<<"
tau aggiornato: "<<tau<<endl;
stato;
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2;
aggiunto;
dJ=b*u-int1d(Th,1)((w-winf)*p);
}
Jplot[n]=J;
48
//quantità monitorate per il test d’arresto
dJnorm=abs(dJ);
dJtol=dJnorm/dJ0norm;
//informazioni sull’iterazione corrente
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//check di sicurezza sul numero di iterazioni
if (n>(nmax-1)) break;
};
cout<<"********RISULTATI FINALI********"<<endl;
cout<<"
u0="<<u0<<endl;
cout<<"
uottimo="<<u<<endl;
cout<<"
J0="<<J0<<endl;
cout<<"
Jottimo="<<J<<endl;
cout<<"
numero iterazioni="<<n<<endl;
//plot del funzionale J in funzione delle iterazioni
real[int] aa(n),bb(n);
int i=0;
for(i=0;i<n;i++){
aa[i]=i;
bb[i]=Jplot[i];
}
plot([aa,bb],wait=1);
//plot della variabile di stato w (temperatura)
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T variabile di stato");
//plot della variabile aggiunta p
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p variabile aggiunta");
49
chip sd tauJott.edp
//**** CHIP - STEEPEST DESCENT - TAUJott ****
//mesh, spazio FE, incognite, funzione test
mesh Th=readmesh("chip.msh");
plot(Th,wait=1);
fespace Vh(Th,P1);
Vh w,p,wtest,ptest;
real u,u0,uold,
J,J0,
dJ,dJold,
dJnorm,dJ0norm,dJtol;
//dati problema
real k=237.0;
real h=5.7;
real q=15.08e3;
real winf=25.0;
real wd=65.0;
real b=0.001;
func f=0.0;
//conducibilità termica alluminio [W/m◦ C]
//convezione naturale [W/m^2◦ C]
//flusso di calore [W/m^2]
//temperatura ambiente [◦ C]
//temperatura desiderata [◦ C]
//coefficiente termine coercivo [◦ C^4m^4/W^2]
//forzante
//problema di stato
problem stato(w,wtest) = int2d(Th)(k*(dx(w)*dx(wtest)+dy(w)*dy(wtest)))
-int2d(Th)(f*wtest)
+int1d(Th,1)((h+u)*w*wtest)
-int1d(Th,1)((h+u)*winf*wtest)
-int1d(Th,3)(q*wtest);
//problema aggiunto
problem aggiunto(p,ptest) = int2d(Th)(k*(dx(p)*dx(ptest)+dy(p)*dy(ptest)))
-int2d(Th)(f*ptest)
+int1d(Th,1)((h+u)*p*ptest)
-int1d(Th,3)((w-wd)*ptest);
//inizializzazione algoritmo iterativo
u=0.0;
u0=u;
real tol=1.e-6;
int n=0;
int nmax=100;
real[int] Jplot(nmax+1);
//controllo iniziale [W/m^2◦ C]
//tolleranza
//contatore algoritmo iterativo
//numero massimo iterazioni
//vettore per plot funzionale J
stato;
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T iniziale");
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2; //J(u0)
J0=J;
Jplot[0]=J0;
aggiunto;
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p iniziale");
dJ=b*u-int1d(Th,1)((w-winf)*p); //J’(u0)
dJnorm=abs(dJ); //|J’(u0)|
dJ0norm=dJnorm;
dJtol=1;
//per far effettuare il while almeno una volta
real Jott=0.773943;
//J(uott)
50
real tau=abs(J0-Jott)/(dJ0norm^2); //parametro di accelerazione iniziale
//informazioni sull’iterazione 0
cout<<"********iterazione n="<<n<<endl;
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//algoritmo iterativo di ottimizzazione: steepest descent
while (dJtol>tol){
n=n+1;
//memorizzo u e dJ al passo precedente
uold=u;
dJold=dJ;
//aggiornamento del controllo u
u=uold-tau*dJold;
//risolvo problema di stato
stato;
//calcolo funzionale costo
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2;
Jplot[n]=J;
//risolvo problema aggiunto
aggiunto;
//calcolo derivata del funzionale costo J
dJ=b*u-int1d(Th,1)((w-winf)*p);
//quantità monitorate per il test d’arresto
dJnorm=abs(dJ);
dJtol=dJnorm/dJ0norm;
//informazioni sull’iterazione corrente
cout<<"********iterazione n="<<n<<endl;
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//aggiornamento parametro di accelerazione tau
tau=abs(J-Jott)/(dJnorm^2);
//check di sicurezza sul numero di iterazioni
if (n>(nmax-1)) break;
};
cout<<"********RISULTATI FINALI********"<<endl;
cout<<"
u0="<<u0<<endl;
cout<<"
uottimo="<<u<<endl;
cout<<"
J0="<<J0<<endl;
51
cout<<"
cout<<"
Jottimo="<<J<<endl;
numero iterazioni="<<n<<endl;
//plot del funzionale J in funzione delle iterazioni
real[int] aa(n),bb(n);
int i=0;
for(i=0;i<n;i++){
aa[i]=i;
bb[i]=Jplot[i];
}
plot([aa,bb],wait=1);
//plot della variabile di stato w (temperatura)
plot(w,nbiso=50,value=1,fill=1,wait=1,cmm="T variabile di stato");
//plot della variabile aggiunta p
plot(p,nbiso=50,value=1,fill=1,wait=1,cmm="p variabile aggiunta");
52
chip cg taumezzi.edp
//**** CHIP - CONJUGATE GRADIENT - TAU/2 ****
//-è possibile scegliere tra 3 diverse condizioni per il ciclo di aggiornamento di tau
//-è possibile scegliere tra 3 diversi valori di gamma (PRP, PRP+, CD)
//mesh,spazio FE, incognite, funzione test
mesh Th=readmesh("chip.msh");
plot(Th,wait=1);
fespace Vh(Th,P1);
Vh w,p,wtest,ptest;
real u,u0,uold,
J,J0,Jold,
dJ,dJold,
d,dold,
dJnorm,dJ0norm,dJtol;
//dati problema
real k=237.0;
real h=5.7;
real q=15.08e3;
real winf=25.0;
real wd=65.0;
real b=0.001;
func f=0.0;
//conducibilità termica alluminio [W/m◦ C]
//convezione naturale [W/m^2◦ C]
//flusso di calore [W/m^2]
//temperatura ambiente [◦ C]
//temperatura desiderata [◦ C]
//coefficiente termine coercivo [◦ C^4m^4/W^2]
//forzante
//problema di stato
problem stato(w,wtest) = int2d(Th)(k*(dx(w)*dx(wtest)+dy(w)*dy(wtest)))
-int2d(Th)(f*wtest)
+int1d(Th,1)((h+u)*w*wtest)
-int1d(Th,1)((h+u)*winf*wtest)
-int1d(Th,3)(q*wtest);
//problema aggiunto
problem aggiunto(p,ptest) = int2d(Th)(k*(dx(p)*dx(ptest)+dy(p)*dy(ptest)))
-int2d(Th)(f*ptest)
+int1d(Th,1)((h+u)*p*ptest)
-int1d(Th,3)((w-wd)*ptest);
//inizializzazione algoritmo iterativo
u=0.0;
u0=u;
real taufix=500.0;
real tau=taufix;
real sigma=1.e-3;
real eta=0.3;
real gamma;
real tol=1.e-6;
int n=0;
int nmax=100;
real[int] Jplot(nmax+1);
//controllo iniziale [W/m^2◦ C]
//parametro di accelerazione iniziale
//parametro per seconda condizione
//parametro per terza condizione
//parametro gradiente coniugato
//tolleranza
//contatore algoritmo iterativo
//numero massimo iterazioni
//vettore per plot funzionale J
stato;
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T iniziale");
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2; //J(u0)
J0=J;
Jplot[0]=J0;
aggiunto;
53
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p iniziale");
dJ=b*u-int1d(Th,1)((w-winf)*p); //J’(u0)
d=-dJ;
//d0=-J’(u0)
dJnorm=abs(dJ); //|J’(u0)|
dJ0norm=dJnorm;
dJtol=1;
//per far effettuare il while almeno una volta
//informazioni sull’iterazione 0
cout<<"********iterazione n="<<n<<endl;
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
d="<<d<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//algoritmo iterativo di ottimizzazione: conjugate gradient
while (dJtol>tol){
n=n+1;
cout<<"********iterazione n="<<n<<endl;
//memorizzo u, J, dJ, d al passo precedente
uold=u;
Jold=J;
dJold=dJ;
dold=d;
tau=taufix;
//aggiornamento del controllo u
u=uold+tau*dold;
//risolvo problema di stato
stato;
//calcolo funzionale costo
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2;
//risolvo problema aggiunto
aggiunto;
//calcolo derivata del funzionale costo J
dJ=b*u-int1d(Th,1)((w-winf)*p);
//calcolo gamma, 3 possibili scelte: PRP, PRP+, CD
gamma=dJ*(dJ-dJold)/dJold^2;
//gamma=max(gamma,0.0);
//gamma=-(dJ^2)/(dold*dJold);
//calcolo direzione d
d=-dJ+gamma*dold;
//aggiornamento parametro di accelerazione tau
//3 diverse possibili condizioni
while (J>Jold){
//while (Jold-J<-sigma*tau*dold*dJold^2){
//while (abs(dold*dJ)>eta*abs(dold*dJold)){
tau=tau/2;
54
u=uold+tau*dold;
cout<<"
tau aggiornato: "<<tau<<endl;
stato;
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2;
aggiunto;
dJ=b*u-int1d(Th,1)((w-winf)*p);
gamma=dJ*(dJ-dJold)/dJold^2;
//gamma=max(gamma,0.0);
//gamma=-(dJ^2)/(dold*dJold);
d=-dJ+gamma*dold;
}
Jplot[n]=J;
//quantità monitorate per il test d’arresto
dJnorm=abs(dJ);
dJtol=dJnorm/dJ0norm;
//informazioni sull’iterazione corrente
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
d="<<d<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//check di sicurezza sul numero di iterazioni
if (n>(nmax-1)) break;
};
cout<<"********RISULTATI FINALI********"<<endl;
cout<<"
u0="<<u0<<endl;
cout<<"
uottimo="<<u<<endl;
cout<<"
J0="<<J0<<endl;
cout<<"
Jottimo="<<J<<endl;
cout<<"
numero iterazioni="<<n<<endl;
//plot del funzionale J in funzione delle iterazioni
real[int] aa(n),bb(n);
int i=0;
for(i=0;i<n;i++){
aa[i]=i;
bb[i]=Jplot[i];
}
plot([aa,bb],wait=1);
//plot della variabile di stato w (temperatura)
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T variabile di stato");
//plot della variabile aggiunta p
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p variabile aggiunta");
55
chip cg tauJott.edp
//**** CHIP - CONJUGATE GRADIENT - TAUJott ****
//-è possibile scegliere tra 3 diversi valori di gamma (PRP, PRP+, CD)
//mesh,spazio FE, incognite, funzione test
mesh Th=readmesh("chip.msh");
plot(Th,wait=1);
fespace Vh(Th,P1);
Vh w,p,wtest,ptest;
real u,u0,uold,
J,J0,Jold,
dJ,dJold,
d,dold,
dJnorm,dJ0norm,dJtol;
//dati problema
real k=237.0;
real h=5.7;
real q=15.08e3;
real winf=25.0;
real wd=65.0;
real b=0.001;
func f=0.0;
//conducibilità termica alluminio [W/m◦ C]
//convezione naturale [W/m^2◦ C]
//flusso di calore [W/m^2]
//temperatura ambiente [◦ C]
//temperatura desiderata [◦ C]
//coefficiente termine coercivo [◦ C^4m^4/W^2]
//forzante
//problema di stato
problem stato(w,wtest) = int2d(Th)(k*(dx(w)*dx(wtest)+dy(w)*dy(wtest)))
-int2d(Th)(f*wtest)
+int1d(Th,1)((h+u)*w*wtest)
-int1d(Th,1)((h+u)*winf*wtest)
-int1d(Th,3)(q*wtest);
//problema aggiunto
problem aggiunto(p,ptest) = int2d(Th)(k*(dx(p)*dx(ptest)+dy(p)*dy(ptest)))
-int2d(Th)(f*ptest)
+int1d(Th,1)((h+u)*p*ptest)
-int1d(Th,3)((w-wd)*ptest);
//inizializzazione algoritmo iterativo
u=0.0;
u0=u;
real gamma;
real tol=1.e-6;
int n=0;
int nmax=100;
real[int] Jplot(nmax+1);
//controllo iniziale [W/m^2◦ C]
//parametro gradiente coniugato
//tolleranza
//contatore algoritmo iterativo
//numero massimo iterazioni
//vettore per plot funzionale J
stato;
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T iniziale");
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2; //J(u0)
J0=J;
Jplot[0]=J0;
aggiunto;
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p iniziale");
dJ=b*u-int1d(Th,1)((w-winf)*p); //J’(u0)
d=-dJ;
//d0=-J’(u0)
dJnorm=abs(dJ); //|J’(u0)|
56
dJ0norm=dJnorm;
dJtol=1;
//per far effettuare il while almeno una volta
real Jott=0.773943;
//J(uott)
real tau=abs(J0-Jott)/(dJ0norm^2); //parametro di accelerazione iniziale
//informazioni sull’iterazione 0
cout<<"********iterazione n="<<n<<endl;
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
d="<<d<<endl;
cout<<"
tau="<<tau<<endl;
cout<<"
u="<<u<<endl;
cout<<"
toll="<<dJtol<<endl;
//algoritmo iterativo di ottimizzazione: conjugate gradient
while (dJtol>tol){
n=n+1;
//memorizzo u, dJ, e d al passo precedente
uold=u;
dJold=dJ;
dold=d;
//aggiornamento del controllo u
u=uold+tau*dold;
//risolvo problema di stato
stato;
//calcolo funzionale costo
J=int1d(Th,3)(0.5*((w-wd)^2))+0.5*b*u^2;
Jplot[n]=J;
//risolvo problema aggiunto
aggiunto;
//calcolo derivata del funzionale costo J
dJ=b*u-int1d(Th,1)((w-winf)*p);
//calcolo gamma, 3 possibili scelte: PRP, PRP+, CD
gamma=dJ*(dJ-dJold)/dJold^2;
//gamma=max(gamma,0.0);
//gamma=-(dJ^2)/(dold*dJold);
//calcolo direzione d
d=-dJ+gamma*dold;
//quantità monitorate per il test d’arresto
dJnorm=abs(dJ);
dJtol=dJnorm/dJ0norm;
//informazioni sull’iterazione corrente
cout<<"********iterazione n="<<n<<endl;
cout<<"
J="<<J<<endl;
cout<<"
J’="<<dJ<<endl;
cout<<"
norma di J’="<<dJnorm<<endl;
cout<<"
d="<<d<<endl;
cout<<"
tau="<<tau<<endl;
57
cout<<"
cout<<"
u="<<u<<endl;
toll="<<dJtol<<endl;
//aggiornamento parametro di accelerazione tau
tau=abs(J-Jott)/(dJnorm^2);
//check di sicurezza sul numero di iterazioni
if (n>(nmax-1)) break;
};
cout<<"********RISULTATI FINALI********"<<endl;
cout<<"
u0="<<u0<<endl;
cout<<"
uottimo="<<u<<endl;
cout<<"
J0="<<J0<<endl;
cout<<"
Jottimo="<<J<<endl;
cout<<"
numero iterazioni="<<n<<endl;
//plot del funzionale J in funzione delle iterazioni
real[int] aa(n),bb(n);
int i=0;
for(i=0;i<n;i++){
aa[i]=i;
bb[i]=Jplot[i];
}
plot([aa,bb],wait=1);
//plot della variabile di stato w (temperatura)
plot(w,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="T variabile di stato");
//plot della variabile aggiunta p
plot(p,nbiso=50,value=1,fill=1,boundary=0,wait=1,cmm="p variabile aggiunta");
58
Riferimenti bibliografici
[Ç98] Y. A. Çengel, Termodinamica e Trasmissione del Calore, McGraw-Hill, Milano,
1998.
[D08] L. Dedè, Adaptive and Reduced Basis Methods for Optimal Control Problems in
Environmental Applications, Tesi di dottorato, Politecnico di Milano, 2008.
[EG04] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, SpringerVerlag, New York, 2004.
[Q08] A. Quarteroni, Modellistica Numerica per Problemi Differenziali, Springer-Verlag
Italia, Milano, 2008.
[QSS00] A. Quarteroni, R. Sacco, F. Saleri, Matematica Numerica, Springer-Verlag
Italia, Milano, 2000.
[S07] S. Salsa, Elementi di Ottimizzazione e Teoria del Controllo, Dispense del corso di
EDP2, Politecnico di Milano, 2007.
[Y01] D. Yang, C++ and Object Oriented Numeric Computing for Scientists ane
Engineers, Springer-Verlag, New York, 2001.
[C++] C++.com, web page, http://www.cplusplus.com.
[FreeFem++] FreeFem++.org, web page, http://freefem.org.
59