POLITECNICO DI MILANO

CORSO DI LAUREA IN INGEGNERIA MATEMATICAORIENTAMENTO 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. 708835FRANCESCA STORTI Matr. 708131

Anno Accademico 2008-2009

Indice

1 Introduzione 3

2 Problema in esame 42.1 Equazioni del modello . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42.2 Parametri del modello . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62.3 Analisi di buona posizione . . . . . . . . . . . . . . . . . . . . . . . . . . 8

3 Problema di controllo ottimo 113.1 Analisi di esistenza e unicita . . . . . . . . . . . . . . . . . . . . . . . . . 153.2 Risoluzione numerica: metodo iterativo . . . . . . . . . . . . . . . . . . 21

4 Implementazione del codice 254.1 Pre-processing e mesh computazionale . . . . . . . . . . . . . . . . . . . 254.2 Integrazione numerica . . . . . . . . . . . . . . . . . . . . . . . . . . . . 284.3 Assemblaggio e risoluzione dei sistemi lineari . . . . . . . . . . . . . . . 314.4 Ciclo iterativo . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 364.5 Post-processing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 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 unaventola, 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 valoridi temperatura: il superamento di tali valori puo causare un malfunzionamento del chipstesso fino a portarlo a rottura. Lo strumento matematico utilizzato per raggiungere talescopo e la teoria del controllo ottimo per equazioni alle derivate parziali. L’andamentodella temperatura viene modellato tramite un’equazione di diffusione con opportunecondizioni al contorno, mentre il parametro di cui si vuole trovare il valore ottimaleviene modellato tramite un funzionale costo ad esso associato.

In particolare, il sistema fisico e costituito da tre componenti: il chip, affiancato daun dissipatore passivo di calore e dalla ventola. Il dissipatore, posto a diretto contattocon il chip, permette di disperdere nell’ambiente circostante il calore ricevuto dal chiptramite flusso conduttivo. Il processo dissipativo avviene pero sia per convezione liberache forzata: il lavoro del dissipatore e infatti coadiuvato dalla ventola, che accelera ilraffredddamento dell’intero sistema.

Trovare il punto di funzionamento ottimale della ventola significa trovare la suapotenza e la sua regolazione ideali, in base alla temperatura di lavoro del chip e alla ca-pacita di dispersione termica del dissipatore. L’EDP considerata modella la temperaturaall’interno del dissipatore passivo; il funzionale costo da minimizzare riguarda invece latemperatura sul chip, che non deve discostarsi troppo dal suo valore di funzionamentoideale; infine, il parametro di cui vogliamo trovare il valore ottimale e il coefficientedi convezione forzata, direttamente collegato al punto di funzionamento ottimale dellaventola.

Nella sezione 2 presentiamo il modello fisico considerato con tutti i parametri ad essoassociati, insieme all’equazione di diffusione che descrive l’andamento della temperaturaed alla sua analisi di buona posizione. Nella sezione 3 presentiamo il vero e proprioproblema di controllo ottimo: descriviamo l’espressione del funzionale costo e ricaviamole equazioni che legano tale funzionale all’EDP per la temperatura; presentiamo inoltrel’analisi di esistenza e unicita della variabile di controllo ottimo. Nella sezione 4 pre-sentiamo diversi metodi numerici necessari per la risoluzione del problema di controlloottimo, in particolare descriviamo il ciclo iterativo che permette il raggiungimento delvalore ottimale della variabile di controllo. Nella sezione 5 presentiamo il codice da noisviluppato tramite il linguaggio di programmazione C++ per la risoluzione del problemadi controllo. Inizialmente abbiamo pero implementato il problema tramite FreeFem++,in modo da poter effettuare un confronto tra i risultati ottenuti con i due diversi pro-grammi; i programmi FreeFem++ sono riportati in Appendice. Infine nella sezione 6riportiamo e commentiamo i principali risultati ottenuti.

3

2 Problema in esame

2.1 Equazioni del modello

Per studiare la propagazione del calore e il campo di temperatura all’interno del dissipa-tore, consideriamo i diversi scambi termici che avvengono nell’intero sistema. All’internodel dissipatore il calore si propaga tramite conduzione termica. La distribuzione di tem-peratura nel dissipatore, in regime stazionario, e quindi descritta dalla seguente legge didiffusione:

−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: laquantita di calore scambiata dipende dalla geometria e dalle caratteristiche del corpo edalla differenza di temperatura tra i due diversi mezzi. In regime stazionario il flusso ter-mico, ovvero la potenza termica per unita di superficie, trasmesso dal chip al dissipatoreper 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. Sipuo distinguere tra due diversi tipi di convezione: convezione naturale se lo scambio dicalore avviene per una semplice differenza di temperatura e convezione forzata se invecelo scambio di calore e favorito dall’azione di dispositivi esterni che generano trasportodi massa, come ad esempio la ventola nel caso in esame. Il flusso termico legato allaconvezione e direttamente proporzionale alla differenza di temperatura, come espressodalla 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 (nelnostro caso aria) a distanza sufficientemente grande dal dissipatore. Il coefficiente hnon e una proprieta del fluido, ma e un parametro determinato sperimentalmente, il cuivalore dipende da tutte le variabili che influenzano il processo di convezione, quali lageometria della superficie di scambio, la natura del moto, le proprieta e la velocita delfluido circostante. Tale parametro puo essere legato alla semplice convezione naturaleoppure puo includere il contributo di un’eventuale convezione forzata (come nel nostrocaso).

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 volumedel dissipatore, con Γ1 la superficie a contatto con l’aria circostante, interessata allaconvezione, con Γ2 la superficie a contatto con il supporto solido del dispositivo, con-siderata adiabatica, cioe che non permette la trasmissione del calore verso l’ambienteesterno (attraverso la quale il flusso di calore q e nullo) e con Γ3 la superficie a contattocon 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ı daaumentare la quantita di calore scambiato per convezione. Il chip e fissato alla superficieorizzontale 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

k∇T · n = 0 su Γ2

k∇T · n = q su Γ3

, (1)

dove T = T (x, y) e la temperatura, funzione a valori scalari delle due variabili spaziali xe y, k e la conducibilita termica dell’alluminio, materiale di cui e costituito il dissipatorein esame, q e il flusso di calore tra chip e dissipatore, n e il vettore normale esterno allediverse parti del dominio, T∞ e la temperatura dell’ambiente circostante, (h + u) e ilcoefficiente di convezione forzata dell’aria (infatti, essendo h, coefficiente di convezionenaturale dell’aria, molto piccolo rispetto ad u, coefficiente di convezione forzata, possi-amo 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 ladistribuzione di temperatura all’interno del dissipatore; la seconda e una condizione diRobin e definisce il flusso di calore scambiato per convezione tra il dissipatore e l’ariacircostante; la terza e la quarta sono condizioni di Neumann e definiscono rispettiva-mente il flusso di calore nullo attraverso la superficie adiabatica del dissipatore e il flussodi 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 buonaposizione (esistenza e unicita della soluzione).

2.2 Parametri del modello

Vediamo ora in dettaglio il calcolo dei parametri fisici presenti nel modello espresso dalsistema (1).

Il chip considerato e il tipo T-5000 della Intel, mentre il dissipatore e una versionesemplificata del tipo HS 6115 della Vemaline Products, Inc. ([C98, pag. 360]): esso e undispositivo anodizzato nero sagomato di alluminio, costituito da 10 alette di lunghezza25 mm, larghezza 76 mm e spessore 3 mm. La lunghezza totale del dispositivo e di102 mm.

Le principali proprieta termodinamiche del sistema chip-dissipatore sono evidenziatedal coefficiente di conducibilita termica del dissipatore k e dal coefficiente di convezionenaturale dell’aria h.La conducibilita termica k e quella propria dell’alluminio, cioe k = 237 W

mC .La convezione naturale h riferita all’aria dipende dalla particolare geometria del dispos-itivo considerato:

h = 1.31karia

Sopt, (2)

dove karia = 0.026 WmC e la conducibilita termica dell’aria e Sopt indica la spaziatura

ottimale tra le alette del dissipatore. Tale spaziatura puo essere calcolata tramite laseguente formula:

Sopt = 2.714L

Ra, (3)

dove L e la lunghezza caratteristica del dissipatore, corrispondente nel nostro casoalla 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

ν2Pr , (4)

dove:

• g = 9.81 ms2 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 sullasuperficie del dispositivo e la temperatura dell’aria sufficientemente lontano dallasuperficie,

• δ e una lunghezza caratteristica, coincidente in questo caso con L,

• ν = 1.57× 10−5 m2

s e la viscosita cinematica dell’aria,

• Pr = 0.712 e il numero di Prandtl (adimensionale).

6

Consideriamo T∞ coincidente con la temperatura ambiente, pari cioe a 20C, e Ts parialla temperatura massima di funzionamento del chip, cioe 90C; γ puo essere calcolatocome l’inverso della temperatura media tra T∞ e Ts misurate in K, cioe nel seguentemodo:

γ =2

(Ts + T∞). (5)

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’ariah, che risulta pari a:

h = 5.7W

m2C.

Il coefficiente di convezione forzata dell’aria u, che e direttamente collegato alle carat-teristiche termodinamiche della ventola, e il parametro che possiamo variare in modoche il chip non si riscaldi troppo; il suo valore e solitamente compreso tra 25 W

m2Ce

250 Wm2C

.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 flussodi calore disperso dal dissipatore Integrando la prima equazione del sistema (1) su tuttoil dominio Ω e sfruttando il teorema della divergenza otteniamo:

−∫

Ω∇ · (k∇T )dΩ = −

∫∂Ωk∇T · ndΓ = 0 ,

che puo essere riscritto separatamente sulle tre diverse porzioni di bordo del dominiocome:

−∫

Γ1

k∇T · ndΓ−∫

Γ2

k∇T · ndΓ−∫

Γ3

k∇T · ndΓ = 0 ,

e sostituendo le condizioni al contorno del sistema (1):∫Γ1

(h+ u)(T − T∞)dΓ−∫

Γ3

qdΓ = 0 .

Quest’ultima equazione puo essere riscritta nel seguente modo:∫Γ3

qdΓ =∫

Γ1

(h+ u)(T − T∞)dΓ ,

e ricordando che le quantita q e (h + u) sono costanti, possiamo ricavare l’espressioneper il calcolo del flusso di calore q:

q =(h+ u)|Γ3|

∫Γ1

(T − T∞)dΓ ,

dove |Γ3| rappresenta la lunghezza del bordo Γ3. Scegliamo a questo punto un valoreragionevole per il coefficiente di convezione forzata pari a (h+ u) = 25 W

m2Ce un valore

7

medio per la temperatura pari a T = 90C: il valore T − T∞ diventa cosı una costante,per cui l’integrale puo essere riscritto nel seguente modo:

q =(h+ u)|Γ3|

(T − T∞) |Γ1| ,

dove |Γ1| rappresenta la lunghezza del bordo Γ1. Sostituendo i valori dei parametripresenti nell’ultima equazione otteniamo infine il valore del flusso termico q:

q = 15080W

m2.

2.3 Analisi di buona posizione

Ci occupiamo ora dell’analisi di buona posizione del problema ellittico di partenzaespresso dal sistema (1), prima di affrontare il problema di controllo ottimo vero eproprio.

Gli spazi funzionali che utilizzeremo sono: H1(Ω) e R, in particolare T ∈ H1(Ω) eu ∈ R. Prima di affrontare il problema di controllo vero e proprio, verifichiamo la buonaposizione del problema ellittico di partenza.

Ricaviamo innanzitutto la formulazione debole del problema (1). Consideriamo quin-di una funzione test v ∈ H1(Ω) e moltiplichiamo l’equazione differenziale del sistema(1) per tale funzione:

−∫

Ω∇ · (k∇T )vdΩ = 0 .

Applichiamo la regola di integrazione per parti, ottenendo:

−∫

∂Ωk∇T · nvdΓ +

∫Ωk∇T · ∇vdΩ = 0 .

Sostituendo le condizioni al contorno del sistema (1) e ricordando che h, u e q sonograndezze costanti, otteniamo:

(h+ u)∫

Γ1

(T − T∞)vdΓ− q

∫Γ3

vdΓ +∫

Ωk∇T · ∇vdΩ = 0 ,

da cui, portando al secondo membro gli integrali che non contengono l’incognita T :∫Ωk∇T · ∇vdΩ + (h+ u)

∫Γ1

TvdΓ = (h+ u)T∞∫

Γ1

vdΓ + q

∫Γ3

vdΓ .

A questo punto definiamo la forma bilineare a : H1(Ω)×H1(Ω) → R:

a(T, v) =∫

Ωk∇T · ∇vdΩ + (h+ u)

∫Γ1

TvdΓ , (6)

e il funzionale lineare F : H1(Ω) → R:

F (v) = (h+ u)T∞∫

Γ1

vdΓ + q

∫Γ3

vdΓ . (7)

8

La formulazione debole del problema (1) risulta quindi essere:

trovare T ∈ H1(Ω) : a(T, v) = F (v) ∀v ∈ H1(Ω) . (8)

A questo punto per dimostrare l’esistenza e l’unicita della soluzione del problema scrittonella forma debole (8), utilizziamo il lemma di Lax-Milgram:

Lemma 2.1. Sia V uno spazio di Hilbert, a(·, ·) : V × V → R una forma bilinearecontinua e coerciva. Sia F (·) : V → R un funzionale lineare e continuo. Allora esisteunica 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 consider-

azione i due membri della forma bilineare separatamente e applichiamo ad entrambi ladisuguaglianza di Schwarz :∣∣∣∣∫

Ωk∇T · ∇vdΩ

∣∣∣∣ ≤ k ‖∇T‖L2(Ω) ‖∇v‖L2(Ω) ,∣∣∣∣(h+ u)∫

Γ1

TvdΓ∣∣∣∣ ≤ (h+ u) ‖T‖L2(Γ1) ‖v‖L2(Γ1) ,

ricordando che i parametri k, h ed u sono costanti e positivi. Utilizziamo poi ladefinizione di norma nello spazio H1(Ω) per il primo membro e la disuguaglianza ditraccia (con costante C) per il secondo:∣∣∣∣∫

Ωk∇T · ∇vdΩ

∣∣∣∣ ≤ k ‖T‖H1(Ω) ‖v‖H1(Ω) ,∣∣∣∣(h+ u)∫

Γ1

TvdΓ∣∣∣∣ ≤ C2(h+ u) ‖T‖H1(Ω) ‖v‖H1(Ω) .

Dalle ultime due disequazioni si ricava quindi:

|a(T, v)| ≤ maxk,C2(h+ u)

‖T‖H1(Ω) ‖v‖H1(Ω) . (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 operatorelineare, continuo e iniettivo e B : X → Z un operatore lineare, continuo e compatto. Seesiste una costante c > 0 tale che c ‖x‖X ≤ ‖Ax‖Y + ‖Bx‖Z , allora esiste una costantec tale che:

c ‖x‖X ≤ ‖Ax‖Y ∀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 ∂Ω sialipschitziano) e sia f : H1(Ω) → R un funzionale lineare tale che f(v) 6= 0 ∀v 6= 0,v ∈ R. Allora esiste una costante positiva cR tale che:

cR ‖v‖H1(Ω) ≤ ‖∇v‖L2(Ω) + |f(v)| . (10)

Dalla (10), detta anche disuguaglianza di Poincare-Friedrichs, elevando al quadrato esfruttando la disuguaglianza di Young si ottiene:

c2R ‖v‖2H1(Ω) ≤ (‖∇v‖L2(Ω) + |f(v)|)2

= ‖∇v‖2L2(Ω) + |f(v)|2 + 2 ‖∇v‖L2(Ω) |f(v)|

≤ 2(‖∇v‖2

L2(Ω) + |f(v)|2). (11)

Per trattare il nostro problema di partenza, ovvero il caso di equazione di diffusione concondizione di Robin e Neumann, scegliamo opportunamente come funzionale f(v) =(h+u)1/2

k1/2 ‖v‖2L2(Γ1). Riscriviamo quindi la (11) utilizzando l’espressione del funzionale:

c2R ‖v‖2H1(Ω) ≤ 2

(‖∇v‖2

L2(Ω) +h+ u

k‖v‖2

L2(Γ1)

),

che equivale alla seguente:

k ‖∇v‖2L2(Ω) + (h+ u) ‖v‖2

L2(Γ1) ≥k

2c2R ‖v‖

2H1(Ω) . (12)

Dalla definizione della forma bilineare (6) ricaviamo che:

a(v, v) = k ‖∇v‖2L2(Ω) + (h+ u) ‖v‖2

L2(Γ1) ,

che combinata con la (12) ci fornisce la seguente:

a(v, v) ≥ k

2c2R ‖v‖

2H1(Ω) . (13)

La forma bilineare e dunque coerciva.Per dimostrare la continuita del funzionale F (·), utilizziamo nell’ordine la disug-

uaglianza triangolare e la disuguaglianza di Cauchy-Schwarz :

|F (v)| ≤∣∣∣∣(h+ u)T∞

∫Γ1

vdΓ∣∣∣∣ +

∣∣∣∣q ∫Γ3

vdΓ∣∣∣∣

≤ (h+ u)T∞ ‖v‖L2(Γ1) + q ‖v‖L2(Γ3) ,

ricordando che i parametri h, u, T∞ e q sono costanti e positivi. Applichiamo infine ladisuguaglianza di traccia (con costanti C e C) e otteniamo:

|F (v)| ≤[(h+ u)T∞C + qC

]‖v‖H1(Ω) . (14)

Il funzionale e dunque continuo.Grazie alle (9), (13) e (14), sono quindi soddisfatte le ipotesi del lemma di Lax-

Milgram, 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 prob-lema (1), e finalizzato al controllo della temperatura raggiunta dal chip, cioe quellarelativa al bordo Γ3 (superficie di contatto tra il chip e il dissipatore). Questa temper-atura non puo infatti superare un certo valore di soglia stabilito a priori, oltre cui none piu garantito il corretto funzionamento del dispositivo. Per tenere conto di questacondizione, introduciamo il seguente funzionale costo, che rappresenta la differenza trala temperatura effettiva T a cui si trova il chip e la sua temperatura di lavoro idealemedia Td:

J(u) =12

∫Γ3

(T − Td)2dΓ +12βu2 , (15)

dove il primo termine valuta appunto il quadrato della differenza tra le due temperatureT e Td sulla superficie Γ3 di contatto tra il dissipatore e il chip, mentre il secondotermine pesato con un opportuno coefficiente β (positivo e sufficientemente piccolo)e contenente la variabile di controllo u, e un termine opportunamente aggiunto pergarantire la coercivita del problema, condizione necessaria per la buona posizione dellostesso.Unendo il sistema di partenza (1) all’espressione del funzionale costo (15) si ottiene unproblema di controllo ottimo. In termini generali per controllo si intende la possibilitadi 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 puoessere 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 ilparametro di convezione forzata;

• le variabili che descrivono lo stato del sistema, dette funzioni di stato, definite inun opportuno spazio funzionale e dipendenti dal valore assunto dal controllo u; nelnostro caso si tratta della temperatura T ∈ H1(Ω);

• l’insieme di equazioni e relative condizioni al bordo che costituiscono il cosiddettoproblema di stato, il quale descrive il sistema fisico che si intende controllare eche deve essere soddisfatto dalle variabili di stato e di controllo; in forma generalel’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 eche agisce sullo stato T , f e il termine noto e B e un altro operatore che introduceil 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 casodalla (15);

11

• la funzione di osservazione, dipendente dalla variabile di stato e quindi a sua voltadal controllo, definita anch’essa in un opportuno spazio funzionale, va confrontatacon un certo valore desiderato stabilito a priori che rappresenta l’obiettivo daraggiungere attraverso la minimizzazione del funzionale costo; nel nostro caso sitratta 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 euno stato T che minimizzino il funzionale J(u) sotto il vincolo dato dalle equazionidel problema di stato (1). Per la risoluzione di tale problema utilizziamo il metododella lagrangiana: attraverso l’uso dei moltiplicatori di Lagrange ci si riconduce ad unproblema di ottimizzazione libera, si cercano cioe i punti critici liberi del funzionalelagrangiano definito in generale nel modo seguente:

L(T, p, u) = J(u) + 〈p, f +Bu−AT (u)〉 , (17)

dove p ∈ H1(Ω) rappresenta il moltiplicatore di Lagrange (o variabile aggiunta), ilsimbolo 〈·, ·〉 denota la dualita tra H1(Ω) e H−1(Ω) e si e fatto riferimento alla for-ma generale dell’equazione di stato (16). Il funzionale lagrangiano equivale quindi alfunzionale costo con l’aggiunta di un termine, formalmente nullo, che rappresenta ilprodotto tra la variabile aggiunta e il vincolo dato dall’equazione di stato espressa informa implicita. Nel nostro caso otteniamo:

L(T, p, u) = J(u) + 〈p,∇ · (k∇T )〉 . (18)

Il problema di ottimizzazione libera risulta quindi essere:

trovare (T, p, u) : ∇L(T, p, u) = 0 , (19)

esprimibile anche come: LT = 0

Lp = 0

Lu = 0

, (20)

dove si e usata la notazione abbreviata LT per indicare la derivata di Gateaux del fun-zionale L rispetto alla variabile T ; analogo significato hanno le notazioni Lp e Lu. Le treequazioni del sistema (20) forniranno rispettivamente l’equazione aggiunta, l’equazionedi stato e la condizione di ottimalita, mediante le quali determineremo lo stato, la vari-abile aggiunta e il controllo ottimi.Vediamo ora, secondo due diversi approcci, come ottenere queste tre equazioni checostituiscono 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 itermini legati al termine noto f e all’operatore B non sono presenti, in quanto nullinel nostro problema. Nel nostro caso il funzionale lagrangiano in forma debole risultaquindi essere:

L(T, p, u) =12

∫Γ3

(T − Td)2dΓ +12βu2

−∫

Ωk∇T · ∇pdΩ− (h+ u)

∫Γ1

(T − T∞)pdΓ + q

∫Γ3

pdΓ . (21)

A questo punto il problema di controllo diventa:

trovare (T, p, u) ∈ (H1(Ω)×H1(Ω)× R) :

∇L(T, p,u)[ϕ, φ, ψ] = 0 ∀(ϕ, φ, ψ) ∈ (H1(Ω)×H1(Ω)× R) ,

traducibile, tramite le derivate di Gateaux, nel sistema di tre equazioni:LT = 0 ∀ϕ ∈ H1(Ω)

Lp = 0 ∀φ ∈ H1(Ω)

Lu = 0 ∀ψ ∈ R

. (22)

Calcolando le tre equazioni del sistema (22), annullando cioe le derivate di Gateauxdel funzionale lagrangiano espresso dalla (21) rispetto alle tre variabili (considerate inquesto caso indipendenti tra loro), otteniamo la formulazione cercata per il problema dicontrollo. L’equazione di stato in forma debole risulta:

Lp[φ] = −∫

Ωk∇T · ∇φdΩ− (h+ u)

∫Γ1

(T − T∞)φdΓ + q

∫Γ3

φdΓ = 0 ∀φ ∈ H1(Ω) .

(23)L’equazione aggiunta in forma debole risulta:

LT [ϕ] = −∫

Ωk∇ϕ · ∇pdΩ− (h+ u)

∫Γ1

ϕpdΓ +∫

Γ3

(T − Td)ϕdΓ = 0 ∀ϕ ∈ H1(Ω) ,

che, cambiando i segni e integrando per parti il primo termine, diventa:

LT [ϕ] = −∫

Ωϕ∇ · (k∇p)dΩ +

∫∂Ωkϕ∇p · ndΓ + (h+ u)

∫Γ1

ϕpdΓ

−∫

Γ3

(T − Td)ϕdΓ = 0 ∀ϕ ∈ H1(Ω) . (24)

Dalla (24) si ricava quindi il problema aggiunto:

−∇ · (k∇p) = 0 in Ω

k∇p · n + (h+ u)p = 0 su Γ1

k∇p · n = 0 su Γ2

k∇p · n− (T − Td) = 0 su Γ3

. (25)

13

Infine l’equazione che esprime il principio di minimo (o condizione di ottimalita) risulta:

Lu[ψ] =βuψ − ψ

∫Γ1

(T − T∞) pdΓ

=(βu−

∫Γ1

(T − T∞)pdΓ)ψ ≥ 0 ∀ψ ∈ R . (26)

Le equazioni (23), (24) e (26) sono quelle che andranno risolte per risolvere il problemadi ottimo.

L’approccio basato sul funzionale lagrangiano puo essere sviluppato anche in unaltro 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 : H1(Ω) × R(Γ1) → H−1(Ω). Per tenere conto di tale vincolo, scriviamo ilfunzionale costo nella seguente forma:

J(T ) =12

∫Γ3

(T − Td)2dΓ +12βu2 +

∫Ωp∇ · (k∇T )dΩ ,

dove abbiamo aggiunto al funzionale costo originale (15) un termine praticamente nullo.La variabile p appare come un moltiplicatore a nostra disposizione, che scegliamo dipen-dente solo dalla variabile di controllo u. Considerando questo approccio, la condizionedi ottimalita J ′(u)(v − u) ≥ 0 diventa:∫

Γ3

(T (u)− Td) (T (v)− T (u)) dΓ+βu(v−u)+∫

Ωp∇·(k∇ (T (v)− T (u))) dΩ ≥ 0 . (27)

Integriamo due volte per parti l’ultimo termine della (27):∫Ωp∇ · (k∇(T (v)− T (u))) dΩ

=∫

∂Ωkp∇(T (v)− T (u)) · n dΓ−

∫Ωk∇p · ∇(T (v)− T (u))dΩ

=∫

∂Ωkp∇(T (v)− T (u)) · n dΓ−

∫∂Ωk(T (v)− T (u))∇p · n dΓ

+∫

Ωk(T (v)− T (u))∆p dΩ . (28)

Consideriamo ora il primo termine della (28) e sostituiamo in esso le condizioni alcontorno del problema di partenza (1):∫

∂Ωkp∇(T (v)− T (u)) · n dΓ =

∫Γ1

h(T (u)− T (v))dΓ +∫

Γ1

T∞(v − u)dΓ

+∫

Γ1

(uT (u)− vT (v))dΓ . (29)

14

Riscriviamo il termine non lineare (l’ultimo della (29)) nel seguente modo:∫Γ1

(uT (u)− vT (v))dΓ =∫

Γ1

(u− v)T (u)dΓ +∫

Γ1

u(T (u)− T (v))dΓ

+∫

Γ1

(v − u)(T (u)− T (v))dΓ , (30)

in modo da separare le variabili u e v da T (u) e T (v), rispettivamente. Raggruppandotutti i termini trovati, la condizione di ottimalita (27) diventa:∫

Γ3

(T (u)− Td)(T (v)− T (u))dΓ + βu(v − u)

+∫

Γ1

h(T (u)− T (v))dΓ +∫

Γ1

T∞(v − u)dΓ

+∫

Γ1

(u− v)T (u)dΓ +∫

Γ1

u(T (u)− T (v))dΓ +∫

Γ1

(v − u)(T (u)− T (v))dΓ

−∫

∂Ω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 su Γ2

k∇p · n− (T (u)− Td) = 0 su Γ3

,

e il principio di minimo:(βu−

∫Γ1

(T − T∞)p dΓ)

(v − u) ≥ 0 ∀v ∈ R .

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:∫

Γ1

J ′(u)(v − u)dΓ ≥ 0 ∀v ∈ R ,

dove J ′(u) indica la derivata del funzionale costo; tale forma e quella che viene effetti-vamente 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 unproblema di minimo vincolato (si veda [S07]), esprimibile nel caso generale come:

min F (x, u)con il vincolo:G(x, u) = 0 ,

(32)

15

dove F : X ×U → R e un funzionale della forma F (x, u) = P (x) +N(u) e G e un altrofunzionale tale che G : X × U → Y . X, Y e U sono spazi di Banach riflessivi e U eun covesso chiuso (quindi anche debolmente chiuso). Diciamo che (x, u) e soluzione delproblema (57) se e un punto di minimo vincolato per F . In particolare, vale il seguenteteorema di esistenza:

Teorema 3.1. Assumiamo che:

1. infx∈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 eG(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 cuilavoriamo sono: X = H1(Ω), U = R(Γ1) e Y = H−1(Ω). La variabile x corrisponde allatemperatura 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 dueparti P (T ) = 1

2

∫Γ3

(T −Td)2dΓ e N(u) = 12βu

2. Il vincolo G(T, u) e l’equazione ellitticada 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 momento che la funzione integranda e elevata al quadrato, il funzionale P (T ) =12

∫Γ3

(T − Td)2dΓ e sempre non negativo, quindi il suo valore minimo e raggiuntoquando la funzione T−Td si annulla, cioe quando T = Td. Di conseguenza abbiamoche:

infT∈H1(Ω)

P = 0 ≥ −∞.

Inoltre, dall’analisi di buona posizione del problema (10), sappiamo che l’equazioneG(T, u) = 0 ammette soluzione unica.

2. La seconda ipotesi richiede la coercivita del funzionale N , facilmente dimostrabilegrazie alla presenza del parametro β > 0:

N(u) = 12β

∫Γ1u2dΓ = β|u|2,

avendo inglobato il fattore 12 all’interno del parametro β e avendo indicato con | · |

la norma in R(Γ1), cioe il valore assoluto. Dall’ultima equazione si vede quindiche vale l’uguaglianza, con γ = 2.

16

3. Consideriamo le due successioni debolmente convergenti:

uk u in R(Γ1)

Tk T in H1(Ω),

dove (Tk, uk) ∈ H1(Ω)× R(Γ1). Vogliamo dimostrare che:

G(Tk, uk) G(T, u) in H−1(Ω),

cioe che:

−∫Ω∇ · (k∇Tk)vdΩ → −

∫Ω∇ · (k∇T )vdΩ ∀v ∈ H1(Ω).

Integriamo per parti due volte:

−∫Ω∇ · (k∇Tk)vdΩ

e sostituiamo le condizioni al contorno:

−∫

Ω∇ · (k∇Tk)vdΩ =

∫Ωk∇Tk∇vdΩ−

∫∂Ωk∂Tk

∂nvdΓ

=∫

∂ΩkTk

∂v

∂ndΓ−

∫ΩkTk∆vdΩ−

∫∂Ωk∂Tk

∂nvdΓ

=∫

∂ΩkTk

∂v

∂ndΓ−

∫ΩkTk∆vdΩ +

∫Γ1

(h+ uk)(Tk − T∞)vdΓ+

−∫

Γ3

qvdΓ .

Consideriamo adesso gli ultimi quattro integrali separatamente:

•∫∂Ω kTk

∂v∂ndΓ →

∫∂Ω kT

∂v∂ndΓ

per la convergenza debole della successione Tk.• −

∫Ω kTk∆vdΩ → −

∫Ω kT∆vdΩ

per la convergenza debole della successione Tk.•

∫Γ1

(h+ uk)(Tk − T∞)vdΓ puo essere scomposto in altri quattro integrali:

(a)∫Γ1hTkvdΓ →

∫Γ1hTvdΓ

per la convergenza debole della successione Tk.(b) −

∫Γ1hT∞vdΓ e una costante e rimane inalterata nel passaggio al limite

per k →∞.(c)

∫Γ1ukTkvdΓ e una forma trilineare nelle variabili u, T e v.

(d) −∫Γ1ukT∞vdΓ e una costante e rimane inalterata nel passaggio al limite

per k →∞.Analizziamo piu in dettaglio la forma trilineare (c); essa puo essere riscrittanel seguente modo:

(c)∫

Γ1

ukTkvdΓ =∫

Γ1

(ukTk ± uTk ± uT )vdΓ

=∫

Γ1

[(uk − u)Tk + (Tk − T )u+ Tu)]vdΓ

17

e, dal momento che uk u in R(Γ1) e Tk T in H1(Ω):∫Γ1

(uk − u)TkvdΓ → 0∫Γ1

(Tk − T )uvdΓ → 0,quindi:∫Γ1TkukvdΓ →

∫Γ1TuvdΓ.

• −∫Γ3qvdΓ e una costante e rimane inalterata nel passaggio al limite per

k →∞.

In conclusione:

−∫

Ω∇ · (k∇Tk)vdΩ →

∫∂ΩkT

∂v

∂ndΓ−

∫ΩkT∆vdΩ +

∫Γ1

hTvdΓ+

−∫

Γ1

hT∞vdΓ +∫

Γ1

uTvdΓ−∫

Γ1

hT∞vdΓ−∫

Γ3

qvdΓ ,

che, controintegrando per parti, significa:

−∫Ω∇ · (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 infx→x0

P (x) + lim infu→u0

N(u)

e F (x, u) si dice debolmente semicontinuo inferiormente in X×U se e debolmentesemicontinuo inferiormente in ogni suo punto (x0, u0) ∈ X × U .

Nel nostro caso, dato (T0, u0) ∈ H1(Ω)× R(Γ1), possiamo dire che:

F (T0, u0) = P (T0) +N(u0)

=12

∫Γ3

(T0 − Td)2dΓ +12βu2

0

=12

∫Γ3

limT→T0

(T − Td)2dΓ +12β lim

u→u0

u2

=12

∫Γ3

lim infT→T0

(T − Td)2dΓ +12β lim inf

u→u0

u2.

A questo punto possiamo applicare il lemma di Fatou al funzionale P (T ) valutatoin T = T0, dal momento che (T − Td)2 e una successione di funzioni misurabili avalori non negativi:

12

∫Γ3

lim infT→T0

(T − Td)2dΓ +12β lim inf

u→u0

u2 ≤

18

≤ 12

lim infT→T0

∫Γ3

(T − Td)2dΓ +12β lim inf

u→u0

u2,

quindi:

F (T0, u0) = P (T0) +N(u0) ≤ lim infT→T0

P (T ) + lim infu→u0

N(u).

Abbiamo percio dimostrato che il nostro funzionale F (T, u) e debolmente semi-continuo inferiormente in un punto (T0, u0) ∈ H1(Ω) × R(Γ1). Tuttavia, (T0, u0)rappresenta un punto qualsiasi dell’insieme H1(Ω)×R(Γ1), quindi F e debolmentesemicontinuo inferiormente in tutto H1(Ω)× R(Γ1).

5. Sia (Tk, uk) ∈ H1(Ω)×R(Γ1) una successione minimizzante tale che |P (Tk)| ≤Me G(Tk, uk) = 0, cioe:

|∫Γ3

(Tk − Td)2dΓ| ≤M

−∇ · (k∇Tk) = 0.

Allora possiamo dire che:

−∫Ω∇ · (k∇Tk)vdΩ = 0 ∀v ∈ H1(Ω),

cioe che:

−∫

Ω∇ · (k∇Tk)vdΩ =

∫Ωk∇Tk∇vdΩ +

∫Γ1

(h+ uk)TkvdΓ +

−∫

Γ1

(h+ uk)T∞vdΓ−∫

Γ3

qvdΓ = 0 .

v e una qualsiasi funzione appartenente allo spazioH1(Ω), quindi possiamo sceglierev = Tk e sostituirla nell’ultima equazione:∫

Ωk∇Tk∇TkdΩ +

∫Γ1

(h+ uk)TkTkdΓ +

−∫

Γ1

(h+ uk)T∞TkdΓ−∫

Γ3

qTkdΓ = 0 ,

che, portando a destra dell’uguale gli ultimi tre termini, risulta:∫Ωk(∇Tk)2dΩ

= −∫

Γ1

(h+ uk)(Tk)2dΓ +∫

Γ1

(h+ uk)T∞TkdΓ +∫

Γ3

qTkdΓ .

Analizziamo gli integrali dell’ultima equazione separatamente:

•∫Ω k(∇Tk)2dΩ = k||∇Tk||2L2(Ω).

19

•∫Γ1

(h+ uk)(Tk)2dΓ = (h+ uk)||Tk||2L2(Γ1),che, applicando la disuguaglianza di traccia con costante C1 > 0, diventa:∫Γ1

(h+ uk)(Tk)2dΓ = (h+ uk)||Tk||2L2(Γ1) ≤ C21 (h+ uk)||Tk||2H1(Ω).

• Applicando la disuguaglianza di Cauchy-Schwarz al terzo integrale, otteni-amo: ∫

Γ1

(h+ uk)T∞TkdΓ ≤ ||(h+ uk)T∞||L2(Γ1)||Tk||L2(Γ1)

= (h+ uk)T∞ |Γ1| ||Tk||L2(Γ1) ,

e, applicando la disuguaglianza di traccia con costante C2 > 0:∫Γ1

(h+ uk)T∞TkdΓ ≤ C2(h+ uk)T∞ |Γ1| ||Tk||H1(Ω).Poiche C2 e una costante arbitraria ma positiva, possiamo sceglierla in modotale da rendere tutto il secondo membro della disequazione maggiore di 1. Intal caso, possiamo scrivere l’ulteriore disuguaglinza:∫Γ1

(h+ uk)T∞TkdΓ ≤ C22 (h+ uk)2T 2

∞ |Γ1|2 ||Tk||2H1(Ω).

• Sappiamo per ipotesi che:|∫Γ3

(Tk − Td)2dΓ| ≤M

quindi, la funzione integranda (Tk − Td)2 e limitata e, in particolare, abbi-amo che ||Tk − Td||2L2(Γ3) e limitata, quindi anche la norma di ||Tk||2L2(Γ3) elimitata da una costante che chiamiamo M ′. Considerando l’ultimo integralee applicando la disuguaglianza di Cauchy-Schwarz, possiamo dire che:∫Γ3qTkdΓ ≤ ||q||2L2(Γ3)||Tk||2L2(Γ3) = q |Γ3| ||Tk||2L2(Γ3) ≤ q |Γ3|M ′.

Raggruppando tutti i termini trovati, otteniamo:

k||∇Tk||2L2(Ω) ≤ −C21 (h+uk)||Tk||2H1(Ω)+C

22 (h+uk)2T 2

∞ |Γ1|2 ||Tk||2H1(Ω)+q |Γ3|M ′,

che, scomponendo la norma H1(Ω), diventa:

k||∇Tk||2L2(Ω) + C21 (h+ uk)||Tk||2L2(Ω) + C2

1 (h+ uk)||∇Tk||2L2(Ω) +

−C22 (h+ uk)2T 2

∞ |Γ1|2 ||Tk||2L2(Ω) − C22 (h+ uk)2T 2

∞ |Γ1|2 ||∇Tk||2L2(Ω) ≤ q |Γ3|M ′ .

Raccogliendo le costanti al primo membro per tornare alla norma in H1(Ω):

maxk + C21 (h+ uk)− C2

2 (h+ uk)2T 2∞ |Γ1|2 , C2

1 (h+ uk)+

− C22 (h+ uk)2T 2

∞ |Γ1|2||Tk||2H1(Ω) ≤

≤ q |Γ3|M ′ ,

e chiamando:

C = maxk+C2(h+uk)−C2(h+uk)2T 2∞ |Γ1|2 , C2(h+uk)−C2(h+uk)2T 2

∞ |Γ1|2,otteniamo:

C||Tk||2H1(Ω) ≤ q |Γ3|M ′,

20

cioe:

||Tk||2H1(Ω) ≤ q |Γ3| M ′

C.

Ribattezzando M0 = q |Γ3| M ′

C, arriviamo alla dimostrazione dell’ultima ipotesi:

||Tk||2H1(Ω) ≤M0.

Nel problema di controllo da studiare, esiste quindi una soluzione (T, u) ∈ H1(Ω)×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 ditipo iterativo (si vedano [Q08], [D08]). Una volta definito l’algoritmo iterativo, sara poinecessario operare una scelta sul tipo di discretizzazione numerica da utilizzare per larisoluzione 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 principiodi minimo, le quali garantiscono che la soluzione trovata sia quella che minimizza ilfunzionale lagrangiano espresso dalla (21). Una volta definito quindi il problema dicontrollo, il metodo iterativo per la ricerca dell’ottimo puo essere schematizzato nelseguente 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 statoT ;

• noto T si calcola il valore del funzionale costo J (o del funzionale lagrangianoL: e possibile utilizzare indifferentemente J o L, dato che dalla definizione diquest’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 variabileaggiunta p;

• dato p, si calcola il valore della derivata del funzionale lagrangiano Lu (da qui inpoi indicata con J ′) tramite il principio di minimo;

• a questo punto, se il valore del funzionale J e della sua derivata J ′ (che indicanola distanza dall’ottimo, rispetto ad una tolleranza scelta) rispettano opportunicriteri di arresto, il ciclo termina e il valore di u all’iterazione corrente coincidecon l’ottimo; in caso contrario si aggiorna opportunamente il valore del controllou e si esegue un’altra iterazione, fino a quando non saranno soddisfatti i criteri diarresto del ciclo.

In particolare abbiamo scelto un criterio di arresto (che indica il raggiungimento dell’ot-timo) basato sulla valutazione della norma della derivata del funzionale costo rispetto

21

ad una tolleranza fissata. Il valore di tale tolleranza viene scelto sufficientemente piccoloin relazione al valore iniziale della quantita monitorata e al livello di vicinanza all’ottimoreale che si vuole ottenere. Il criterio di arresto per ogni iterazione k si traduce quindinella seguente: ∥∥J ′(uk)

∥∥‖J ′(u0)‖

≤ 10−6 . (33)

Per quanto riguarda l’aggiornamento della variabile di controllo u e possibile utilizzareun metodo di tipo gradiente, cioe un algoritmo iterativo che data la soluzione al passok calcola quella al passo successivo k + 1 nel modo seguente (si utilizza la notazionevettoriale in quanto piu generale):

uk+1 = uk + τkdk k = 0, 1, 2, ... , (34)

dove τk indica un parametro reale e strettamente positivo detto parametro di acceler-azione (o di rilassamento), mentre dk rappresenta la direzione di discesa, cioe tale percui vale:

dk T · ∇J(uk) < 0 se ∇J(uk) 6= 0.

A seconda delle possibili scelte per la direzione dk, si ottengono diversi metodi iteratividi 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 peril problema di controllo in esame abbiamo scelto di utilizzare gli ultimi due metodi so-praelencati.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 − τkJ ′(uk) k = 0, 1, 2, ... . (35)

La valutazione della derivata del funzionale J ′ in un dato punto del dominio di control-lo fornisce un’indicazione della sensitivita del funzionale J in quel punto rispetto allevariazioni della variabile di controllo u: questo fa sı che il metodo di discesa piu ripidasia, tra gli altri possibili, il metodo piu utile alla comprensione del ruolo giocato da J ′

e quindi dalla variabile aggiunta p.Il metodo del gradiente coniugato (CG, conjugate gradient) invece corrisponde allaseguente scelta per la direzione di discesa dk:

dk =

−∇J(uk) per k = 0,−∇J(uk) + γkdk−1 per k ≥ 1,

dove γk e uno scalare da scegliere in modo tale che dk T · dk−1 = 0. Questo parametropuo essere calcolato in uno dei possibili modi seguenti:

γkPRP =

∇J(uk)T(∇J(uk)−∇J(uk−1)

)∇J(uk−1)T ∇J(uk−1)

(Polak-Ribiere-Poylak),

γkPRP+ = max

γk

PRP , 0

(Polak-Ribiere-Poylak modificato),

γkCD = −∇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 traducein:

uk+1 = uk + τkdk k = 0, 1, 2, ... ,

con dk =

−J ′(uk) per k = 0,−J ′(uk) + γkdk−1 per k ≥ 1,

(36)

e γk puo essere calcolato in uno dei seguenti modi:

γkPRP =

J ′(uk)(J ′(uk)− J ′(uk−1)

)‖J ′(uk−1)‖2 , (37)

γkPRP+ = max

γk

PRP , 0, (38)

γkCD = −

∥∥J ′(uk)∥∥2

dk−1 J ′(uk−1). (39)

Entrambi i metodi presentati comportano la valutazione di un appropriato valore delparametro di accelerazione τk. La scelta di τk deve essere in grado di garantire laconvergenza del funzionale costo al minimo, deve dunque valere ad ogni passo del cicloiterativo 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 dasoddisfare ad ogni iterazione la condizione di discesa (40); per alcuni sempliciproblemi 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); inaltre parole ad ogni iterazione τk viene posto pari al valore fissato τ , si calcolanoil controllo uk+1 e il corrispondente valore del funzionale J(uk+1) e si verifica lacondizione di discesa: se la (40) e soddisfatta si procede con il ciclo iterativo, incaso contrario, ovvero se uk+1 fornisce una soluzione peggiore della precedente, siritorna 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 euna versione particolare e semplificata del metodo di Armijo, che rientra nellacategoria dei metodi line search;

3. nel caso in cui risulti noto il valore finito del funzionale costo in corrispondenzadell’ottimo, Jopt = infu∈U J(u) ≥ 0, si puo scegliere il parametro di accelerazionecome τk =

∣∣J(uk)− Jopt

∣∣ /∥∥∇J(uk)∥∥2; in questo caso pero non e piu garantita

23

la monotonia della successioneJ(uk)

, il che significa che puo accadere che la

condizione (40) non sia rispettata ad ogni iterazione; per questo motivo il valoreattribuito a Jopt deve essere il piu preciso possibile, in quanto valori anche pocodistanti da quello esatto possono rendere molto lenta la convergenza dell’algoritmoiterativo e in alcuni casi anche portare a divergenza .

La scelta del parametro τk, oltre a soddisfare la (40), la quale assicura la convergenzadel funzionale al minimo, deve anche garantire che la convergenza venga raggiunta il piurapidamente possibile. E possibile quindi introdurre criteri piu severi della (40) nellascelta dei possibili valori di τk, con il fine di evitare due diversi tipi di problemi: unavelocita di decrescita troppo bassa e l’uso di passi troppo piccoli. La prima difficolta sipuo ovviare richiedendo che:

J(uk)− J(uk+1) ≥ −στkdk T · ∇J(uk) , (41)

dove σ ∈ (0, 1/2). Questa condizione equivale a richiedere che la velocita media didecrescita lungo dk di J in uk+1 sia almeno pari ad una frazione assegnata della velocitadi decrescita iniziale in uk. La seconda difficolta invece puo essere superata richiedendoche la velocita di decrescita in uk+1 non sia inferiore ad una frazione assegnata di quellain uk, attraverso la seguente:∣∣∣dk T · ∇J(uk+1)

∣∣∣ ≤ η∣∣∣dk T · ∇J(uk)

∣∣∣ , (42)

dove η ∈ (σ, 1), cosı da garantire che valga contemporaneamente anche la (41). Inparticolare per il nostro problema di controllo le condizioni precedenti si scrivono come:

J(uk)− J(uk+1) ≥ −στkdkJ ′(uk)∣∣∣dkJ ′(uk+1)∣∣∣ ≤ η

∣∣∣dkJ ′(uk)∣∣∣ .

Nella pratica si scelgono σ ∈ [10−5, 10−1] e η ∈ [0.1, 0.5]. Dunque per il calcolo delparametro τk si scegliera uno dei metodi sopraelencati compatibilmente con il rispettodelle condizioni (40), (41) e (42).

24

4 Implementazione del codice

In questa sezione diamo una descrizione generale del codice da noi implementato, senzasoffermarci troppo sui dettagli, per i quali rimandiamo alla documentazione Doxygen.Sviluppiamo il codice per le simulazioni numeriche tramite il linguaggio di program-mazione C++. La risoluzione del problema in esame richiede l’implementazione delleseguenti 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’e-quazione 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 daincludere (librerie o altri file .hpp) e le dichiarazioni delle classi e/o delle funzioni, chevengono poi definite nel secondo file corrispondente con estensione .cpp. Questa suddi-visione 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, tracui vector, list, set e map. Inoltre per la risoluzione dei sistemi lineari si e scelto diutilizzare la libreria UMFPACK, che implementa un metodo di risoluzione diretto moltoefficace 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 utilizzatoGetPot: un parser che permette di processare un file di testo (chiamato datafile) checontiene tutte le informazioni utili per l’esecuzione del programma. L’enorme vantaggiofornito dall’impiego di un file di questo tipo per il passaggio dei parametri consiste nellapossibilita 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 unoggetto di classe GetPot, che verra poi passato ai costruttori delle classi che necessitanodi leggere parametri dall’esterno. Scegliamo di creare la mesh computazionale tramiteFreeFem++, salvandola in un file di estensione .msh. Il nostro programma leggera poitale file (tramite la libreria fstream), memorizzando la mesh in una opportuna strutturadati. 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 costi-tuiscono la griglia.

class Pointprivate:

double _x;double _y;int _id_node;int _id_boundary;

...

Ogni nodo viene memorizzato tramite: le coordinate x e y, il numero identificativoid 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 creareun 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 definiscequindi il singolo elemento della mesh e contiene: un vettore di tre puntatori ad oggettidi classe Punto, element (si e preferito creare un vettore di puntatori e non di oggettidi classe Punto perche, dal momento che nella struttura dati della mesh ci sara gia unvettore contenente tutti i nodi della griglia, sarebbe inutile memorizzarli una secondavolta) ; il numero identificativo di ogni elemento, id element. I metodi di questaclasse permettono di creare un certo triangolo e di ricavare o stampare a video tutte leinformazioni relative all’elemento stesso o ai nodi ad esso associati.

class Triangleprivate:

int _id_element;std::vector <Point*> _element;std::vector <double> _cof_jacob;

...

Dal momento che il calcolo degli integrali bidimensionali, necessari per la costruzionedei sistemi lineari associati al problema di stato e al problema aggiunto, verra effettuatopassando all’elemento di riferimento (il triangolo di vertici (0, 0), (0, 1) e (1, 0)), e neces-sario conoscere i coefficienti della matrice Jacobiana associata alla trasformazione lineareaffine tra elemento generico della griglia ed elemento di riferimento. Tali coefficienti sipossono 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 cal-colo 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 deisistemi lineari richiedono il calcolo di integrali monodimensionali sui lati di bordo. Perquesto motivo e necessario creare una classe che contenga tutte le informazioni relativeai lati di bordo.

class Edgeprivate:

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 delbordo a cui appartiene ( id bound), i numeri identificativi dei nodi che rappresentanol’estremo sinistro ( id sx) e l’estremo destro ( id dx), tutte le informazioni sui dueestremi 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’elementodi riferimento (il segmento di lunghezza unitaria), e necessario memorizzare la lunghezzadi ogni lato di bordo ( length).

Le tre classi appena descritte, Point, Triangle e Edge, vengono utilizzate dallaclasse MeshFF per memorizzare la griglia computazionale, leggendo tutte le informazioninecessarie dal file di estensione .msh.

class MeshFFprivate:

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 elementidi classe Triangle; i nodi della griglia sono salvati invece nel vettore nodes, che contieneelementi 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 lospazio per i tre vettori corrispondenti nodes, grid e edges. Inoltre la classe contienele 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 quellibidimensionali, 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 memorizzatoattraverso near, un puntatore ad una lista di interi, che verra istanziato in modo taleche per ogni nodo contenga la lista di tutti i suoi vicini, se stesso compreso. I vicini sonomemorizzati in ordine crescente. I metodi della classe Mesh permettono di costruire lastruttura dati contenente la griglia; consentono inoltre di trovare e salvare la lista deivicini per ogni nodo e di stampare a video tutte le altre informazioni relative alla meshcomputazionale.

Tutta la struttura dati per la mesh appena illustrata e pensata per essere costruitaa partire da un file di tipo .msh fornito da FreeFem++, che noi abbiamo scelto comegeneratore di griglia per il nostro problema. Il programma e comunque in grado disupportare l’utilizzo di una generica classe Mesh, diversa da MeshFF, in quanto tutto ilcodice e template su Mesh: e dunque possibile costruire e utilizzare una struttura datidiversa da quella presentata, che possa ad esempio leggere le informazioni a partire daun 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’implementazionedi due integratori, uno bidimensionale e uno monodimensionale. In entrambi i casi, ilcalcolo viene riportato sull’elemento di riferimento.

Gli integrali bidimensionali da calcolare sono associati ai termini diffusivi dell’e-quazione di stato e dell’equazione aggiunta e assumono entrambi la seguente forma:

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:

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 definitelocalmente su ogni singolo elemento. Per riportare il calcolo di ogni singolo integralesull’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 triangolodi vertici (0, 0), (0, 1) e (1, 0)), x indica il sistema di coordinate sull’elemento genericocorrente e la matrice J e lo Jacobiano della trasformazione affine. In particolare, lamatrice J e il vettore c sono definiti nel seguente modo:

J =[x2 − x1 x3 − x1

y2 − y1 y3 − y1

], c =

[x1

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 gradientidelle funzioni di base che, riportati sull’elemento di riferimento, assumono per il verticei 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’ele-mento di riferimento, l’integrale da calcolare (43) diventa:

k

∫Ω

(∇T

xφ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 funzionidi base sull’elemento di riferimento e creare una classe che implementi una formuladi integrazione numerica di tipo Gaussiano. I file base.hpp e base.cpp contengonorispettivamente le dichiarazioni e le definizioni dei prodotti scalari tra i gradienti dellefunzioni di base, necessari per il calcolo dell’espressione (44). La classe Quadraturacontiene 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 direzionex e y rispettivamente; weight contiene il valore del peso di quadratura; pfn indica unpuntatore ad una funzione di due variabili ed e definito nel file typedef.hpp: integrandcontiene la funzione da integrare, cioe una opportuna funzione definita nel file base.cpp.Il metodo integrate implementa la formula di quadratura gaussiana sommando il valoredella funzione integrand nei nodi di quadratura e moltiplicando il valore finale per ilpeso di quadratura.

Gli integrali monodimensionali da calcolare servono sia per la valutazione del fun-zionale costo e della sua derivata, sia per includere le condizioni al contorno dei sistemiassociati all’equazione di stato e all’equazione aggiunta. Le condizioni al bordo vanno inparte aggiunte alla matrice del sistema e in parte al termine noto. Tutti questi terminipossono essere ridotti alla valutazione di due diversi integrali:

∫Γψi(x)ψj(x)dΓ ,∫

Γψi(x)dΓ ,

dove le funzioni di base ψi sono la restrizione delle funzioni di base φi sui bordi deldominio di integrazione. Per passare all’elemento di riferimento (segmento di lughezzaunitaria: estremo sinistro x1 = 0 ed estremo destro x2 = 1) e necessario il cambio divariabili:

x = (x2 − x1)x+ a ,

dove x1 e x2 sono riferiti all’elemento corrente, x indica la coordinata dell’elementocorrente 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 baseassociata all’estremo destro. Il calcolo degli integrali monodimensionali richiede quindila valutazione dei seguenti termini:

(x2 − x1)∫

Γψ1(x)ψ2(x)dΓ ,

(x2 − x1)∫

Γψ1(x)ψ1(x)dΓ ,

(x2 − x1)∫

Γψ2(x)ψ2(x)dΓ ,

(x2 − x1)∫

Γψ1(x)dΓ ,

(x2 − x1)∫

Γψ2(x)dΓ ,

dove Γ indica l’elemento di riferimento monodimensionale. Grazie alla simmetria dellefunzioni di base e a quella dell’elemento di riferimento, il valore del secondo integralecoincide 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 definitanella 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 contieneil valore dei pesi di quadratura per ogni nodo; pfn1d indica un puntatore a una funzionedi una variabile ed e definito nel file typedef.hpp: integrand contiene la funzione daintegrare, cioe una opportuna funzione definita nel file base.cpp. Il metodo integrateimplementa la formula di quadratura sommando il valore della funzione integrand neinodi di quadratura e moltiplicando il valore finale per il peso di quadratura: la formuladi 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 prob-lema aggiunto, il calcolo del funzionale costo e il calcolo della sua derivata, richiedono l’u-

31

tilizzo di tutte le strutture definite precedentemente. E’ inoltre necessaria una appositastruttura dati dove poter memorizzare le matrici associate ai sistemi lineari.

I sistemi lineari vengono risolti tramite la libreria UMFPACK, che richiede in ingressouna matrice nel formato CSC (Compressed Sparse Column). Pertanto, tutte le matriciandranno salvate in tale formato. La struttura dati per le matrici in formato CSC consistein tre puntatori: matrix contiene gli elementi non nulli della matrice per colonna, rowscontiene gli indici di riga degli elementi di matrix e columns contiene nella i-esimacomponente 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 dirighe (uguale al numero di colonne dal momento che la matrice e quadrata) e nz indica ilnumero totale di elementi non nulli. Una volta creata la matrice, e possibile accedere aisingoli elementi e attribuirgli un valore diverso da zero tramite il metodo fill_Matriceoppure modificarli aggiungendo un certo valore a quello gia presente tramite il metodomodify_Matrice.

La classe Problem permette di calcolare tutte le incognite del problema tramitela risoluzione del sistema di stato e il sistema aggiunto e tramite la valutazione delfunzionale costo. Si tratta di una classe template su Mesh. Si noti inoltre che in questocaso tutti gli attributi sono dichiarati protected per implementare successivamentel’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 delsistema aggiunto e rimane inalterata ad ogni passo del ciclo iterativo dal momento chenon contiene la variabile di controllo u. La matrice _Astiffness viene percio creatauna sola volta per entrambi i problemi dal costruttore della classe Problem, chiamandola 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 essicalcola i nove integrali associati alle combinazioni dei tre vertici di ogni triangolo: gliintegrali bidimensionali devono essere valutati per tutti i possibili prodotti dei gradientidelle funzioni di base associate ad ogni nodo, come indicato dalla formula (44). Gliintegrali calcolati per ogni elemento vengono poi inseriti all’interno della matrice tramiteil metodo fill_matrice che riceve in ingresso gli indici di riga e colonna dell’elementoda riempire e il valore dell’integrale riportato sull’elemento corrente.

33

Tornando alla classe Problem, il metodo stato costruisce e risolve il sistema associatoall’equazione di stato. Tramite la funzione bc definita nel file assembla.hpp si aggiungealla matrice A_stiffness il termine incognito che deriva dalle condizioni al contorno:si tratta di calcolare degli integrali monodimensionali sul bordo Γ1 (tramite un ciclosui lati di bordo) e inserirli nella matrice tramite il metodo modify_Matrice. Tramitela funzione termine_noto_stato, sempre definita nel file assembla.hpp, si costruisceinvece 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 inmodo totalmente analogo al metodo stato. Tramite la funzione bc e la funzione

34

termine_noto_aggiunto, entrambe definite nel file assembla.hpp, il termine incog-nito derivante dalle condizioni al contorno viene aggiunto alla matrice A_stiffness eviene costruito il termine noto associato al sistema aggiunto.Sia il metodo stato che il metodo aggiunto utilizzano il metodo solve, che permettedi risolvere i sistemi lineri associati al problema di stato e al problema aggiunto tramitel’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’inte-gratore monodimensionale e la soluzione T del problema di stato _sol_stato. La partecontenente la variabile di controllo u non ha bisogno di integrazione dal momento chetale 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 necessarioeffettuare 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 costituitadalla 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 metodie tutti gli attributi di tale classe. In aggiunta a questi, la classe Ciclo contiene anchetutti 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 daldatafile, e che indicano il metodo iterativo che si intende utilizzare: _flag_metodopermette 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 delparametro di accelerazione τ (si vedano le scelte 2 e 3 nella sezione 3.2); _flag_gammapermette di scegliere tra i tre possibili modi di calcolare il parametro γ del metodo delgradiente coniugato ((37), (38), (39)); _flag_cond consente di scegliere quale delle trecondizioni (40), (41), (42) usare nel ciclo di aggiornamento di τ , nel caso si sia sceltal’opzione 2 per il calcolo di questo parametro. Il metodo inizializza a partire daun valore iniziale del controllo u0, calcola il valore iniziale della temperatura T , dellavariabile aggiunta p e il valore iniziale del funzionale costo e della sua derivata. Il meto-do itera compie l’iterazione: tramite _flag_metodo si sceglie se chiamare il metodoitera_SD (che implementa il metodo steepest descent) o itera_CG (che implementail metodo conjugate gradient). Il ciclo si arresta quando risulta essere soddisfatta lacondizione (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, chepermette di salvare i risultati ottenuti in opportuni file per la visualizzazione.Tramite la libreria fstream, la distribuzione finale della temperatura T e della variabileaggiunta p vengono scritti su due file; viene inoltre creato il file postproc.edp, checontiene le istruzioni necessarie per leggere i due file contenenti le variabili T e p, e chepermette la lora visualizzazione tramite FreeFem++; il metodo inoltre lancia direttamentel’esecuzione di FreeFem++ tramite il comando system.Il metodo postprocess consente inoltre di scrivere su file .dat i valori di alcune variabilidel problema, salvate in opportuni vettori durante l’esecuzione del ciclo iterativo. Inparticolare vengono salvati in questo modo i valori del funzionale costo J , della suaderivata J ′ e della variabile di controllo u. Vengono salvati inoltre su un altro file .dati valori della temperatura T corrispondenti ai nodi sul bordo Γ3, per poi visualizzarnel’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 simulazionieffettuate con il codice C++ illustrato nella sezione 4.

Consideriamo inizialmente il metodo di discesa piu ripida, con un valore iniziale dellavariabile 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 correntee a quella precedente. Il valore della tolleranza che permette l’uscita dal ciclo whileprincipale e fissato pari a 10−6, mentre il numero massimo di iterazioni e posto paria 100. Decidiamo inoltre di fissare il valore del parametro di coercivita β (si veda ladefinizione del funzionale costo (15)) pari a 10−3. Il valore ottimale del parametro dicontrollo u e del funzionale costo J , raggiunti dopo n iterazioni, sono i seguenti:

numero totale di iterazioni = 11valore ottimale di u = 38.2395valore finale del funzionale J = 0.773907

Il miglior funzionamento del dissipatore di calore per il chip e quindi garantito con uncoefficiente di convezione forzata dell’aria (h + u) pari a 43.94 W

m2C, ricordando che il

coefficiente di convezione naturale h nel nostro caso e pari a 5.7 Wm2C

.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 comeil valore di J decresca ad ogni iterazione e che il minimo e raggiunto in corrispondenzadell’ultima iterazione, associata appunto al valore ottimale di u (figura 2). La derivatadel funzionale costo J ′ assume un valore dell’ordine di −102 alla prima iterazione, cor-rispondente al caso con controllo nullo, ma dalla seconda iterazione in poi si mantienemolto piccolo fino a convergere a 0 (figura 3).La distribuzione della temperatura sul dissipatore (figura 4 in alto) presenta valori piuelevati in corrispondenza della fonte di calore (il chip) che si trova a contatto con la su-perficie in basso e valori che diminuiscono gradualmente allontanandosi da esso, grazieal lavoro di dispersione termica effettuato passivamente dal dissipatore e alla presenzadella ventola. La distribuzione della variabile aggiunta ha lo stesso andamento della tem-peratura, ma i suoi valori sono molto vicini allo zero (figura 4 in basso). La temperaturadesiderata di funzionamento del chip e pari a 65C: i valori di temperatura sul bordoΓ3, ovvero la superficie inferiore del dissipatore, che si trova a contatto con il chip, sonodi poco superiori a tale valore (figura 5). Il dissipatore e la ventola permettono quindidi mantenere la temperatura sul chip intorno al valore ideale desiderato, garantendo ilcorretto funzionamento del chip stesso. L’andamento della temperatura sul bordo Γ3

e influenzato in modo evidente dalla geometria del dissipatore: si nota infatti che latemperatura risulta piu bassa nella zona centrale del bordo, al di sopra della quale sitrova il tronco centrale del dissipatore dove avviene la convezione, mentre rimane piuelevata nelle regioni laterali, piu vicine alla parte di bordo adiabatica e dunque menoinfluenzata 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 ′ in funzione delle iterazioni.

39

Figura 4: Distribuzione della temperatura T e della variabile aggiunta p sul dissipatore incorrispondenza 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 delmetodo del gradiente coniugato. Per esempio, scegliendo il metodo PRP (37) per ilcalcolo del parametro γ e utilizzando l’opzione 2 per l’aggiornamento di τ abbinatoalla condizione (40) per il ciclo, si ottengono gli andamenti del funzionale J e della suaderivata J ′ mostrati in figura 6. In questo caso il numero di iterazioni per giungere aconvergenza 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 trediverse condizioni (40), (41) e (42) per il ciclo), oppure con calcolo di τ tramite Jott; ilmetodo conjugate gradient con le tre diverse scelte (PRP, PRP+ e CD) per il calcolodel parametro γ, abbinato alle diverse opzioni per il calcolo di τ . Come si puo notaredalla tabella, i valori di u e J in corrispondenza dell’ottimo che si ottengono sono tuttimolto simili. La principale differenza tra i diversi metodi consiste nel numero di iterazioninecessarie per raggiungere l’ottimo: tale valore rimane comunque molto limitato, essendocompreso tra un minimo di 9 ed un massimo di 25, dunque si puo affermare che tutti imetodi utilizzati hanno una velocita di convergenza confrontabile. Per quanto riguarda ilconfronto diretto tra il metodo steepest descent e il gradiente coniugato, come si osservadai risultati riportati nella tabella, non vi sono differenze significative, o perlomeno essenon si presentano nel particolare caso qui studiato. I risultati relativi alla scelta dellacondizione (42) per l’aggiornamento di τ non vengono riportati, dal momento che conquesta condizione non si raggiunge il valore ottimo: questo potrebbe essere dovuto al

41

Figura 6: Metodo CG: funzionale costo J e derivata J ′ in funzione delle iterazioni.

42

fatto che tale condizione risulta piu restrittiva delle altre due condizioni possibili e difatto ostacola l’avanzamento del ciclo. Occorre far notare inoltre, come gia accennatonella sezione 3.2, che il metodo per il calcolo di τ che utilizza il valore ottimale delfunzionale, Jott, funziona correttamente solo nel caso in cui viene fornito il valore esattopreciso di Jott: per valori che si discostano da quello esatto, anche solo dell’ordine del10%, il metodo impiega molte piu iterazioni per giungere a convergenza e in certi casipuo anche accadere che non converga affatto.

Metodo τ γ condizione n u JSD τ/2 / 1 11 38.2395 0.773907SD τ/2 / 2 10 38.2397 0.773907SD Jott / / 17 38.2229 0.773908CG τ/2 PRP 1 15 38.229 0.773907CG τ/2 PRP 2 9 38.2408 0.773907CG τ/2 PRP+ 1 17 38.2395 0.773907CG τ/2 PRP+ 2 16 38.2397 0.773907CG τ/2 CD 1 11 38.2368 0.773907CG τ/2 CD 2 9 38.2345 0.773907CG Jott PRP / 25 38.2206 0.773909CG Jott PRP+ / 17 38.2229 0.773908CG 2Jott CD / 12 38.2216 0.773908

Tabella 1: Tabella con i risultati: per ogni metodo e riportato il numero di iterazioni n necessarie perla convergenza, il controllo ottimo u trovato e il corrispondente valore del funzionale costo J . Le scelte1 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++, epossibile effettuare un confronto tra i risultati ottenuti con i due diversi programmi. Perogni 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 leggermenteda quelli gia calcolati. Tali differenze sono pero dell’ordine di 10−5 e possono quindi, inrelazione ai valori ottenuti, essere considerate trascurabili; anche la distribuzione dellatemperatura e della variabile aggiunta presentano solo differenze trascurabili.

Per i risultati illustrati finora, si e utilizzato un parametro di coercivita β pari a10−3: in questo modo si ottiene la coercivita del problema introducendo solo una piccolaperturbazione. 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 difficiletrovare il minimo: per esempio nel caso gia considerato con β = 10−3 bastano 11iterazioni per raggiungere il valore ottimo di u, mentre per β = 10−5 sono necessarieben 598 iterazioni. Questo comportamento e dovuto al fatto che al diminuire del valoredi β il funzionale costo J assume una forma sempre piu appiattita ed e quindi piu difficiletrovarne il minimo. In figura 7 e possibile notare il diverso andamento del funzionalecosto rispetto ad u, nel caso β = 10−2 e β = 10−3. Inoltre si osserva che utilizzandovalori di β sempre piu piccoli, il valore ottimale di u aumenta, mentre diminuiscono ilvalore del funzionale costo J e i valori di temperatura T , che si avvicinano sempre dipiu 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 ottimoad un caso reale: il mantenimento della temperatura su di un chip ad un certo livellodesiderato, in modo da evitare il malfunzionamento o la rottura del chip stesso. Perraggiungere 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’equazionealle derivate parziali (equazione di diffusione con condizioni al contorno di Neumann eRobin) a cui e stato associato un funzionale costo per tenere conto sia della temperaturadesiderata sul chip, sia del parametro di controllo. Dopo un’analisi matematica teori-ca del modello formulato, la sua discretizzazione e successiva implementazione tramiteFreeFem++ e C++, hanno permesso di trovare il valore ottimale del parametro di controlloche permette il corretto funzionamento del chip.

Per la risoluzione numerica del problema sono stati utilizzati diversi metodi di ot-timizzazione: diverse versioni del metodo di steepest descent e del metodo del gradi-ente coniugato. Gli schemi implementati hanno fornito risultati confrontabili e in buonaccordo con la fisica del problema.

Inoltre, questo lavoro prevedeva l’apprendimento del linguaggio di programmazioneC++, conseguito tramite lo sviluppo del codice a elementi finiti necessario per risolvere ilproblema fornito dalla teoria del controllo ottimo. L’analisi approfondita di un problemacompleto, partendo dall’analisi teorica fino ad arrivare all’implementazione del codiceper le simulazioni, ci ha permesso di comprendere tutte le fasi in cui viene affrontato unproblema 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; //conducibilita termica alluminio [W/mC]real h=5.7; //convezione naturale [W/m^2C]real q=15.08e3; //flusso di calore [W/m^2]

real winf=25.0; //temperatura ambiente [C]real wd=65.0; //temperatura desiderata [C]real b=0.001; //coefficiente termine coercivo [C^4m^4/W^2]func f=0.0; //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; //controllo iniziale [W/m^2C]u0=u;

real taufix=500.0; //parametro di accelerazione iniziale

real tau=taufix;

real sigma=1.e-3; //parametro per seconda condizione

real eta=0.5; //parametro per terza condizione

real tol=1.e-6; //tolleranza

int n=0; //contatore algoritmo iterativo

int nmax=100; //numero massimo iterazioni

real[int] Jplot(nmax+1); //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; //conducibilita termica alluminio [W/mC]real h=5.7; //convezione naturale [W/m^2C]real q=15.08e3; //flusso di calore [W/m^2]

real winf=25.0; //temperatura ambiente [C]real wd=65.0; //temperatura desiderata [C]real b=0.001; //coefficiente termine coercivo [C^4m^4/W^2]func f=0.0; //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; //controllo iniziale [W/m^2C]u0=u;

real tol=1.e-6; //tolleranza

int n=0; //contatore algoritmo iterativo

int nmax=100; //numero massimo iterazioni

real[int] Jplot(nmax+1); //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<<" 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,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; //conducibilita termica alluminio [W/mC]real h=5.7; //convezione naturale [W/m^2C]real q=15.08e3; //flusso di calore [W/m^2]

real winf=25.0; //temperatura ambiente [C]real wd=65.0; //temperatura desiderata [C]real b=0.001; //coefficiente termine coercivo [C^4m^4/W^2]func f=0.0; //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; //controllo iniziale [W/m^2C]u0=u;

real taufix=500.0; //parametro di accelerazione iniziale

real tau=taufix;

real sigma=1.e-3; //parametro per seconda condizione

real eta=0.3; //parametro per terza condizione

real gamma; //parametro gradiente coniugato

real tol=1.e-6; //tolleranza

int n=0; //contatore algoritmo iterativo

int nmax=100; //numero massimo iterazioni

real[int] Jplot(nmax+1); //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; //conducibilita termica alluminio [W/mC]real h=5.7; //convezione naturale [W/m^2C]real q=15.08e3; //flusso di calore [W/m^2]

real winf=25.0; //temperatura ambiente [C]real wd=65.0; //temperatura desiderata [C]real b=0.001; //coefficiente termine coercivo [C^4m^4/W^2]func f=0.0; //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; //controllo iniziale [W/m^2C]u0=u;

real gamma; //parametro gradiente coniugato

real tol=1.e-6; //tolleranza

int n=0; //contatore algoritmo iterativo

int nmax=100; //numero massimo iterazioni

real[int] Jplot(nmax+1); //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<<" 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;

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

[C98] Y. A. Cengel, Termodinamica e Trasmissione del Calore, McGraw-Hill, Milano,1998.

[D08] L. Dede, Adaptive and Reduced Basis Methods for Optimal Control Problems inEnvironmental Applications, Tesi di dottorato, Politecnico di Milano, 2008.

[EG04] A. Ern, J.-L. Guermond, Theory and Practice of Finite Elements, Springer-Verlag, New York, 2004.

[Q08] A. Quarteroni, Modellistica Numerica per Problemi Differenziali, Springer-VerlagItalia, Milano, 2008.

[QSS00] A. Quarteroni, R. Sacco, F. Saleri, Matematica Numerica, Springer-VerlagItalia, Milano, 2000.

[S07] S. Salsa, Elementi di Ottimizzazione e Teoria del Controllo, Dispense del corso diEDP2, Politecnico di Milano, 2007.

[Y01] D. Yang, C++ and Object Oriented Numeric Computing for Scientists aneEngineers, Springer-Verlag, New York, 2001.

[C++] C++.com, web page, http://www.cplusplus.com.

[FreeFem++] FreeFem++.org, web page, http://freefem.org.

59