Sistemi lineari particolari..
Ciao a tutti,
da ex studente di ingegneria vi pongo un problema che mi assale in questo ultimo periodo legato al tema in oggetto.
Dato il classico sistema lineare: Ax=b
dove A è una matrice nxn fatta in un modo particolare:
presenta solo numeri Reali positivi sulla diagonale e solo -1 in tutti le altre posizioni. Matematicamente parlando l'elemento aij dove i diverso da j è uguale a -1 mentre l'elemento aij dove i è uguale a j appartiene a R+.
Un esempio di A (3x3):
[3 -1 -1
-1 1,2 -1
-1 -1 75]
il problema per me, data la matrice A, è trovare i vettori b e x (di dimensioni n) tali per cui:
Tutti i singoli elementi xi siano >0
Tutti i singoli elementi bi siano >=0
e che:
La Sommatoria dei singoli elementi xi per i che va da 1 a n sia MINIMA
e contemporaneamente
La Sommatoria i singoli elementi bi per i che va da 1 a n sia MASSIMA
Spero di essermi spiegato..
Ovviamente i requisiti richiesti non saranno soddisfabili per tutte le possibili matrici A. Ad esempio quando molti elementi sulla diagonale sono uguali diventa un pò difficile matematicamente far tornare i conti.
Vorrei comunque condizionare meglio il problema per chiarirmi un pò le idee.
Qualcuno può aiutarmi?
Grazie a tutti per la disponibilità!
Saluti
D.
da ex studente di ingegneria vi pongo un problema che mi assale in questo ultimo periodo legato al tema in oggetto.
Dato il classico sistema lineare: Ax=b
dove A è una matrice nxn fatta in un modo particolare:
presenta solo numeri Reali positivi sulla diagonale e solo -1 in tutti le altre posizioni. Matematicamente parlando l'elemento aij dove i diverso da j è uguale a -1 mentre l'elemento aij dove i è uguale a j appartiene a R+.
Un esempio di A (3x3):
[3 -1 -1
-1 1,2 -1
-1 -1 75]
il problema per me, data la matrice A, è trovare i vettori b e x (di dimensioni n) tali per cui:
Tutti i singoli elementi xi siano >0
Tutti i singoli elementi bi siano >=0
e che:
La Sommatoria dei singoli elementi xi per i che va da 1 a n sia MINIMA
e contemporaneamente
La Sommatoria i singoli elementi bi per i che va da 1 a n sia MASSIMA
Spero di essermi spiegato..
Ovviamente i requisiti richiesti non saranno soddisfabili per tutte le possibili matrici A. Ad esempio quando molti elementi sulla diagonale sono uguali diventa un pò difficile matematicamente far tornare i conti.
Vorrei comunque condizionare meglio il problema per chiarirmi un pò le idee.
Qualcuno può aiutarmi?
Grazie a tutti per la disponibilità!
Saluti
D.
Risposte
Ci provo.
Prima di tutto bisogna notare che il tuo è un problema di ottimizzazione multiobiettivo, perché vuoi minimizzare $\sum_i x_i$ e contemporaneamente massimizzare $\sum_i b_i$. Le possibilità sono due: o riesci a formulare il problema come un ottimo di Pareto, ma qui non ti so aiutare, oppure accetti di assegnare dei pesi ai due ottimi in modo da ridurti ad averne uno solo: $\min (W_x \sum_i x_i - W_b \sum_i b_i)$.
In questo caso io proverei a procedere come segue:
prima di tutto, visto che sia le $x$ che le $b$ sono incognite, evidenzia il fatto che si tratta di un sistema omogeneo sottodeterminato riscrivendolo come $[[3,-1,-1,-1,0,0],[-1,1.2,-1,0,-1,0],[-1,-1,75,0,0,-1]] [[x_1],[x_2],[x_3],[b_1],[b_2],[b_3]]=0$
Questa matrice per com'è fatta ha rango 3. Scomponila in valori singolari (SVD) come $U \cdot W \cdot V^T$, con $U$ matrice 3 x 6, $W$ matrice diagonale dei valori singolari e $V$ matrice 6 x 6.
Tre dei valori singolari saranno nulli (entro la precisione macchina). Le corrispondenti colonne della matrice $V$, e siano $\mathbf{V}_1,\mathbf{V}_2,\mathbf{V}_3$, costituiscono una base del nullspace della matrice, quindi tutte le soluzioni si possono esprimere come combinazione lineare $[[x_1],[x_2],[x_3],[b_1],[b_2],[b_3]]=\lambda_1 \mathbf{V}_1 + \lambda_2 \mathbf{V}_2 + \lambda_3 \mathbf{V}_3$. I coefficienti della combinazione lineare, $\lambda_1, \lambda_2, \lambda_3$, sono le nuove variabili.
Il problema ora è diventato
$ min_{ \lambda_1, \lambda_2, \lambda_3} \lambda_1 (W_x(V_{11}+V_{21}+V_{31})-W_b(V_{41}+V_{51}+V_{61}))+\lambda_2 (W_x(V_{12}+V_{22}+V_{32})-W_b(V_{42}+V_{52}+V_{62}))$
$+\lambda_3 (W_x(V_{13}+V_{23}+V_{33})-W_b(V_{43}+V_{53}+V_{63}))$
(dove i coefficienti di $\lambda_1, \lambda_2, \lambda_3$ sono costanti)
soggetto a
$\lambda_1 V_{11}+\lambda_2 V_{12}+\lambda_3V_{13}>0$
$\lambda_1 V_{21}+\lambda_2 V_{22}+\lambda_3V_{23}>0$
$\lambda_1 V_{31}+\lambda_2 V_{32}+\lambda_3V_{33}>0$
$\lambda_1 V_{41}+\lambda_2 V_{42}+\lambda_3V_{43}>=0$
$\lambda_1 V_{51}+\lambda_2 V_{52}+\lambda_3V_{53}>=0$
$\lambda_1 V_{61}+\lambda_2 V_{62}+\lambda_3V_{63}>=0$
Se rinunci alla diseguaglianza stretta per le prime tre condizioni, questo è un problema di programmazione semidefinita. Ci sono diversi pacchetti software che lo risolvono; ad esempio SDPA.
PS ...ma da dove salta fuori quella matrice (se è lecito chiedere)?
Prima di tutto bisogna notare che il tuo è un problema di ottimizzazione multiobiettivo, perché vuoi minimizzare $\sum_i x_i$ e contemporaneamente massimizzare $\sum_i b_i$. Le possibilità sono due: o riesci a formulare il problema come un ottimo di Pareto, ma qui non ti so aiutare, oppure accetti di assegnare dei pesi ai due ottimi in modo da ridurti ad averne uno solo: $\min (W_x \sum_i x_i - W_b \sum_i b_i)$.
In questo caso io proverei a procedere come segue:
prima di tutto, visto che sia le $x$ che le $b$ sono incognite, evidenzia il fatto che si tratta di un sistema omogeneo sottodeterminato riscrivendolo come $[[3,-1,-1,-1,0,0],[-1,1.2,-1,0,-1,0],[-1,-1,75,0,0,-1]] [[x_1],[x_2],[x_3],[b_1],[b_2],[b_3]]=0$
Questa matrice per com'è fatta ha rango 3. Scomponila in valori singolari (SVD) come $U \cdot W \cdot V^T$, con $U$ matrice 3 x 6, $W$ matrice diagonale dei valori singolari e $V$ matrice 6 x 6.
Tre dei valori singolari saranno nulli (entro la precisione macchina). Le corrispondenti colonne della matrice $V$, e siano $\mathbf{V}_1,\mathbf{V}_2,\mathbf{V}_3$, costituiscono una base del nullspace della matrice, quindi tutte le soluzioni si possono esprimere come combinazione lineare $[[x_1],[x_2],[x_3],[b_1],[b_2],[b_3]]=\lambda_1 \mathbf{V}_1 + \lambda_2 \mathbf{V}_2 + \lambda_3 \mathbf{V}_3$. I coefficienti della combinazione lineare, $\lambda_1, \lambda_2, \lambda_3$, sono le nuove variabili.
Il problema ora è diventato
$ min_{ \lambda_1, \lambda_2, \lambda_3} \lambda_1 (W_x(V_{11}+V_{21}+V_{31})-W_b(V_{41}+V_{51}+V_{61}))+\lambda_2 (W_x(V_{12}+V_{22}+V_{32})-W_b(V_{42}+V_{52}+V_{62}))$
$+\lambda_3 (W_x(V_{13}+V_{23}+V_{33})-W_b(V_{43}+V_{53}+V_{63}))$
(dove i coefficienti di $\lambda_1, \lambda_2, \lambda_3$ sono costanti)
soggetto a
$\lambda_1 V_{11}+\lambda_2 V_{12}+\lambda_3V_{13}>0$
$\lambda_1 V_{21}+\lambda_2 V_{22}+\lambda_3V_{23}>0$
$\lambda_1 V_{31}+\lambda_2 V_{32}+\lambda_3V_{33}>0$
$\lambda_1 V_{41}+\lambda_2 V_{42}+\lambda_3V_{43}>=0$
$\lambda_1 V_{51}+\lambda_2 V_{52}+\lambda_3V_{53}>=0$
$\lambda_1 V_{61}+\lambda_2 V_{62}+\lambda_3V_{63}>=0$
Se rinunci alla diseguaglianza stretta per le prime tre condizioni, questo è un problema di programmazione semidefinita. Ci sono diversi pacchetti software che lo risolvono; ad esempio SDPA.
PS ...ma da dove salta fuori quella matrice (se è lecito chiedere)?
Facendo alcune prove ho notato che affinché il vettore x risulti con tutte le componenti maggiori o uguali a 0 è necessario che il Determinante della matrice A sia maggiore di zero. Infatti, indipendentemente dalle dimensioni della matrice A (3x3, 10x10, 26x26, ecc) se det A è negativo, indipendentemente dal vettore b: che potrebbe essere fissato per comodità a tutti 1, ritrovo tutti gli elementi xi negativi.
Ciò significa che, tra tutte le possibili matrici A, ai fini della soluzione del problema è necessario individuare quelle con Det A maggiore di 0. Considerando la tipicità di tale matrice, cioè tutto -1 salvo gli elementi sulla diagonali che sono numeri reali positivi, come faccio a capire più o meno velocemente quando tale matrice si presta a realizzare un vettore x positivo ?
Esiste un metodo che mi potrebbe aiutare?
Grazie mille ancora a chi vorrà rispondermi.
d.
Ciò significa che, tra tutte le possibili matrici A, ai fini della soluzione del problema è necessario individuare quelle con Det A maggiore di 0. Considerando la tipicità di tale matrice, cioè tutto -1 salvo gli elementi sulla diagonali che sono numeri reali positivi, come faccio a capire più o meno velocemente quando tale matrice si presta a realizzare un vettore x positivo ?
Esiste un metodo che mi potrebbe aiutare?
Grazie mille ancora a chi vorrà rispondermi.
d.
Mi sto perdendo. Quali sono le variabili su cui si può agire? Perché se fissi i $b$ tutti a 1, e visto che mi sembra di capire che la matrice $A$ è data (e non è singolare), allora le $x$ vengono quel che vengono, non c'è margine di manovra.
Ci sono varie matrici A.
Per alcune, quelle con Det A positivo, provando a fissare bi ad esempio a 1 trovo il vettore xi, come dici tu, quel che l'é. In questo caso vale (credo) il metodo di ottimizzazione multiobiettivo per cui ci dovrebbero essere margini di manovra sui due vettori x e b.
Per altre matrici A, per quanto si possa provara a pacioccare sui termini xi e bi, se Det A è negativo non se ne esce.. Se poi, anche in questo caso, riesci a trovarmi i due vettori x e b che soddisfano le condizione richieste (xi>0 e bi>=0) e risolvono il sistema Ax=b allora sei un drago!!
Ti torna?
Per alcune, quelle con Det A positivo, provando a fissare bi ad esempio a 1 trovo il vettore xi, come dici tu, quel che l'é. In questo caso vale (credo) il metodo di ottimizzazione multiobiettivo per cui ci dovrebbero essere margini di manovra sui due vettori x e b.
Per altre matrici A, per quanto si possa provara a pacioccare sui termini xi e bi, se Det A è negativo non se ne esce.. Se poi, anche in questo caso, riesci a trovarmi i due vettori x e b che soddisfano le condizione richieste (xi>0 e bi>=0) e risolvono il sistema Ax=b allora sei un drago!!
Ti torna?
Più che il determinante c'entrano gli autovalori; la tua matrice di esempio ha autovalori 75.0271, 3.44181, 0.731135, tutti e tre positivi, e quindi è definita positiva; quindi $\mathbf{x}^T \cdot A \mathbf{x} >0$ per ogni $\mathbf{x}$ non nullo.
Se come $\mathbf{x}$ prendiamo la soluzione del sistema $A \mathbf{x}=\mathbf{b}$ allora risulta $\mathbf{x}^T A \mathbf{x}= \mathbf{x}^T \mathbf{b} >0$, e poiché $\mathbf{b}=[[1],[1],[1]]$, è come dire $\sum_i x_i >0$.
Questo non garantisce ancora che le singole componenti $x_i$ siano $>0$, però.
Poi il determinante è il prodotto degli autovalori, e quindi è positivo se tutti gli autovalori sono positivi. Ma è positivo anche se un numero pari di autovalori sono negativi.
Quindi secondo me queste tue osservazioni empiriche non ti portano da nessuna parte.
Valuta se ti può andar bene 'pesare' i due ottimi e, se sì, posta una matrice significativa; con un momento di calma (che potrebbe voler dire anche qualche giorno) provo a farti i conti.
Se poi tu potessi dirci da dove viene il problema, non è escluso che lo si possa riformulare in un modo più simpatico.
Se come $\mathbf{x}$ prendiamo la soluzione del sistema $A \mathbf{x}=\mathbf{b}$ allora risulta $\mathbf{x}^T A \mathbf{x}= \mathbf{x}^T \mathbf{b} >0$, e poiché $\mathbf{b}=[[1],[1],[1]]$, è come dire $\sum_i x_i >0$.
Questo non garantisce ancora che le singole componenti $x_i$ siano $>0$, però.
Poi il determinante è il prodotto degli autovalori, e quindi è positivo se tutti gli autovalori sono positivi. Ma è positivo anche se un numero pari di autovalori sono negativi.
Quindi secondo me queste tue osservazioni empiriche non ti portano da nessuna parte.
Valuta se ti può andar bene 'pesare' i due ottimi e, se sì, posta una matrice significativa; con un momento di calma (che potrebbe voler dire anche qualche giorno) provo a farti i conti.
Se poi tu potessi dirci da dove viene il problema, non è escluso che lo si possa riformulare in un modo più simpatico.
Ecco qui un esempio di matrice A significativa. E' di dimensioni 27x27.
Di seguito gli elementi sulla diagonale principale. Tutti gli altri li ometto perchè come sai sono tutti -1.
5,5
7
7,5
12
11
24
24
22
36
74
14
26
17
179
28
47
239
269
119
139
9
7,5
15
69
239
10
Attendo tue news..
Grazie
d.
Di seguito gli elementi sulla diagonale principale. Tutti gli altri li ometto perchè come sai sono tutti -1.
5,5
7
7,5
12
11
24
24
22
36
74
14
26
17
179
28
47
239
269
119
139
9
7,5
15
69
239
10
Attendo tue news..
Grazie
d.
Hai dato solo 26 elementi.
va bene anche solo con quelli
Ho giocato un po' con il tuo problema e ti anticipo già che la conclusione non ti piacerà; ma andiamo con ordine.
Non è necessario ricorrere a SDP, come avevo detto; il problema si riconduce alla vecchia buona programmazione lineare (LP).
Un pacchetto open source per la programmazione lineare è GLPK; io me lo sono installato su Linux, ma credo che esista anche per Windows.
GLPK gestisce anche vincoli di eguaglianza stretta, e questo fa cadere la necessità di lavorare nel sottospazio vincolare per rendere implicito il vincolo: $A\mathbf{x}-\mathbf{b}=\mathbf{0}$ può essere gestito direttamente da GLPK assieme alle condizioni sulle 'variabili strutturali' $x_{i} >= 0, b_{i} >= 0$. Quindi non serve la decomposizione in valori singolari.
A parte questo, ho seguito la falsariga delineata mel mio primo post. Considero variabili sia $\mathbf{x}$ che $\mathbf{b}$ e quindi amplio il sistema in $[[5.5,-1,...,-1,-1,0,...,0],[...,...,...,...,...,...,...,...],[-1,...,-1,10,0,...,0,-1]] [[x_1],[...],[b_{26}]]=[[0],[...],[0]]$. Introduco una unica funzione obiettivo pesata sulle due funzioni obiettivo che ti interessano.
Qui c'è il codice c++ che ho scritto:
Puoi compilarlo con
dopo aver installato la libreria glpk (io uso ubuntu 12.04 e l'ho trovata nei repository). E' ampiamente commentato e dovrebbe risultarti facile cambiare la matrice.
Qual è il problema, che non ti piacerà? Che con i vincoli $x_i >= 0, b_i >= 0$ l'ottimo è banale: $\mathbf{x}=\mathbf{0}, \mathbf{b}=\mathbf{0}$.
Non è necessario ricorrere a SDP, come avevo detto; il problema si riconduce alla vecchia buona programmazione lineare (LP).
Un pacchetto open source per la programmazione lineare è GLPK; io me lo sono installato su Linux, ma credo che esista anche per Windows.
GLPK gestisce anche vincoli di eguaglianza stretta, e questo fa cadere la necessità di lavorare nel sottospazio vincolare per rendere implicito il vincolo: $A\mathbf{x}-\mathbf{b}=\mathbf{0}$ può essere gestito direttamente da GLPK assieme alle condizioni sulle 'variabili strutturali' $x_{i} >= 0, b_{i} >= 0$. Quindi non serve la decomposizione in valori singolari.
A parte questo, ho seguito la falsariga delineata mel mio primo post. Considero variabili sia $\mathbf{x}$ che $\mathbf{b}$ e quindi amplio il sistema in $[[5.5,-1,...,-1,-1,0,...,0],[...,...,...,...,...,...,...,...],[-1,...,-1,10,0,...,0,-1]] [[x_1],[...],[b_{26}]]=[[0],[...],[0]]$. Introduco una unica funzione obiettivo pesata sulle due funzioni obiettivo che ti interessano.
Qui c'è il codice c++ che ho scritto:
#include <iostream>
#include <glpk.h>
using namespace std;
int main ()
{
// Parametri del problema
// Numero di equazioni del sistema
const size_t m = 26;
// Elementi sulla diagonale
double diagonal [m] = { 5.5, 7., 7.5, 12., 11., 24., 24., 22., 36., 74., 14.,
26., 17., 179., 28., 47., 239., 269., 119., 139., 9., 7.5, 15., 69., 239., 10. };
// Numero di variabili del sistema (sia x che b)
const size_t n = 2 * m;
// Peso del sottoobiettivo "somma delle x" nell'obiettivo globale.
const double weight_x = 1.;
// Peso del sottoobiettivo "somma delle b" nell'obiettivo globale.
const double weight_b = 1.;
// Alloco un oggetto problema di glp.
glp_prob * lp = glp_create_prob ();
glp_set_obj_dir (lp, GLP_MIN); // minimizzazione
// Vincoli delle 'variabili strutturali' e loro coefficienti nella funzione obiettivo
glp_add_cols (lp, n);
for (size_t j = 0; j < n; j++) {
glp_set_col_bnds (lp, 1 + j, GLP_LO, 0., 0.);
if (j < m) glp_set_obj_coef (lp, 1 + j, weight_x);
else glp_set_obj_coef (lp, 1 + j, - weight_b);
}
// Ogni equazione è un vincolo di eguaglianza.
glp_add_rows (lp, m);
for (size_t i = 0; i < m; i++) {
glp_set_row_bnds (lp, 1 + i, GLP_FX, 0., 0.);
}
// Matrice dei vincoli
size_t constraint_matrix_size = m * m + m;
int constraint_matrix_row [1 + constraint_matrix_size];
int constraint_matrix_col [1 + constraint_matrix_size];
double constraint_matrix_element [1 + constraint_matrix_size];
size_t element_index = 0;
for (size_t i = 0; i < m; i++) {
for (size_t j = 0; j < m; j++) {
element_index ++;
constraint_matrix_row[element_index] = 1 + i;
constraint_matrix_col[element_index] = 1 + j;
if (i == j) constraint_matrix_element[element_index] = diagonal[i];
else constraint_matrix_element[element_index] = -1.;
}
element_index++;
constraint_matrix_row[element_index] = 1 + i;
constraint_matrix_col[element_index] = 1 + m + i;
constraint_matrix_element[element_index] = -1.;
}
// Carico la matrice dei vincoli nel problema
glp_load_matrix (lp, element_index, constraint_matrix_row, constraint_matrix_col,
constraint_matrix_element);
// Risolvo
int outcome = glp_simplex (lp, 0);
if (outcome != 0) {
cout << "glp_simplex ha ritornato " << outcome << " (non buono)." << endl;
}
cout << "Valore ottimo dell'obiettivo = " << glp_get_obj_val (lp) << endl;
cout << endl;
for (int i = 0; i < m; i++) {
cout << "Valore ottimo di x_" << 1 + i << " = "
<< glp_get_col_prim (lp, 1 + i) << endl;
}
for (int i = 0; i < m; i++) {
cout << "Valore ottimo di b_" << 1 + i << " = "
<< glp_get_col_prim (lp, 1 + m + i) << endl;
}
return 0;
}
Puoi compilarlo con
g++ -o nome_eseguibile nome_sorgente.cpp -lglpk -lm
dopo aver installato la libreria glpk (io uso ubuntu 12.04 e l'ho trovata nei repository). E' ampiamente commentato e dovrebbe risultarti facile cambiare la matrice.
Qual è il problema, che non ti piacerà? Che con i vincoli $x_i >= 0, b_i >= 0$ l'ottimo è banale: $\mathbf{x}=\mathbf{0}, \mathbf{b}=\mathbf{0}$.
si, dopo qualche analisi ho verifica che il sistema potrebbe avere una soluzione solo se la somma degli inversi degli elementi sulla diagonale è minore di 1..
Accedi a tutti gli appunti
Tutor AI: studia meglio e in meno tempo