Applicazione della Proper Orthogonal Decomposition (POD) al Problema di Minima

Compliance

Pierluigi Colombo

Abstract In molte applicazioni, le classiche procedure di discretizzazione possono portare ad alti

costi computazionali, a causa del grande numero di gradi di libertà che la

discretizzazione comporta. Un modo efficiente di ridurre la complessità del problema

discretizzato è la Proper Orthogonal Decomposition (POD). Questo metodo consiste

nel catturare la dinamica dominante del fenomeno fisico attraverso delle soluzioni già

calcolate, e la costruzione a partire da queste di un sottospazio di dimensione ridotta

con una base ortonormale. Come vedremo questa tecnica è in grado di ridurre in

maniera considerevole i costi computazionali del calcolo di una nuova soluzione. In

questo lavoro si tratterà della POD, del problema di ottimizzazione topologica di

minima compliance e di come la POD vi è stata implementata. L’algoritmo di

ottimizzazione topologica unito alla POD è stato poi applicato ad un problema concreto

ed i risultati sono stai poi comparati con quelli dell’algoritmo di ottimizzazione

topologica classico.

1 POD e relazione con SVD In questo capitolo si presenta l’obiettivo della POD e la sua relazione con la Singular

Value Decomposition (SVD). L’obiettivo è di trovare una base ortonormale dello

spazio generato da 𝑛 vettori 𝑦1, 𝑦2, …, 𝑦𝑛 con 𝑦𝑧 ∈ ℝ𝑚 con 𝑧 = 1, … , 𝑛.

Sia 𝑌 ∈ ℝ𝑚×𝑛 la matrice le cui colonne sono i vettori 𝑦1, 𝑦2,…, 𝑦𝑛. La SVD ci

garantisce che esistono due matrici ortogonali 𝑈 ∈ ℝ𝑚×𝑚 e 𝑉 ∈ ℝ𝑛×𝑛 ed una

matrice diagonale Σ ∈ ℝ𝑚×𝑛 tale che:

𝑌 = 𝑈Σ𝑉𝑇 (1)

In particolare esistono 𝑟 valori positivi 𝜎1 ≥ 𝜎2 ≥ ⋯ ≥ 𝜎𝑟 tale che:

Σ = (𝐷 00𝑇 𝑂

) (2)

Con 𝐷 = 𝑑𝑖𝑎𝑔(𝜎1, 𝜎2, … , 𝜎𝑟). Le colonne {ϕ𝑗}𝑗=1

𝑚 di 𝑈 e {ψ𝑖}𝑖=1

𝑛 di 𝑉 sono dette

rispettivamente vettori singolari sinistri e vettori singolari destri di 𝑌, ed i valori

{𝜎𝑘}𝑗=1𝑟 sono detti valori singolari di 𝑌

Per le matrici 𝑌𝑌𝑇 e 𝑌𝑇𝑌, simmetriche e definite positive, valgono le seguenti

relazioni:

𝑌𝑌𝑇 = (𝑈Σ𝑉𝑇)(𝑈Σ𝑉𝑇)𝑇 = 𝑈ΣΣT𝑈𝑇 (3)

𝑌𝑇𝑌 = (𝑈Σ𝑉𝑇)𝑇(𝑈Σ𝑉𝑇) = 𝑉ΣTΣ𝑉𝑇 (4)

quindi {(ϕ𝑗 , 𝜎𝑗2)}

𝑗=1

𝑟 e {(ψ𝑖 , 𝜎𝑖

2)}𝑖=1𝑟 sono coppie autovalore-autovettori per le matrici

𝑌𝑌𝑇 E 𝑌𝑇𝑌 rispettivamente. Le rimanti coppie {(ϕ𝑗 , 0)}𝑗=𝑟+1

𝑚(se 𝑚 > 𝑛) o

{(ψ𝑖 , 0)}𝑖=𝑟+1𝑛 (se 𝑛 > 𝑚) identificano una base per 𝐾𝑒𝑟(𝑌𝑌𝑇), 𝐾𝑒𝑟(𝑌𝑇𝑌).

Essendo in generale 𝑟 < min (𝑚, 𝑛), la SVD di 𝑌 può essere scritta in termini di matrici

ridotte:

𝑌 = 𝑈𝑟D𝑉𝑟𝑇 (5)

Con

𝑈𝑟 = [(ϕ1, ϕ2, … , 𝜙𝑟] ∈ ℝ𝑚×𝑟 (6)

𝑉𝑟 = [(ψ1, ψ2, … , ψ𝑟] ∈ ℝ𝑛×𝑟

(7)

Data una matrice 𝐴 ∈ ℝ𝑛×𝑛 simmetrica e positiva (𝐴 = 𝐴𝑇), e sia (𝜆, 𝑣𝜆) una coppia

autovalore-autovettore, Allora 𝐴𝑣𝜆 = 𝜆𝑣𝜆. Quindi moltiplicando per 𝐴𝑇 si ottiene

𝐴𝑇𝐴𝑣𝜆 = 𝜆𝐴𝑇𝑣𝜆 = 𝜆2𝑣𝜆. Quindi 𝜆2 è un autovalore della matrice 𝐴𝑇𝐴, che è a sua

volta il quadrato del valore singolare 𝜎𝜆 = 𝜆. Essendo per ipotesi la matrice singolare

𝐴 positiva si ottiene che gli autovalori di 𝐴 coincidono con suoi valori singolari.

Un’altra proprietà fondamentale della SVD è che le prime 𝑟 colonne della matrice

ortogonale 𝑈 sono una base per lo spazio delle colonne di 𝑌. Quindi le colonne

𝑦1, 𝑦2, … 𝑦𝑛 ammettono una espansione di Fourier (finita) rispetto agli autovalori

singolari {𝜙𝑗}𝑗=1

𝑟. Infatti sapendo che (𝑈𝑟)𝑇𝑈𝑟 = 𝐼𝑟 (𝑈𝑟 è ortogonale), si ha che:

𝑦𝑗 = ∑ (𝐷(𝑉𝑟)𝑇)𝑘𝑗𝑈∙,𝑘 =𝑟𝑘=1

∑ ((𝑈𝑟)𝑇𝑈𝑟𝐷(𝑉𝑟)𝑇)𝑘𝑗𝑈∙,𝑘 =𝑟𝑘=1 ∑ ((𝑈𝑟)𝑇𝑈𝑟𝐷(𝑉𝑟)𝑇)𝑘𝑗𝑈∙,𝑘 =𝑟

𝑘=1

∑ ((𝑈𝑟)𝑇𝑌)𝑘𝑗𝑈∙,𝑘 =𝑟𝑘=1 ∑ (𝑦𝑗 , 𝜙𝑘)

ℝ𝑚𝜙𝑘𝑟𝑘=1

Quindi

𝑦𝑗 = ∑(𝑦𝑗 , 𝜙𝑘)ℝ𝑚𝜙𝑘

𝑟

𝑘=1

(7)

Formuliamo ora il problema della POD. Si vuol trovare la base {𝜙𝑗}𝑗=1

𝑟 tale che l’errore

di approssimazione di ogni vettore 𝑦𝑖 ∈ ℝ𝑚 con la sua proiezione sui vettori

{𝜙𝑗}𝑗=1

𝑟sia minimizzato. Se consideriamo il caso di 𝑟 = 1 il problema è formulato nel

seguente modo:

𝑚𝑎𝑥 {∑|(𝑦𝑖 , �̃�)|2

, ϕ̃ ∈ ℝ𝑚: ‖�̃�‖ℝ𝑚

2𝑛

𝑖=1

= 1}

(P1)

𝑃1 è un problema di ottimizzazione vincolato che può essere risolto con il metodo dei

moltiplicatori di Lagrange. Qui (∙,∙) identifica il prodotto scalare in ℝ𝑚 e ‖∙‖ℝ𝑚 indica

la norma euclidea in ℝ𝑚. Introducendo la funzione 𝑒: ℝ𝑚 → ℝ definita

𝑒(𝜙) = 1 − ‖𝜙‖2, 𝜙 ∈ ℝ𝑚 (8)

il vincolo del problema 𝑃1 può essere espresso come 𝑒(𝜙) = 0. Da notare come

∇𝑒(𝜙) = −2𝜙, che è diverso da zero per ogni 𝜙 tale che 𝑒(𝜙) = 0, quindi la nostra

soluzione ottima è un punto regolare per il vincolo del problema 𝑃1. Sia ora

𝐿: ℝ𝑚 × ℝ → ℝ la lagrangiana associata al problema 𝑃1:

𝐿(𝜙, 𝜈) = ∑|(𝑦𝑖 , 𝜙)|2 +

𝑛

𝑖=1

𝜈(1 − ‖𝜙‖2)

(9)

Sia 𝜙 un ottimo locale del problema 𝑃1. Essendo 𝜙 un punto regolare per il vincolo

𝑒(𝜙) = 0 allora esiste un unico moltiplicare di Lagrange 𝜈 tale che ∇𝐿(𝜙, 𝜈) = 0.

Calcolando il gradiente di 𝐿 rispetto a 𝜙 si ottiene:

𝜕𝐿

𝜕𝜙𝑖

(𝜙, 𝜈) =𝜕

𝜕𝜙𝑖(∑ |∑ 𝑌𝑘𝑖𝜙𝑘

𝑚

𝑘=1

|

2

+ 𝜈(1 − ∑ 𝜙𝑘2)

𝑚

𝑘=1

𝑛

𝑗=1

) =

= 2 ∑ (∑ 𝑌𝑘𝑗𝜙𝑘𝑌𝑖𝑗

𝑚

𝑘=1

)

𝑛

𝑗=1

− 2𝜈𝜙𝑖 = 2(𝑌𝑌𝑇𝜙)𝑖 − 2𝜈𝜙𝑖

Quindi la soluzione 𝜙 del problema 𝑃1 risolve il problema agli autovalori

𝑌𝑌𝑇𝜙 = 𝜈𝜙

(10)

In più il valore della funzione obiettivo valutata in 𝜙 è uguale all’autovalore associato

in equazione 10. Prendendo il massimo degli autovalori della matrice 𝑌𝑌𝑇(che è

semidefinita positiva) si può dimostrare che questo è anche il massimo assunto dalla

funzione obiettivo del problema 𝑃1 e tale autovalore coincide con il valore singolare

più alto della matrice 𝑌. La soluzione quindi del problema 𝑃1 è il vettore singolare

sinistro associato al valore singolare più grande.

Questi ragionamenti possono essere estesi nel caso in cui non si cerchi un solo vettore

�̃� che massimizzi il problema 𝑃1 ma 𝑙 vettori 𝜙1, 𝜙2, … , 𝜙𝑙 tale che massimizzino il

problema 𝑃2 qui esposto:

𝑚𝑎𝑥 {∑ ∑|(𝑦𝑖 , 𝜙𝑘)|2

𝑛

𝑖=1

𝑙

𝑘=1

, (𝜙𝑝, 𝜙𝑞) = 𝛿𝑝𝑞 ∀𝑝, 𝑞: 1 ≤ 𝑝, 𝑞 ≤ 𝑙}

(P2)

La soluzione in questo caso è composta dai vettori che costituiscono le prime 𝑙 colonne

della matrice dei vettori singolari sinistri. In tale soluzione il valore del funzionale è

∑ 𝜎𝑖𝑙𝑖=1 . Un valore estremamente importante è la somma dei valori associati ai vettori

singolari sinistri che non sono usati nella creazione del sottospazio. Infatti tale valori

equivale all’errore commesso proiettando i vettori {𝑦𝑖}𝑖=1𝑛 nel sottospazio generato dai

vettori {𝜙𝑘}𝑘=1𝑙 . Quindi la POD consiste nel risolvere il problema agli autovalori della

matrice 𝑌𝑌𝑇 e troncare da un certo punto in poi i valori singolari della matrice 𝑌.

In genere però in molte applicazioni può succedere che il numero di gradi di liberta 𝑚

dei vettori 𝑦𝑖 sia di molto inferiore (diversi ordini di grandezza) del numero 𝑛 di vettori

𝑦𝑖 disponibili. In tal caso la matrice 𝑌𝑌𝑇 ∈ ℝ𝑚×𝑚 può essere estremamente grande.

Per poter risolvere il problema della POD però può essere sfruttata la matrice

𝑌𝑇𝑌 ∈ ℝ𝑛×𝑛, molto più piccola della precedente. L’idea è calcolare gli autovettori e

gli autovalori della matrice 𝑌𝑇𝑌. Infatti se 𝜓𝑖 è l’autovettore associato all’autovalore

𝜆𝑖, quindi

𝑌𝑇𝑌𝜓𝑖 = 𝜆𝑖𝜓𝑖

Allora

𝜙𝑖 =1

√𝜆𝑖

𝑌𝜓𝑖 , ∀ 𝑖: 1 ≤ 𝑖 ≤ 𝑙

Infatti se 𝜓𝑖 è un autovettore della matrice 𝑌𝑇𝑌 allora:

𝑌𝑇𝑌𝜓𝑖 = 𝜆𝑖𝜓𝑖 → (𝑌𝑌𝑇)𝑌𝜓𝑖 = 𝜆𝑖𝑌𝜓𝑖

Quindi 𝑌𝜓𝑖è autovettore della matrice 𝑌𝑌𝑇. Il fattore 1

√𝜆𝑖 è una normalizzazione.

2 Applicazione

La POD è stata sfruttata per risolvere un problema minima compliance, un problema

di ottimizzazione topologica, che ha l’obiettivo di costruire la più rigida struttura

utilizzando una quantità limitata di materiale. Per risolvere il problema di

ottimizzazione è stato utilizzato il S.I.M.P (Solid Isotropic Material with Penalization),

un metodo che introduce una appropriata funzione di densità per descrivere la

soluzione ottima.

In particolare, dato un corpo Ω ⊂ ℝdim (𝑑𝑖𝑚 = 2,3) con una forza 𝑓 applicata su Γ𝑁 ⊂

𝜕Ω, l’obiettivo è calcolare il sottoinsieme Ω𝑚𝑎𝑡 ⊂ Ω che minimizza il lavoro delle

forze esterne ∫ 𝑓 ∙ �⃗⃗�ΓN

𝑑𝛾 con �⃗⃗� il campo di spostamento.

2.1 Problema elastico

Qui è presentato il problema elastico per un materiale con coefficiente di Poisson 𝜇 e

modulo di Young 𝐸. Sia Ω un dominio, Γ𝐷 ⊂ 𝜕Ω la porzione del bordo in cui è fissato

uno spostamento, Γ𝑁 ⊂ 𝜕Ω la porzione del bordo in cui è imposto il carico 𝑓 e Γ𝐹 =

𝜕Ω\(Γ𝐷 ∪ Γ𝑁). Il problema dell’equilibrio elasto-statico consiste nel trovare il campo

di spostamento �⃗⃗� tale che:

{

−∇∙σ(�⃗⃗⃗�)=0 𝑖𝑛 Ω𝑢=0 𝑖𝑛 Γ𝐷

σ(�⃗⃗⃗�)�⃗⃗�=𝑓 𝑖𝑛 Γ𝑁

σ(�⃗⃗⃗�)�⃗⃗�=0 𝑖𝑛 Γ𝐹

(11)

Con �⃗⃗� il vettore unitario uscente ortogonale a 𝜕Ω e 𝜎 il tensore di stress. Con l’ipotesi

di piccoli spostamenti il tensore 𝜎 nella teoria elastica è:

𝜎(�⃗⃗�) = 2𝜇𝜖(�⃗⃗�) + 𝜆 𝑡𝑟(𝜖(�⃗⃗�))𝐼 (12)

con i coefficienti di Lamè definiti come:

𝜆 =𝐸𝜈

(1+𝜈)(1−𝜈) e 𝜇 =

𝐸

2(1+𝜈)

𝜖(�⃗⃗�) definito come:

𝜖(�⃗⃗�) =∇�⃗⃗� + ∇�⃗⃗�𝑇

2

ed 𝐼 il tensore identità. In accordo con la legge di Hook lo stress il tensore di stress può

essere scritto come

𝜎(�⃗⃗�) = 𝐸0𝜖(�⃗⃗�) (13)

𝐸0 è detto tensore di elasticità e dipende dalle proprietà del materiale 𝜆 e 𝜇.

2.2 Ottimizzazione topologica

Per inserire nel problema 11 le informazioni della struttura del dominio, oltre a quelle

del materiale, si definisce il tensore

𝐸𝑚𝑎𝑡 = 𝜒𝑚𝑎𝑡(𝑥)𝐸0 (14)

con 𝜒𝑚𝑎𝑡(𝑥) la funzione caratteristica del sottoinsieme Ω𝑚𝑎𝑡 ⊂ Ω definita:

𝜒𝑚𝑎𝑡(𝑥) = {1 𝑥 ∈ Ω𝑚𝑎𝑡

0 𝑥 ∉ Ω𝑚𝑎𝑡 (15)

in questo caso ogni dominio ottimo Ω𝑚𝑎𝑡sarà identificato univocamente dalla sua

funzione caratteristica 𝜒𝑚𝑎𝑡. Con l’introduzione del tensore 𝐸𝑚𝑎𝑡 è possibile ridefinire

il tensore 𝜎 come:

𝜎𝑚𝑎𝑡(�⃗⃗�) = 𝐸𝑚𝑎𝑡𝜖(�⃗⃗�) (16)

La formulazione debole quindi del problema (11), con all’interno anche la struttura del

dominio definita da (15) è: trovare �⃗⃗� ∈ 𝑈 tale che

𝑎(�⃗⃗�, �⃗�) = 𝑙(�⃗�) ∀�⃗� ∈ 𝑈 (17)

con

𝑎(�⃗⃗�, �⃗�) = ∫ 𝜎𝑚𝑎𝑡(�⃗⃗�): 𝜖(�⃗�)𝑑ΩΩ

(18)

𝑙(�⃗⃗�) = ∫ 𝑓 ∙ �⃗⃗�Γ𝑁

𝑑𝛾 (19)

con 𝑓definita come in (11), �⃗� funzione test e

𝑈 = {�⃗⃗� ∈ [𝐻1(Ω)]𝑑𝑖𝑚, 𝑑 = 2,3 ∶ �⃗⃗� = 0 𝑠𝑢 Γ𝐷} (20)

Il problema di minima compliance è quindi definito: trovare Ω𝑚𝑎𝑡 ⊂ Ω tale che

{

min(𝑙(�⃗⃗�)) 𝑐𝑜𝑛 �⃗⃗�(𝐸𝑚𝑎𝑡) ∈ 𝑈

𝑎(�⃗⃗�, �⃗�) = 𝑙(�⃗�) ∀�⃗� ∈ 𝑈

∫ 𝜒𝑚𝑎𝑡𝑑Ω = 𝑉𝑜𝑙(Ω𝑚𝑎𝑡) ≤ 𝛼𝑉0Ω

(21)

𝛼 è detto volume fraction ed è compreso nell’intervallo (0,1), 𝑉0 è il volume della

struttura iniziale. In genere la funzione binaria 𝜒𝑚𝑎𝑡 viene sostituita da una funzione

𝜌: Ω → [0,1]. Il tensore elastico quindi diventa

𝐸(𝑥) = ℎ(𝜌(𝑥))𝐸0 (22)

𝜌 è considerata una densità, nel senso che ∫ 𝜌𝑑ΩΩ

è uguale al volume della struttura.

Nel modello S.I.M.P si ha

ℎ(𝜌(𝑥)) = 𝜌𝑝 (23)

con 𝑝 > {2

1−𝜈,

4

1+𝜈}. In questo modo valori di 𝑝 ≠ {0,1} vengono penalizzati dal

modello poiché contribuiscono all’integrale ∫ 𝜌𝑑ΩΩ

non fornendo però rigidezza alla

struttura. Il problema (21) quindi diventa: trova 𝜌: Ω → [0,1] ∈ 𝐿∞ tale che

{

min(𝑙(�⃗⃗�)) 𝑐𝑜𝑛 �⃗⃗�(𝜌) ∈ 𝑈

𝑎𝜌(�⃗⃗�, �⃗�) = 𝑙(�⃗�) ∀�⃗� ∈ 𝑈

∫ 𝜌𝑑Ω = 𝑉𝑜𝑙(Ω𝑚𝑎𝑡) ≤ 𝛼𝑉0Ω 𝑐𝑜𝑛 𝜌min ≤ 𝜌 ≤ 𝜌

(24)

2.3 Algoritmo di ottimizzazione

Un approccio per risolvere il problema di ottimizzazione è:

Imposta una funzione di densità iniziale

While (criterio di converegenza non verificato)

Calcola lo spostamento con la densità in memoria;

Aggiorna la densità con un algoritmo di ottimizzazione;

end

L’ottimizzazione è stata svolta con la libreria NLopt (Versione 2.4.2). La proiezione

nel sottospazio costruito dalla POD viene utilizzata nel calcolo dello spostamento con

l’attuale funzione di densità. Definendo Ωℎ = {𝐾𝑗} una partizione del dominio Ω, si

definisce lo spazio 𝑋ℎ1 come:

𝑋ℎ1 = {𝑣ℎ ∈ 𝐶0(Ω̅): 𝑣ℎ|𝐾𝑗

∈ 𝑃1(𝐾𝑗) ∀ 𝐾𝑗 ∈ Ωℎ} (25)

e lo spazio 𝑉ℎ ⊂ 𝑉 lo spazio finito dimensionale degli spostamenti ammissibili come:

𝑉ℎ = {𝑣ℎ ∈ [𝑋ℎ1]𝑑𝑖𝑚: 𝑣ℎ = 0 𝑠𝑢 Γ𝐷} (26)

lo spostamento numerico all’iterazione 𝑖 del ciclo �⃗⃗�𝑛𝑢𝑚,𝑖 è soluzione del problema

𝐴�⃗⃗�𝑛𝑢𝑚,𝑖 = �⃗� (27)

con 𝐴 matrice il cui generico elemento (𝑖, 𝑗)è definito come

𝐴𝑖𝑗 = 𝑎(𝜙𝑖⃗⃗ ⃗⃗ , �⃗⃗�𝑗) (28)

e �⃗� vettore il cui generico elemento 𝑖 è definito

�⃗�𝑖 = 𝑙(�⃗⃗�𝑗) (29)

con �⃗⃗�𝑗 elemento della base 𝑉ℎ. Se la dimensione di 𝑋ℎ1 è 𝑁𝑒, la matrice

𝐴 ∈ ℝ𝑑𝑖𝑚∙𝑁𝑒×𝑑𝑖𝑚∙𝑁𝑒 ed �⃗� ∈ ℝ𝑑𝑖𝑚∙𝑁𝑒. Nella iterazione dell’algoritmo di ottimizzazione

non viene risolto il sistema (27), ma viene risolto una opportuna proiezione del sistema

(27) attraverso una matrice 𝑃 ∈ ℝ(𝑑𝑖𝑚∙𝑁𝑒)×𝑟, con 𝑟 ≪ 𝑑𝑖𝑚 ∙ 𝑁𝑒 ricavata dalla

decomposizione di Tucker in equazione (17). In questo caso viene risolto il sistema

𝑃𝑇𝐴𝑃�⃗⃗�𝑟𝑖𝑑,𝑖 = 𝑃𝑇�⃗� (30)

La matrice di cui viene svolta la SVD è la matrice 𝑼 i cui vettori sono gli spostamenti

ottimi 𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

della configurazione con volume fraction 𝛼𝑖. 𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

è il

campo di spostamento ottenuto con la densità ottima.

3 Implementazione

3.1 Librerie utilizzate e classe ElasticProblemSolver

Le librerie utilizzate sono 3:

• fenics/2016.1.0 per la discretizzazione del problema 17

• Eigen/3.2.8 per il calcolo dello spostamento 𝑢𝑛𝑢𝑚,𝑖 ad ogni iterazione di

ottimizzazione e per il calcolo della SVD (necessaria per la POD).

• NLopt/2.4.2 per l’ottimizzazione.

• Paraview versione 5.4.1 per la visualizzazione dei risultati.

A queste librerie si aggiunge una classe ElasticProblemSolver, creata da Leonardo

Locatelli e Francesco Clerici, e da me modificata per utilizzare le matrici sparse ed il

solver di Eigen e la POD. La classe contiene gli attributi necessari per la

discretizzazione del problema 28:

• Costanti fisiche: 𝜇 e 𝜆

• Vincoli: condizioni al bordo di Dirichlet e volume fraction imposta

• Discretizzazioni: spazi funzionali e mesh e vettori singolari sinistri per la POD

• Soluzioni: �⃗⃗�𝑛𝑢𝑚 e 𝜌𝑛𝑢𝑚

I principali metodi della classe sono:

• Un solver ElasticProblemSolver::Solve(). Il solver calcola lo spostamento �⃗⃗�𝑛𝑢𝑚

senza la POD usando:

1. Costanti fisiche: 𝜇 e 𝜆

2. Condizione al bordo di Dirichlet

3. Gli spazi funzionali

4. La densità 𝜌𝑛𝑢𝑚

Per il calcolo della soluzione �⃗⃗�𝑛𝑢𝑚 viene usato il solver di Fenics.

• Un metodo ElasticProblemSolver::get_compl() per il calcolo della compliance

usando:

1. Costanti fisiche: 𝜇 e 𝜆

2. Le soluzioni numeriche del problema 28 𝜌𝑛𝑢𝑚 e �⃗⃗�𝑛𝑢𝑚

• Un solver ElasticProblemSolver::Solve_eigen(). Il solver calcola lo spostamento

�⃗⃗�𝑛𝑢𝑚 senza la POD usando:

5. Costanti fisiche: 𝜇 e 𝜆

6. Condizione al bordo di Dirichlet

7. Gli spazi funzionali

8. La densità 𝜌𝑛𝑢𝑚

Per il calcolo della soluzione �⃗⃗�𝑛𝑢𝑚 viene usato il solver di Eigen. Il metodo

utilizza la funzione solver_eigen_sparso del file funzioni.cpp

• ElasticProblemSolver::Solve_proiettato(). Il solver calcola lo spostamento

�⃗⃗�𝑛𝑢𝑚 con la POD usando:

1. Costanti fisiche: 𝜇 e 𝜆

2. Condizione al bordo di Dirichlet

3. Gli spazi funzionali

4. La densità 𝜌𝑛𝑢𝑚

5. I vettori singolari sinistri.

Per il calcolo della soluzione �⃗⃗�𝑛𝑢𝑚 viene usato il solver di Eigen. Il metodo utilizza la

funzione solver_eigen_proiettato del file funzioni.cpp. Nel file funzioni.cpp vengono

definite altre funzioni utili al progetto quali

• void lettura(int (&dimensioni)[2],Eigen::Matrix<double, Eigen::Dynamic,

Eigen::Dynamic> &matrice_spostamenti_online). La funzione legge da file i

dati necessari per la POD

• void stampa_dati(const Eigen::Matrix<double, Eigen::Dynamic,

Eigen::Dynamic> &matrice). La funzione stampa su file i dati necessari per la

POD

• void dolfin2eigen(const dolfin::Vector &vector_dolfin, Eigen::Matrix<double,

Eigen::Dynamic, 1> &vettore_eigen,const dolfin::Matrix

&matrice_dolfin,Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>

&matrice_eigen). La funzione è utilizzata per trasferire i dati da una matrice di

Dolfin ad una matrice di Eigen

• void dolfin2eigen_sparsa(const dolfin::Vector &vector_dolfin,

Eigen::Matrix<double, Eigen::Dynamic, 1> &vettore_eigen,const

dolfin::Matrix &matrice_dolfin,Eigen::SparseMatrix<double>

&matrice_eigen). La funzione fa la stessa cosa della funzione precedente ma la

funzione eigen che riceve i dati è una matrice sparsa. È stata utilizzata questa

funzione per sfruttare la struttura sparsa delle matrice di Dolfin generate.

• 3 funzioni, utilizzate per i metodi solver della classe ElasticProblemSolver

precedentemente citata:

a) void solver_eigen(dolfin::Function &u, const Eigen::Matrix<double,

Eigen::Dynamic, Eigen::Dynamic> &matrice, const

Eigen::Matrix<double, Eigen::Dynamic, 1> &vettore_destro). Solver

che usa le matrici dense di Eigen

b) void solver_eigen_sparso(dolfin::Function &u, const

Eigen::SparseMatrix<double> &matrice, const Eigen::Matrix<double,

Eigen::Dynamic, 1> &vettore_destro, const Eigen::MatrixXd

&vettori_singolari_sinistri). Solver che usa le matrici sparse di Eigen

c) void solver_eigen_proiettato(dolfin::Function &u, const

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> &matrice,

const Eigen::Matrix<double, Eigen::Dynamic, 1> &vettore_destro,

const Eigen::MatrixXd &vettori_singolari_sinistri). Solver che usa la

POD

A questo file si aggiunge un file myfun.hh con le funzioni obiettivo utilizzate

nell’ottimizzazione.

3.2 Fase Offline e Fase Online

Il Progetto comprende una fase Offline ed una fase Online separate. La fase Offline è

stata utilizzata per generare le 𝑆 soluzioni ottime {𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑒(𝑖)

}𝑖=1𝑆 associate al

problema 24 con la volume fraction{𝛼𝑖}𝑖=1𝑆 . Le soluzioni ottime sono state calcolate

senza l’utilizzo della POD. Le {𝑢𝑛𝑢𝑚(𝑖)

}𝑖=1𝑆 calcolate vengono quindi stampate su un file

txt di nome “soluzione_spostamento.txt”. La fase Online invece prevede la lettura del

file “soluzione_spostamento.txt” e il calcolo di una nuova soluzione 𝑢𝑛𝑢𝑚𝑃𝑂𝐷 e 𝜌𝑃𝑂𝐷 per

un valore di valume fraction 𝛼𝑃𝑂𝐷 diverso dagli 𝑆 valori di volume fraction {𝛼𝑖}𝑖=1𝑆

utilizzati per risolvere il problema della fase Offline. 𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜𝑃𝑂𝐷 e 𝜌𝑜𝑡𝑡𝑖𝑚𝑜

𝑃𝑂𝐷 sono

calcolati sfruttando la proiezione nel sottospazio generato dalla SVD della matrice le

cui colonne sono le soluzioni {𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

}𝑖=1𝑆 ad ogni iterazione del ciclo di

iterazione.

3.2.1 Fase Offline

Avviata l’esecuzione del codice verrà chiesto se si vuole procedere con la fase Offline

(premendo 0) o con la fase Online (premendo un qualsiasi altro numero). Avviata la

fase Offline verrà chiesto di inserire la volume fraction 𝛼𝑚𝑖𝑛𝑖𝑚𝑎, la volume fraction

𝛼𝑚𝑎𝑠𝑠𝑖𝑚𝑎 e per quante volume fraction si vuole avere la soluzione (il programma

quindi calcolerà la soluzione per quel numero di volume fraction uniformemente

distribuite tra 𝛼𝑚𝑖𝑛𝑖𝑚𝑎 e 𝛼𝑚𝑎𝑠𝑠𝑖𝑚𝑎)La fase Offline è stata eseguita calcolando gli

spostamenti {𝑢𝑛𝑢𝑚(𝑖)

}𝑖=130 per 30 configurazioni ottime con 𝛼𝑚𝑖𝑛𝑖𝑚𝑎 = 0.25 e

𝛼𝑚𝑎𝑠𝑠𝑖𝑚𝑎 = 0.627. I valori di volume fraction usati sono illustrati in Figura 1

Figura 1: Volume fraction usate per la fase Offline

All’inizio della fase Offline viene definito:

• il dominio da ottimizzare (un Cantilever, Figura 2)

• Le costanti Fisiche 𝐸 e 𝜈, usati per calcolare i coefficienti di Lamè 𝜇 e 𝜆 come

𝜆 =𝐸𝜈

(1+𝜈)(1−𝜈) , 𝜇 =

𝐸

2(1+𝜈)

• Gli spazi funzionali

• Le condizioni al bordo

• I valori di volume fraction {𝛼𝑖}𝑖=130 per i quali viene risolto il problema 24

Figura 2: Cantilever ottimizzato

Una volta definito il problema viene creato un oggetto della classe

ElasticProblemSolver. Per risolvere il problema di ottimizzazione viene creato un

oggetto della classe nlopt::opt. Tramite il costruttore della classe nlopt::opt può essere

scelto l’algoritmo di ottimizzazione. L’algoritmo usato è il “local optimization,

derivative/gradient-based (Method of Moving Asymptotes)”. Tramite il metodo

set_min_objective della classe nlopt::opt si imposta la funzione obiettivo, che nel

nostro caso è la compliance, tramite la funzione myvfunc. myvfunc riceve in ingresso

in input una funzione di densità ed il vincolo di volume fraction per il quale si sta

calcolando la configurazione ottima, calcola il campo di spostamento associato alla

funzione di densità in ingresso tramite il solver ElasticProblemSolver::Solve() e la

direzione di discesa tramite il metodo ElasticProblemSolver::get_grad_compl(); fatto

questo calcola il valore di compliance tramite il metodo

ElasticProblemSolver::get_compl() e lo ritorna. Si impone il vincolo sulla volume

fraction tramite il metodo della classe nlopt::opt add_inequality_constraint. Il risultato

dell’ottimizzazione si ottiene infine attraverso il metodo optimize della classe

nlopt::opt. Da notare come nella fase Offline si calcoli il campo di spostamento ad ogni

iterazione di ottimizzazione tramite il metodo ElasticProblemSolver::Solve(), che non

prevede l’uso della POD. L’ottimizzazione quindi è svolta per ogni valore di volume

fraction {𝛼𝑖}𝑖=1𝑆 e calcola la distribuzione di densità ottima 𝜌𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑎

(𝑖) associata a 𝛼𝑖;

attraverso poi il solver della classe ElasticProblemSolver si ricavare il campo di

spostamento 𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

associato a 𝜌𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

. Il campo di spostamento 𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

viene quindi stampato sul file “soluzione_spostamento.txt”

3.2.2 Fase Online

Nella fase Online vengono usate le informazioni raccolte nel file

“soluzione_spostamento.txt” dalla fase Offline. Nella fase Online verrà chiesto il

valore di Volume fraction 𝛼 per il quale voler fare l’ottimizzazione, e poi si passerà

alla lettura della file “soluzione_spostamento.txt”, i cui dati verranno salvati nella

matrice matrice_spostamenti_online, un oggetto della classe:

Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>

La matrice matrice_spostamenti_online avrà nella colonna 𝑖 il campo di spostamento

𝑢𝑛𝑢𝑚,𝑜𝑡𝑡𝑖𝑚𝑜(𝑖)

associato alla densità ottima 𝜌𝑖 calcolato con il valore di volume fraction

𝛼𝑖.

Per fare la SVD viene creato un oggetto 𝑠𝑣𝑑 della classe:

Eigen::JacobiSVD<Eigen::MatrixXd>

nel cui costruttore viene passato la matrice matrice_spostamenti_online. Per estrarre

poi la matrice dei vettori singolari sinistri si usa il metodo matrixU() della classe

Eigen::JacobiSVD<Eigen::MatrixXd>. I vettori singolari sinistri sono salvati

nell’oggetto vettori_singolari_sinistri della classe Eigen::MatrixXd. Il prossimo

passaggio è scegliere la dimensione del sottospazio nel quale si proietterà il sistema in

equazione 27. Per decidere la dimensione verrà stampata una tabella a schermo del tipo

Nella prima colonna della tabella vengono mostrati i valori singolari in ordine

decrescente. Nella seconda colonna viene mostrata la 𝑑𝑒𝑠𝑐𝑟𝑖𝑧𝑖𝑜𝑛𝑒 𝑑𝑖 𝑣𝑎𝑟𝑖𝑎𝑛𝑧𝑎. La

descrizione di varianza in riga 𝑗 è uguale alla somma dei primi 𝑗 valori singolari fratto

la somma totale. 1 − 𝑑𝑒𝑠𝑐𝑟𝑖𝑧𝑖𝑜𝑛𝑒 𝑑𝑖 𝑣𝑎𝑟𝑖𝑎𝑛𝑧𝑎 è uguale all’errore di proiezione delle

soluzione storiche nel sottospazio di dimensione 𝑗. Il valore in riga 𝑗 della terza colonna

indica la grandezza della dimensione dello spazio di proiezione per avere la

𝑑𝑒𝑠𝑐𝑟𝑖𝑧𝑖𝑜𝑛𝑒 𝑑𝑖 𝑣𝑎𝑟𝑖𝑎𝑛𝑧𝑎 nela riga 𝑗. Grazie alla tabella si può scegliere la dimensione

più appropriata per lo spazio di proiezione.

I passaggi successivi definiscono il problema da risolvere, come nel caso Offline:

• il dominio da ottimizzare (un Cantilever, Figura x1)

• Le costanti Fisiche 𝐸 e 𝜈, usati per calcolare i coefficienti di Lamè 𝜇 e 𝜆 come

𝜆 =𝐸𝜈

(1+𝜈)(1−𝜈) , 𝜇 =

𝐸

2(1+𝜈)

• Gli spazi funzionali

• Le condizioni al bordo

• Il valore di volume fraction 𝛼𝑖 per il quale si vuole risolvere il problema 24.

E si passa quindi a definire il problema di ottimizzazione tramite la classe nlopt::opt.

La funzione obiettivo in questo caso però non è descritta dalla funzione myvfunc

utilizzata nella fase Offline, ma dalla funzione myvfunc_proiettata. myvfunc_proiettata

è analoga a myvfunc ma calcola il campo di spostamento associato ad una funzione di

densità (ad ogni iterazione del ciclo di ottimizzazione) tramite il solver

ElasticProblemSolver::Solve_proiettato(). La soluzione di densità calcolata dopo il

ciclo di ottimizzazione è salvata nel file 3DOpt_online.pvd

4 Risultati In questa sezione verranno mostrati i risultati dell’ottimizzazione topologica con la

POD e paragonati, in termini di tempo e soluzione, con l’algoritmo di ottimizzazione

classico. Nell’algoritmo classico ad ogni iterazione 𝑖 di ottimizzazione viene

calcolato il campo di spostamento 𝑢𝑛𝑢𝑚,𝑖 tramite la risoluzione di un sistema lineare

con una matrice con 126630009 valori, 11253 righe e 11253 colonne, con 450459

elementi non nulli. In Figura 3 è mostrato il risultato con l’algoritmo di

ottimizzazione classico con un valore di volume fraction pari a 0.57. Il tempo di

calcolo è 4 minuti e 38 secondi, la risoluzione del sistema lineare ad ogni ciclo è di

0.02 secondi, con una Compliance del risultato finale pari a 0.153855

Figura 3

In figura 4 invece è mostrata la soluzione con l’algoritmo di ottimizzazione usando la

POD. La dimensione del sottospazio usata è di 30 vettori, la dimensione massima

disponibile con 30 soluzioni storiche, ed ad ogni iterazione il tempo impiegato per la

soluzione del sistema lineare per il calcolo dello spostamento 𝑢 è di 0.0005 secondi,

40 volte di meno del precedente. Il tempo di calcolo di calcolo totale è di 2 minuti e 7

secondi, e la compliance finale è 0.120192

Figura 4

In Figura 5 viene mostrato invece il caso di una ottimizzazione topologica con POD

usando però con una base di dimensione 15. Con una base di questa dimensione

l’errore, riferito alla proiezione delle soluzioni storiche in questo sottospazio, è meno

dello 0.5%. Il tempo di calcolo dello spostamento 𝑢 ad ogni iterazione è di 0.0002

secondi, circa la metà del caso con una dimensione dello spazio di proiezione di 30. Il

tempo di calcolo totale è di 1 minuto e 39 secondi, meno del precedente.

Figura 5

Pur avendo scelto una dimensione dello spazio di proiezione metà del primo caso i

risultati sono molto simili. Questo è dovuto al fatto da una dimensione pari a 15 ad una

pari a 30 si passa da un errore dello 0.5% ad uno dello 0%.

Con una dimensione del sottospazio pari a 5 l’errore di proiezione delle soluzioni

storiche in tale sottospazio è del 2%. Il tempo per risolvere il sistema lineare è un terzo

del caso in cui il sottospazio abbia dimensione 15 e il tempo di calcolo globale è di 1

minuto e 35 secondi. Il risultato è mostrato in Figura 6.

Figura 6

In Figura 7 invece è mostrato il caso con una dimensione del sottospazio pari ad 1.

L’errore di proiezione delle soluzioni storiche in tale sottospazio è del 5%, il tempo di

risoluzione del sistema lineare (equazione) è di 3 10−5 secondi e il tempo totale

dell’algoritmo è di 1 minuto e 7 secondi.

Figura 6

5 Conclusioni Nel progetto è stata introdotto la POD e il problema di ottimizzazione topologica a minima

compliance, è stato introdotto l’algoritmo implementato per poter usare la POD nel ciclo di

ottimizzazione per il calcolo dello spostamento 𝑢 e sono stati presentati i risultati. I risultati

hanno mostrato un calo del tempo di calcolo dell’ottimizzazione con POD rispetto a quello

classico, è una buona scalabilità al variare della dimensione del sottospazio. In particolare al

diminuire della dimensione si aveva un guadagno in tempo con una soluzione praticamente

identica al caso di dimensioni del sottospazio più alte. Questo è dovuto alla distribuzione del

valori dei valori singolari sinistri.

6 Creazione ambiente, Compilazione ed avvio Il progetto è stato testato ed i risultati sono stati prodotti su una macchina virtuale

VirtualBox. L’installazione del sistema operativo è:

• Installazione di vagrant e git

• Con git nella cartella lanciare “vagrant init centos/7”

• Lanciare “vagrant up”

• Apri il file C:\Users\{username}\.vagrant.d\boxes\centos-VAGRANTSLASH-

7\1611.01\virtualbox\Vagrantfile e lo modifichi in modo che compaia

config.vm.synced_folder ".", "/vagrant", type: "virtualbox"

• Lanciare “vagrant plugin install vagrant-vbguest” nella cartella del punto 2

• Lanciare “vagrant up”

• Lanciare “vagrant ssh”

• Installazione dei moduli con “curl -O "http://www1.mate.polimi.it/PACS/gcc5-

glibc.tgz" ”

• sudo tar xvzf gcc5-glibc.tgz -C /

cat - >> ~/.bashrc <<EOF

export mkPrefix=/u/sw

source \$mkPrefix/etc/profile

module load gcc-glibc/5 octave

EOF

La macchina virtuale è pronta. Per avviarla bisogna lanciare:

• vagrant up

• vagrant ssh

• source /u/sw/etc/profile

• module load gcc-glibc/5 octave

• module load fenics

Al primo avvio bisogna scaricare ed installare NLopt-2.4.2.tar.gz dal sito:

http://nlopt.readthedocs.io/en/latest/

Di default, libnlopt.so sarà installato in usr/local/lib e nlopt.h in /usr/local/include. Se

le cartelle sono diverse bisogna correggere il cmake inserito del progetto, nelle righe

• set(NLOPTINC /usr/local/include)

• set(NLOPTLIB /usr/local/lib)

Per la compilazione ed esecuzione lanciare

• cmake .

• make

• ./3Delasticty

L’output del programma è 3DOpt_online.pvd, da aprire con paraview.

IMPORTANTE: Nella cartella ci sono due file: soluzioni_spostamento.txt

e dimensioni.txt . Questi file contengono le informazioni delle 30 simulazioni offline.

Se si lancia una nuova fase offline vengono sovrascritti. Per evitare che succeda basta

rinominarli prima che si lanci una nuova fase offline, in questo modo nel caso si voglia

lanciare una fase online con le stesse soluzioni storiche di questo paper è sufficiente

rinominarli con il nome giusto.

7 References S. Volkwein: Proper Orthogonal Decomposition: Theory and Reduced-Order

Modelling

M.P. Bendsoe, O. Sigmund: Topology optimization Theory, Methods and

Applications