[Matlab] Résoudre un système linéaire itératif

Bonjour à tous,

La semaine dernière je vous demandais comment résoudre un système non-linéaire relativement simple, cette semaine c’est un système linéaire itératif (car il est implicite). Pour faire simple, je dois résoudre un système assez gros (j’ai pas compté mais ça doit faire une quarantaines d’équations) mais de façon itérative… C’est-à-dire que je fais une estimation d’un paramètre (ici, la température) et je recalcule tous les paramètres jusqu’à avoir convergence (mon critère est: abs(T_old - T_new) < 0.1).

J’ai pas encore tenté d’écrire le code mais j’ai déjà toutes mes équations (je pense) et donc je me demandais quelle était l’approche la plus simple et plus efficace pour le faire ? Comme vous le savez, j’y connais très peu en programmation / Matlab donc ça m’aiderait d’avoir une explication de comment le coder!

Je mettrai mon code après parce qu’il risque de toute façon d’avoir des pépins

Merci d’avance!

Edit:

D’un point de vue de la structure du code ça devrait je pense ressembler à ça (conceptuellement)

BLOC AVEC LE SYSTEME D’EQ LINEAIRE (44) Ax = B

À LA FIN DU SYSTÈME, CALCULER DES VARIABLES AVEC LES SOLUTIONS DE CE SYSTEME

CALCULER DE NOUVELLES VARIABLES (POUR OBTENIR MA TEMPERATURE)

CONDITION "IF" T_ANCIEN - T_NOUVEAU < 0.1 => OK, LES SOLUTIONS SONT LES BONNES

SINON ("ELSE") -> RECALCULER

(J’ai mis en majuscule pour montrer que je vois ça sous forme de différents blocs / étapes).

Compliqué à mettre en oeuvre?

Salut,

je dois résoudre un système assez gros (j’ai pas compté mais ça doit faire une quarantaines d’équations) mais de façon itérative

C’est ridiculement petit, résoudre ça de façon directe serait beaucoup plus malin. Il t’est imposé d’utiliser une méthode itérative ? Est-ce que ta matrice est particulière (symétrique et/ou creuse par exemple) ?

Sinon, les méthodes itératives sont déjà implémentées, voir lsqr, gmres, minres pour ne citer que les plus générales (hésite pas à piocher dans "See Also" à la fin des docs de ces fonctions). Remarque qu’elles signalent toutes qu’elles sont pensées pour des matrices grandes (40, c’est pas grand) et creuses (avec plein de zéros).

Ça m’est imposé car c’est un problème d’ingénierie et je dois estimer une température dans une colonne de distillation sauf que je sais pas si la bonne mais comme je connais le flux de produit que je veux sortir, je dois faire en sorte d’avoir les bonnes températures.

J’ai pas complètement écrit la matrice la matrice mais elle est assez creuse (pas symétrique par contre). Je pourrais la copier ici mais je sais pas si c’est une bonne idée au vu de la taille de celle-ci

Je vais regarder les fonctions que tu proposes!

Ça m’est imposé car c’est un problème d’ingénierie et je dois estimer une température dans une colonne de distillation sauf que je sais pas si la bonne mais comme je connais le flux de produit que je veux sortir, je dois faire en sorte d’avoir les bonnes températures.

Je suis en train de me dire qu’on ne parle peut être pas la même langue, alors je vais essayer de mieux comprendre ton problème.

Est-ce que tu dois résoudre un système du genre $MT=R$ avec $M$ une matrice constante, $T$ ton vecteur d’inconnues et $R$ un vecteur quelconque constant ? Si oui, je ne vois pas en quoi le problème concret que tu tentes de résoudre va imposer le choix d’une méthode itérative.

Ou est-ce que $M$ et $R$ dépendent de $T$ auquel cas ton problème n’est plus linéaire et le choix d’une méthode itérative serait plus compréhensible.

Plutôt que copier la matrice, t’aurais pas une forme générale pour les équations vérifiées par la température ?

Je vais expliquer le problème d’abord pour être plus clair sinon ça reste très flou pour vous. Je dois résoudre un problème de chimie industrielle pour la production d’un polymère (styrène). J’ai deux réactions qui ont lieu dans le réacteur et mes trois derniers éléments (Flash, Dist 1 et Dist 2) sont des unités de séparations où je ne connais pas la température. Cependant je connais la production de styrène que je veux obtenir à la fin. Vu que j’ai cette infos (et d’autres comme la pression etc.), je peux écrire pleins d’équations de conservation de masse etc.

Voici quelques équations (je les mets pas toutes mais elles se ressemble concernant la forme):

$\mu _0^{EB} = 5\mu _0^B$

$\sum\limits_i {\mu _0^iM(i) = } \sum\limits_i {\mu _{3.1}^iM(i)} + \sum\limits_i {\mu _{4.2}^iM(i) + \sum\limits_i {\mu _{5.1}^iM(i)} }$

(pour 5 "i" différents car j’ai 5 espèces au total, i = EB, CY, H, ST, H ; M(i) c’est des constantes [masse molaire]).

$\mu _2^{EB} = \mu _1^{EB}(1 - {X_A})$

$\mu _2^H = \mu _1^H + \mu _1^{EB}{X_A} - 3\mu _1^H{X_B}$

(X_A et X_B sont des constantes connues)

Puis, dans la première unité de séparation [3. FLASH sur le schéma] on enlève "H". Donc, faut imaginer qu’après 3. on a plus que quatre espèces, ie.

$\mu _{4.1}^H = \mu _{4.2}^H = 0$ (pareil pour le flux 5)

Finalement, on a les équations qui vont faire intervenir "indirectement" la température au niveau des unités de séparation (3, 4 et 5). Voici un exemple et c’est pareil pour toutes les espèces:

$\mu _{4.1}^i = \mu _{3.2}^i\xi _{4}^i(T_{4}^D,T_{4}^B)$

où $\xi (...)$ est une valeur qui est comprise entre 0 et 1 (c’est simplement le pourcentage de ce qui va se retrouver en haut de l’unité de séparation pour l’espèce $i$). Ce $\xi (...)$ dépend de la volatilé relative et donc de la température. Ces expressions ne sont plus linéaires par contre.

On définit une volatilité relative

${\overline \alpha _{ij}} = \sqrt {\alpha _{ij}^{bas}\alpha _{ij}^{haut}} = \sqrt {\frac{{{P_{s,i}}({T^{bas}}){P_{s,i}}({T^{haut}})}}{{{P_{s,j}}({T^{bas}}){P_{s,j}}({T^{haut}})}}}$

Ce que je peux encore écrire en fonction des températures sachant que

$\log {P_s} = A - \frac{B}{{T + C}}$

où $A$, $B$ et $C$ sont des constantes que je connais.

Là, on imagine qu’on a fait une itération… Puis je peux calculer mes compositions (fractions molaires car je connais tous les flux, c’est bêtement relié) et calculer une nouvelle température ${T_{new}} = \sum\limits_i {{x_i}{P_{s,i}}}$ .

L’exercice me donne un critère de convergence $\left| {{T_{5.1,old}} - {T_{5.1,new}}} \right| < 0.1[K]$

Et mes inconnues sont tous les flux $\mu$ (sauf pour ST car je connais la quantité que je veux produire).

J’espère que c’est un peu plus clair. Je sais que c’est assez long et pourtant j’ai réduit au maximum le nombre d’équations… Ça m’a pris pas mal de temps à fixer toutes ces équations, maintenant reste à l’implémenter dans Matlab ce qui a l’air plus complexe que je ne le pensais

Je vais être honnête, je comprends pas grand chose à ces équations ni à tes notations. Mais si j’essaye de résumer, tu as des équations liant les $\mu$ (et non les $T$ comme je croyais au début…) avec des coefficients qui sont des fonctions des $T$, elles même fonctions des $\mu$. Donc ton système n’est pas linéaire, en fait, mais de la forme $M(\mu)\mu = R$.

Si je comprends bien, c’est dans un cadre scolaire puisque tu parles d’exercice. Vu la façon dont tu présentes les choses, j’ai l’impression que ton système se comporte de façon sympathique et la solution $\mu$ est un point fixe de la suite définit par récurrence par $M(\mu_n)\mu_{n+1} = R_n$. Donc probablement que ce qui est attendu est juste de construire la suite des $\mu_n$ à coup de linsolve jusqu’à ce que $T(\mu_{n+1}) - T(\mu_n)$ converge à la précision souhaitée.

Exactement! Désolé si c’était très ambigu. Je vais tenter de faire un linslove avec la documentation MATLAB et revenir ici si j’ai des soucis. Mais du coup, si c’est un linsolve c’est bien un système linéaire, non?

Est-ce qu’un fsolve ferait l’affaire aussi ? C’est simplement pour éviter de mettre les équations sous forme de matrice car ça risque d’être compliqué vu le nombre d’équations et de paramètres.

Merci @adri1.

Je voyais pas trop comment créer une matrice avec toutes ces variables mais j’ai écrit le code avec $fsolve$. Il manque encore des données (les pressions de saturation que je dois chercher avec les coefficients d’Antoine) et le calcul à la fin de mes fractions molaires. Cependant, je vois pas du tout comment faire une boucle / résoudre de manière itérative ce système.

Je laisse quand même toutes les équations pour les curieux même si elles sont pas vraiment utile pour répondre à la question (je suppose). J’ai 45 équations et 45 variables à trouver. Mes paramètres qui dépendent de $T$ sont les $\xi$ (enfin, indirectement car ça dépend de la volatilité relative qui elle dépend de la pression saturante donc de $T$). J’ai plusieurs températures (4 car deux colonnes et 2 températures par colonnes) mais l’exercice me demande de mettre comme critère de convergence uniquement $\left| {{T_{5.1,old}} - {T_{5.1,new}}} \right| < 0.1$ ce qui me perturbe aussi…

J’espère que c’est pas trop bordélique mais tout commentaire est bon à prendre pour que je m’améliore :-)

function fval = root11d(x)
%B = Benzene, ST = Styrene, EB = Ethylbenzene, H = Hydrogen, CY =
%Cyclohexane

%Molar masses in kg/mol
MW_B = (78.11)/1000; 
MW_ST = (104.15)/1000;
MW_EB = (106.17)/1000;
MW_H = (2.01)/1000;
MW_CY = (84.16)/1000;

%Conversion rates of rxn 

X_A = 0.5;
X_B = 0.95;

%Splitting ratio xi 

%Unit 4
xi_4_lk = 0.99; %EB
xi_4_hk = 0.01; %ST 
P_s_4_lk_B = 10^(-;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_4_lk_V = ;
P_s_4_hk_B = ;
P_s_4_hk_V = ;
P_s_4_B_B = ; 
P_s_4_B_V = ; 
P_s_4_CY_B = ;
P_s_4_CY_V =;
P_4 = 0.05; %pressure of the col 4 [bar]

alpha_4_lkhk = ((P_s_4_lk_B*P_s_4_lk_V)/(P_s_4_hk_B*P_s_4_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_4_min = log(xi_4_lk*(1-xi_4_hk)/xi_4_hk*(1-xi_4_lk))/log(alpha_4_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((P_s_4_B_B*P_s_4_B_V)/(P_s_4_hk_B*P_s_4_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_4_CYhk = ((P_s_4_CY_B*P_s_4_CY_V)/(P_s_4_hk_B*P_s_4_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_4_B = ((alpha_4_Bhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_Bhk-1)*xi_4_hk);
xi_4_ST = 0.01
xi_4_EB = 0.99
xi_4_CY = ((alpha_4_CYhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_CYhk-1)*xi_4_hk);

%Unit 5

xi_5_lk = 0.99; %CY
xi_5_hk = 0.01; %EB
P_s_5_lk_B = ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_5_lk_V = ;
P_s_5_hk_B = ;
P_s_5_hk_V = ;
P_s_5_B_B = ; 
P_s_5_B_V = ; 
P_s_5_ST_B = ;
P_s_5_ST_V =;
P_5 = 1; %pressure of the col 5 [bar]

%Variables Def / Notation: u_flow_spcies, eg. u_0_B is the flow rate "0" of B (benzene)
u_0_B = x(1); 
u_0_ST = x(2); 
u_0_EB = x(3); 
u_0_H = x(4); 
u_0_CY = x(5);

u_1_B = x(6); 
u_1_ST = x(7); 
u_1_EB = x(8); 
u_1_H = x(9); 
u_1_CY = x(10);

u_2_B = x(11); 
u_2_ST = x(12); 
u_2_EB = x(13); 
u_2_H = x(14); 
u_2_CY = x(15);

u_31_B = x(16); 
u_31_ST = x(17); 
u_31_EB = x(18); 
u_31_H = x(19); 
u_31_CY = x(20);

u_32_B = x(21); 
u_32_ST = x(22); 
u_32_EB = x(23); 
u_32_H = x(24); 
u_32_CY = x(25);

u_41_B = x(26); 
u_41_ST = x(27); 
u_41_EB = x(28); 
u_41_H = x(29); 
u_41_CY = x(30);

u_42_B = x(31); 
u_42_ST = x(32); 
u_42_EB = x(33); 
u_42_H = x(34); 
u_42_CY = x(35);

u_51_B = x(36); 
u_51_ST = x(37); 
u_51_EB = x(38); 
u_51_H = x(39); 
u_51_CY = x(40);

u_52_B = x(41); 
u_52_ST = x(42); 
u_52_EB = x(43); 
u_52_H = x(44); 
u_52_CY = x(45);


%Define sys of eq 

fval(1,1) = 5*u_0_B - mu_0_EB;
fval(2,1) = u_42_B*MW_B + u_32_B*MW_B + u_51_B*MW_B - u_0_B*MW_B ...
    + u_42_ST*MW_ST + u_32_ST*MW_ST + u_51_ST*MW_ST - u_0_ST*MW_ST ...
    + u_42_EB*MW_EB + u_32_EB*MW_EB + u_51_EB*MW_EB - u_0_EB*MW_EB ...
    + u_42_H*MW_H + u_32_H*MW_H + u_51_H*MW_H - u_0_H*MW_H ...
    + u_42_CY*MW_CY + u_32_CY*MW_CY + u_51_CY*MW_CY - u_0_CY*MW_CY;
fval(3, 1) = u_0_B + u_52_B - u_1_B; 
fval(4, 1) = u_0_ST + u_52_ST - u_1_ST;
fval(5, 1) = u_0_EB + u_52_EB - u_1_EB;
fval(6, 1) = u_0_H + u_52_H - u_1_H;
fval(7, 1) = u_0_CY + u_52_CY - u_1_CY;
fval(8, 1) = u_1_EB*(1-X_A)-u_2_EB;
fval(9, 1) = u_1_ST + u_1_EB*X_A -u_2_ST;
fval(10, 1) = u_1_B*X_B-u_2_CY+u_1_H;
fval(11, 1) = u_1_EB*X_A -3*u_1_B*X_B -u_2_H + u_1_H; 
%Flash Eq 
fval(12, 1) = u_31_B; 
fval(13, 1) = u_31_ST;
fval(14, 1) = u_31_EB; 
fval(15, 1) = u_31_CY;
fval(16, 1) = u_2_H - u_31_H;
fval(17, 1) = u_2_B - u_3_B; 
fval(18, 1) = u_2_ST - u_3_ST;
fval(19, 1) = u_2_EB - u_3_EB; 
fval(20, 1) = u_2_CY - u_3_CY;
%Distilation Column 1
fval(21, 1) = u_32_B*xi_4_B - u_41_B; %xi to be defined!!
fval(22, 1) = u_32_ST*xi_4_ST - u_41_ST;
fval(23, 1) = u_32_EB*xi_4_EB - u_41_EB;
fval(24, 1) = u_32_CY*xi_4_CY - u_41_CY;
fval(25, 1) = u_41_H; 
fval(26, 1) = u_42_H; 
fval(27, 1) = u_32_B*(1 - xi_4_B) - u_42_B; 
fval(28, 1) = u_32_ST*(1 - xi_4_ST) - u_42_ST - 2000/MW_ST;
fval(29, 1) = u_32_EB*(1 - xi_4_EB) - u_42_EB;
fval(30, 1) = u_32_CY*(1 - xi_4_CY) - u_42_CY;
fval(31, 1) = u_51_H;
fval(32, 1) = u_52_H;
fval(33, 1) = u_41_B*xi_5_B - u_51_B;
fval(34, 1) = u_41_ST*xi_5_ST - u_51_ST;
fval(35, 1) = u_41_EB*xi_5_EB - u_51_EB;
fval(36, 1) = u_41_CY*xi_5_CY - u_51_CY;
fval(37, 1) = u_41_B*(1 - xi_5_B) - u_52_B;
fval(38, 1) = u_41_ST*(1 - xi_5_ST) - u_52_ST;
fval(39, 1) = u_41_EB*(1 - xi_5_EB) - u_52_EB;
fval(40, 1) = u_41_CY*(1 - xi_5_CY) - u_52_CY;
%Feed constraints 
fval(41, 1) = u_O_ST; 
fval(42, 1) = u_0_CY;
fval(43, 1) = u_0_H;
fval(44, 1) = u_0_B + u_52_B - u_1_B;
fval(45, 1) = u_0_EB + u_52_EB - u_1_EB;

%Mole fractions
x_blabla = u_jrrjrk / ejjeor; 
%Calculate new bubble temperature T_51
T_51 = x_blabla*P_s_blabla + ..;
% If old(T_51_new - T_51_old) < 0.1 -> solution is OK
% Else -> Solve again the system of equations with T_51_new

PS: Considérez les équations comme correctes évidemment

Salut,

Je comprends tout à faire le concept mais pas tellement comment l’écrire… J’ai transcrit tout ça dans la matrice qui est donc une $45 \times 45$. J’ai pas encore mis les température $T_{51}$ et $T_{52}$ mais bon pour ce qui en est des équations on en est presque là (ça sera pareil que pour $T_{41}$ et $T_{42}$). À part ça, c’est sûrement plus optimal d’utiliser linsolve mais ça m’a pris du temps vu la taille de la matrice

Ce que je comprends pas c’est quand est-ce que je dois intégrer "mu_old" dans mes équations… Est-ce que c’est au-moment où j’ai une dépendence sur les température ? (donc où j’ai des "xi")

Aussi, je ne vois pas super bien où je pourrais intégrer la boucle $while$. Est-ce que je met tout le système ? J’aurai aussi des termes non-linéaires vu que la température apparaît dans les équations d’Antoine pour le calcul des pressions saturantes (P_s dans le code).

Merci pour ton aide!

%B = Benzene, ST = Styrene, EB = Ethylbenzene, H = Hydrogen, CY =
%Cyclohexane


%Molar masses in kg/mol
MW_B = (78.11)/1000; 
MW_ST = (104.15)/1000;
MW_EB = (106.17)/1000;
MW_H = (2.01)/1000;
MW_CY = (84.16)/1000;

%Conversion rates of rxn 

X_A = 0.5;
X_B = 0.95;

%Antoine parameters 

A_B = 4.73;
B_B = 1316.55;
C_B = - 35.58;
A_EB = 4.08;
B_EB = 1419.30;
C_EB = -60.54;
A_ST = 4.22; 
B_ST = 1525.06; 
C_ST = -56.40;
A_CY = 4.14;
B_CY = 1316.55;
C_CY = - 35.58;

%Splitting ratio xi

%Unit 4
xi_4_lk = 0.99; %EB
xi_4_hk = 0.01; %ST 
P_s_4_lk_B = 10^(-(A_EB - B_EB/(T_42 + C_EB))) ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_4_lk_V = 10^(-(A_EB - B_EB/(T_41 + C_EB)));
P_s_4_hk_B = 10^(-(A_ST - B_ST/(T_42 + C_ST)));
P_s_4_hk_V = 10^(-(A_ST - B_ST/(T_41 + C_ST)));
P_s_4_B_B = 10^(-(A_B - B_B/(T_42 + C_B))); 
P_s_4_B_V = 10^(-(A_B - B_B/(T_41 + C_B))); 
P_s_4_CY_B = 10^(-(A_CY - B_CY/(T_42 + C_CY)));
P_s_4_CY_V = 10^(-(A_CY - B_CY/(T_41 + C_CY)));
P_4 = 0.05; %pressure of the col 4 [bar]

%Define K-value (ie. P_sat / P_column)
K_4_lk_B = P_s_4_lk_B / P_4;
K_4_lk_V = P_s_4_lk_V / P_4;
K_4_hk_B = P_s_4_hk_B / P_4;
K_4_hk_V = P_s_4_hk_V / P_4;
K_4_B_B = P_s_4_B_B / P_4;
K_4_B_V = P_s_4_B_V / P_4;
K_4_CY_B = P_s_4_B_CY_B / P_4;
K_4_CY_V = P_s_4_CY_V / P_4;


alpha_4_lkhk = (K_4_lk_B*K_4_lk_V/(K_4_hk_B*K_4_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_4_min = log(xi_4_lk*(1-xi_4_hk)/xi_4_hk*(1-xi_4_lk))/log(alpha_4_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((K_4_B_B*K_4_B_V)/(K_4_hk_B*K_4_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_4_CYhk = ((K_4_CY_B*K_4_CY_V)/(K_4_hk_B*K_4_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_4_B = ((alpha_4_Bhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_Bhk-1)*xi_4_hk);
xi_4_ST = 0.01
xi_4_EB = 0.99
xi_4_CY = ((alpha_4_CYhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_CYhk-1)*xi_4_hk);


M = 0;
R = 0;

% First equation
M(1, 1) = 5;
M(1, 2) = -1;

% Second equation
M(2, 1) = -MW_B;
M(2, 2) = -MW_EB;
M(2, 3) = -MW_ST;
M(2, 4) = -MW_CY;
M(2, 5) = -MW_H;
M(2, 21) = MW_B;
M(2, 22) = MW_EB;
M(2, 23) = MW_ST;
M(2, 24) = MW_CY;
M(2, 25) = MW_H;
M(2, 31) = MW_B;
M(2, 32) = MW_EB;
M(2, 33) = MW_ST;
M(2, 34) = MW_CY;
M(2, 35) = MW_H;
M(2, 36) = MW_B;
M(2, 37) = MW_EB;
M(2, 38) = MW_ST;
M(2, 39) = MW_CY;
M(2, 40) = MW_H;

% Third equation
M(3, 1) = 1;
M(3, 6) = -1;
M(3, 36) = 1; 

% Fourth Eq

M(4, 2) = 1;
M(4, 7) = -1;
M(4, 37) = 1; 

% Fifth Eq

M(5, 3) = 1;
M(5, 8) = -1;
M(5, 38) = 1; 

% Eq 6

M(6, 4) = 1;
M(6, 9) = -1;
M(6, 39) = 1;

% Eq 7

M(7, 5) = 1;
M(7, 10) = -1;
M(7, 40) = 1;

% Eq 8 

M(8, 7) = (1-X_A); 
M(8, 12) = -1; 

% Eq 9

M(9, 7) = X_A;
M(9, 8) = 1; 
M(9, 13) = -1;

% Eq 10

M(10, 6) = X_B;
M(10, 10) = 1;
M(10, 14) = -1;

% Eq 11

M(11, 6) = -3*X_B;
M(11, 7) = X_A; 
M(11, 10) = 1;
M(11, 15) = -1;

% Eq 12

M(12, 16) = -1; 

% Eq 13

M(13, 17) = -1;

% Eq 14

M(14, 18) = -1;

% Eq 15

M(15, 19) = -1:

% Eq 16

M(16, 15) = 1;
M(16, 20) = -1;

% Eq 17

M(17, 11) = 1; 
M(17, 21) = -1;

% Eq 18

M(18, 12) = 1; 
M(18, 22) = -1;

% Eq 19

M(19, 13) = 1; 
M(19, 23) = -1;


% Eq 20

M(20, 14) = 1; 
M(20, 24) = -1;

% Eq 21

M(21, 21) = xi_4_B;
M(21, 26) = -1;

% Eq 22

M(22, 22) = xi_4_EB;
M(22, 27) = -1;


% Eq 23

M(23, 23) = xi_4_ST;
M(23, 28) = -1;


% Eq 24

M(24, 24) = xi_4_CY;
M(24, 29) = -1;

% Eq 25

M(25, 30) = 1; 

% Eq 26 

M(26, 35) = 1;

% Eq 27

M(27, 21) = (1 - xi_4_B);
M(27, 31) = -1; 

% Eq 28

M(28, 22) = (1 - xi_4_EB);
M(28, 32) = -1; 

% Eq 29

M(29, 23) = (1 - xi_4_ST);
M(29, 33) = -1;
R(29) = 2000/MW_ST;

% Eq 30

M(30, 24) = (1 - xi_4_CY);
M(30, 34) = -1;

% Eq 31 & 32  

M(31, 40) = 1; 
M(32, 45) = 1;

% Eq 33

M(33, 26) = xi_5_B;
M(33, 36) = -1;

% Eq 34

M(34, 27) = xi_5_EB;
M(34, 37) = -1;

% Eq 35

M(35, 28) = xi_5_ST;
M(35, 38) = -1;

% Eq 36

M(36, 29) = xi_5_CY;
M(36, 39) = -1;

% Eq 37

M(37, 26) = (1 - xi_5_B);
M(37, 41) = -1;

% Eq 38

M(38, 27) = (1 - xi_5_EB);
M(38, 42) = -1;

% Eq 39

M(39, 28) = (1 - xi_5_ST);
M(39, 43) = -1;


% Eq 40

M(40, 29) = (1 - xi_5_CY);
M(40, 44) = -1;

% Eq 41, 42, 43, 44 & 45 - Feed constraints 

M(41, 3) = 1;
M(42, 4) = 1; 
M(43, 5) = 1;
M(44, 1) = 1;
M(44, 6) = -1;
M(44, 41) = 1;
M(45, 2) = 1;
M(45, 7) = -1;
M(45, 42) = 1;

% Solve for mu_new
mu_new = linsolve(M, R);

Je comprends tout à faire le concept […] Ce que je comprends pas c’est quand est-ce que je dois intégrer "mu_old" dans mes équations… Est-ce que c’est au-moment où j’ai une dépendence sur les température ? (donc où j’ai des "xi")

Aussi, je ne vois pas super bien où je pourrais intégrer la boucle while. Est-ce que je met tout le système ? J’aurai aussi des termes non-linéaires vu que la température apparaît dans les équations d’Antoine pour le calcul des pressions saturantes (P_s dans le code).

Ces questions me font dire que tu n’as pas compris l’idée derrière l’algorithme. mu_old va en effet apparaître à travers certains des coefficients (les xi qui dépendent des T qui eux même dépendent de mu).

Pour la boucle while, je l’ai écrit dans ce message. Le but, c’est d’itérer tant que les températures n’ont pas convergées. Et oui tu vas avoir des termes non linéaires en mu_old dans la matrice, c’est le but de toute la méthode : transformer une équation non linéaire en une limite d’un processus itératif linéaire.

Si je reprends le code dans mon message précédent, j’aurai donc un while partir de la ligne 33 jusqu’à la ligne 300? Boucle de la forme:

while abs(T51(mu_old) - T51(mu_new)) <0.1

end

Et par exemple, à la ligne 194 où j’ai une dépendance sur T, j’aurai xi_5_EB(mu_old) ?

D’accord. Par contre, je vois pas comment est-ce que je peux récupérer une solution.. Est-ce que c’est par exemple munew(2) qui va me renvoyer ${\mu {0,EB}}$ ?

Je demande car j’en ai besoin pour calculer les fractions molaire que je calculerais donc juste après la ligne mu_new = linsolve(M, R). Et du coup, ma boucle ira jusqu’après le calcul des nouvelles températures (ie. juste après les fractions molaires…)

Par contre, j’avoue que je vois pas du tout comment je pourrais écrire la valeur de xi_5_EB en fonction de mu_old avec du code…

Est-ce que c’est par exemple munew(2) qui va me renvoyer ${\mu {0,EB}}$ ?

Ben ça je peux pas le deviner à ta place, je ne sais pas à quelle position dans le vecteur $\mu$ tu as mis ${\mu {0,EB}}$…

Je demande car j’en ai besoin pour calculer les fractions molaire que je calculerais donc juste après la ligne mu_new = linsolve(M, R). Et du coup, ma boucle ira jusqu’après le calcul des nouvelles températures (ie. juste après les fractions molaires…)

Hmmm, t’aurais pas intérêt à mettre ces calculs là au début de la boucle plutôt qu’à la fin ? Ça t’évite d’avoir à les initialiser.

Par contre, j’avoue que je vois pas du tout comment je pourrais écrire la valeur de xi_5_EB en fonction de mu_old avec du code…

Ben… Ça aussi c’est pas un truc que je peux deviner, c’est toi qui est censé savoir comment ces coefficients se calculent. C’est toi qui m’avait dit que pour une valeur de $\mu$, tu pouvais calculer les $T$ et en déduire les $\xi$. Ben il suffit d’écrire juste ça.

Au fait, je viens de voir mais là :

while abs(T51(mu_old) - T51(mu_new)) <0.1

ta condition est dans le mauvais sens. while, ça veut dire tant que, pas jusqu’à ce que.

Est-ce que c’est par exemple munew(2) qui va me renvoyer ${\mu {0,EB}}$ ?

Ben ça je peux pas le deviner à ta place, je ne sais pas à quelle position dans le vecteur $\mu$ tu as mis ${\mu {0,EB}}$…

Oui oui, je veux dire du point de vue conceptuel. En Matlab, quand c’est un vecteur à une dimension il suffit de faire mu_new(2) et pas mu_new(2,1) ou pas?

Hmmm, t’aurais pas intérêt à mettre ces calculs là au début de la boucle plutôt qu’à la fin ? Ça t’évite d’avoir à les initialiser.

Sûrement Le soucis si je les mets au début c’est que la première fois que la boucle va tourner elle sera perdue car elle aura pas encore de valeurs pour mes mu_new ?

Ben… Ça aussi c’est pas un truc que je peux deviner, c’est toi qui est censé savoir comment ces coefficients se calculent. C’est toi qui m’avait dit que pour une valeur de $\mu$, tu pouvais calculer les $T$ et en déduire les $\xi$. Ben il suffit d’écrire juste ça.

Ah, ça je l’ai fait avant de lister toutes les équations (lignes 35–60 pour l’unité 4 avant d’initialiser M et R à zéro):

%Splitting ratio xi


%Unit 4
xi_4_lk = 0.99; %EB
xi_4_hk = 0.01; %ST 
P_s_4_lk_B = 10^(-(A_EB - B_EB/(T_42 + C_EB))) ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_4_lk_V = 10^(-(A_EB - B_EB/(T_41 + C_EB)));
P_s_4_hk_B = 10^(-(A_ST - B_ST/(T_42 + C_ST)));
P_s_4_hk_V = 10^(-(A_ST - B_ST/(T_41 + C_ST)));
P_s_4_B_B = 10^(-(A_B - B_B/(T_42 + C_B))); 
P_s_4_B_V = 10^(-(A_B - B_B/(T_41 + C_B))); 
P_s_4_CY_B = 10^(-(A_CY - B_CY/(T_42 + C_CY)));
P_s_4_CY_V = 10^(-(A_CY - B_CY/(T_41 + C_CY)));
P_4 = 0.05; %pressure of the col 4 [bar]

%Define K-value (ie. P_sat / P_column)
K_4_lk_B = P_s_4_lk_B / P_4;
K_4_lk_V = P_s_4_lk_V / P_4;
K_4_hk_B = P_s_4_hk_B / P_4;
K_4_hk_V = P_s_4_hk_V / P_4;
K_4_B_B = P_s_4_B_B / P_4;
K_4_B_V = P_s_4_B_V / P_4;
K_4_CY_B = P_s_4_CY_B / P_4;
K_4_CY_V = P_s_4_CY_V / P_4;


alpha_4_lkhk = (K_4_lk_B*K_4_lk_V/(K_4_hk_B*K_4_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_4_min = log(xi_4_lk*(1-xi_4_hk)/xi_4_hk*(1-xi_4_lk))/log(alpha_4_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((K_4_B_B*K_4_B_V)/(K_4_hk_B*K_4_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_4_CYhk = ((K_4_CY_B*K_4_CY_V)/(K_4_hk_B*K_4_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_4_B = ((alpha_4_Bhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_Bhk-1)*xi_4_hk);
xi_4_ST = 0.01
xi_4_EB = 0.99
xi_4_CY = ((alpha_4_CYhk)^(N_4_min)*xi_4_hk)/(1+(alpha_4_CYhk-1)*xi_4_hk);


%Unit 5
xi_5_lk = 0.99; %CY
xi_5_hk = 0.01; %EB 
P_s_5_lk_B = 10^(-(A_CY - B_CY/(T_52 + C_CY))) ;%Antoine Eq for saturation pressure (@ bottom and top of the col)
P_s_5_lk_V = 10^(-(A_CY - B_CY/(T_51 + C_CY)));
P_s_5_hk_B = 10^(-(A_EB - B_EB/(T_52 + C_EB)));
P_s_5_hk_V = 10^(-(A_EB - B_EB/(T_51 + C_EB)));
P_s_5_B_B = 10^(-(A_B - B_B/(T_52 + C_B))); 
P_s_5_B_V = 10^(-(A_B - B_B/(T_51 + C_B))); 
P_s_5_ST_B = 10^(-(A_ST - B_ST/(T_52 + C_ST)));
P_s_5_ST_V = 10^(-(A_ST - B_ST/(T_51 + C_ST)));
P_5 = 1; %pressure of the col 4 [bar] 

%Define K-value (ie. P_sat / P_column)
K_5_lk_B = P_s_5_lk_B / P_5;
K_5_lk_V = P_s_5_lk_V / P_5;
K_5_hk_B = P_s_5_hk_B / P_5;
K_5_hk_V = P_s_5_hk_V / P_5;
K_5_B_B = P_s_5_B_B / P_5;
K_5_B_V = P_s_5_B_V / P_5;
K_5_ST_B = P_s_5_ST_B / P_5;
K_5_ST_V = P_s_5_ST_V / P_5;


alpha_5_lkhk = (K_5_lk_B*K_5_lk_V/(K_5_hk_B*K_5_hk_V))^(1/2) % relative volatility between lk and hk (geometric mean). B = Bottom of the col , V = Vapour (ie top of the col)
N_5_min = log(xi_5_lk*(1-xi_5_hk)/xi_5_hk*(1-xi_5_lk))/log(alpha_5_lkhk); %Min number of stages (to calculate splitting ratios of the other comp.)
alpha_4_Bhk = ((K_5_B_B*K_5_B_V)/(K_5_hk_B*K_5_hk_V))^(1/2);% rel. vol. between B and ST (geom. mean)
alpha_5_SThk = ((K_5_ST_B*K_5_ST_V)/(K_5_hk_B*K_5_hk_V))^(1/2):% rel. vol. between CY and ST (geom. mean)
xi_5_B = ((alpha_5_Bhk)^(N_5_min)*xi_5_hk)/(1+(alpha_5_Bhk -1)*xi_5_hk);
xi_5_EB = 0.01
xi_5_CY = 0.99
xi_5_ST = ((alpha_5_SThk)^(N_5_min)*xi_5_hk)/(1+(alpha_5_SThk -1)*xi_5_hk);

Là, j’ai déjà défini ce qu’était ces $\xi$. Mon but juste avant ces quelques lignes est de donner des températures (notées T_52, T_51, … dans le code) approximatives pour les faire converger par la suite vers les températures réelles. Mais dans tout ça j’ai jamais dans mon code un "mu_old" ce qui me perturbe énormément pour l’écriture de la boucle.

Ah merci! Je dois donc changer le sens de l’inégalité :-)

Le soucis si je les mets au début c’est que la première fois que la boucle va tourner elle sera perdue car elle aura pas encore de valeurs pour mes mu_new ?

Rien compris, lors du premier tour de boucle, mu_new n’a pas encore été calculé une seule fois, comment diable compte tu te servir de sa valeur ? Au deuxième tour de boucle, le mu_old sera le mu_new du tour précédent et sera utilisé à ce moment là.

Ah, ça je l’ai fait avant de lister toutes les équations (lignes 35–60 pour l’unité 4 avant d’initialiser M et R à zéro):

Elles sortent d’où T_41 et T_42 ? Et les valeurs des différents xi ? De ce que j’avais compris, c’était ça qui dépendait de $\mu$ (et donc dans ta boucle, de mu_old). Si ça dépend pas de $\mu$, alors j’ai absolument rien compris à ton problème depuis le départ et j’ai joyeusement écrit des $M(\mu)$ alors que ta matrice ne dépend pas de $\mu$…

Rien compris, lors du premier tour de boucle, mu_new n’a pas encore été calculé une seule fois, comment diable compte tu te servir de sa valeur ? Au deuxième tour de boucle, le mu_old sera le mu_new du tour précédent et sera utilisé à ce moment là.

Mais… alors, faut bien que je le mette à la fin vu que j’ai aucune idée des "mu_old"… J’aurai une idée de leur valeur uniquement après la première itération.

T41 et T_42 (idem pour T_51 et T_52) sont des températures de ces fameuses unités de séparation. C’est pour ces paramètres que je "devine" les valeurs… Les $\xi _i$ (xi) ce sont des formules de mon cours pour ces unités de séparation. Ces $\xi _i$ ne dépendent pas des flux $\mu$ mais uniquement du paramètre $N{min}$ (calculé juste avant par une autre formule) et de la volatilité relative (ie. de la température -> T_41 et T_42 pour l’unité 4)

Pour répéter, tous les $\mu$ sont mes inconnues mais les quatre températures $T$ le sont aussi (même si je les devine au début, à la fin elles auront des valeurs différents en fonction de la convergence).

J’espère que c’est plus clair…

Mais… alors, faut bien que je le mette à la fin vu que j’ai aucune idée des "mu_old"… J’aurai une idée de leur valeur uniquement après la première itération.

Ben… Non, tu te donnes un mu_old pour commencer l’itération. Enfin, vu ce que tu dis ensuite, j’ai probablement toujours pas compris correctement le problème.

T_41 et T_42 (idem pour T_51 et T_52) sont des températures de ces fameuses unités de séparation. C’est pour ces paramètres que je "devine" les valeurs… Les $\xi _i$ (xi) ce sont des formules de mon cours pour ces unités de séparation.

Je sais pas ce que sont ces "fameuses unités de séparation". Encore une fois, je comprends pas ton problème chimique, mais si tu me donnes un problème mathématiquement bien posé on arrivera à avancer. Quand tu parles de "formules" dans tous les sens, ça m’avance pas parce que ça me dit pas qui dépend de qui. J’ai juste besoin de savoir qui est fonction de quoi dans ton problème, et si ces fonctions sont linéaires ou non.

Ces $\xi _i$ ne dépendent pas des flux $\mu$ mais uniquement du paramètre $N_{min}$ (calculé juste avant par une autre formule) et de la volatilité relative (ie. de la température -> T_41 et T_42 pour l’unité 4)

Pareil, tu parles de formule mais sans me dire qui dépend de quoi je suis dans le flou le plus complet.

Pour répéter, tous les $\mu$ sont mes inconnues mais les quatre températures $T$ le sont aussi (même si je les devine au début, à la fin elles auront des valeurs différents en fonction de la convergence).

Comme t’as dit trois trucs différents au moins, je me suis un peu perdu en route.

Si j’essaye de re-résumer, tu as des flux (un vecteur de $\mu$) inconnus reliés par plein d’équations qui peuvent se mettre sous forme matricielle. Les coefficients de cette matrice dépendent des coefficients $\xi$ qui je comprends bien sont eux-même des fonctions des températures qui sont aussi des inconnues du problème. Jusque là j’ai bon, ou il me manque un truc ?

Maintenant, ce qui me pose question, c’est comment tu calcules ces températures, de quelles autres variables elles dépendent. Sans lien entre les températures et les flux je vois mal comment tu comptes t’en sortir.

Ben… Non, tu te donnes un mu_old pour commencer l’itération. Enfin, vu ce que tu dis ensuite, j’ai probablement toujours pas compris correctement le problème.

Je peux mettre n’importe quoi ? Je peux pas initialiser à zéro car je vais devoir faire un ratio pour obtenir une fraction (molaire).

Je sais pas ce que sont ces "fameuses unités de séparation". Encore une fois, je comprends pas ton problème chimique, mais si tu me donnes un problème mathématiquement bien posé on arrivera à avancer. Quand tu parles de "formules" dans tous les sens, ça m’avance pas parce que ça me dit pas qui dépend de qui. J’ai juste besoin de savoir qui est fonction de quoi dans ton problème, et si ces fonctions sont linéaires ou non.

OK. Désolé. Du coup, sans parler chimie mais simplement mathématiquement voici:

• Toutes les équations que j’ai écris dans le code M(.., ..) sont linéaires et sont simplement des conservations de masse;

• Cependant, dans ces équation on retrouve un paramètre qu’il faut évaluer avec des équations non-linéaire. C’est les fameux $\xi$ (y en a plusieurs mais ils sont évalués avec les mêmes équations - simplement les valeurs changent).

• Ces fameuses équations pour trouver $\xi$ sont les suivantes:

Étape 1: on cherche un paramètre appelé N_min

${N_{\min }} = \ln \left[ {\frac{{{\xi _{lk}}(1 - {\xi _{hk}})}}{{{\xi _{hk}}(1 - {\xi _{lk}})}}} \right]/\ln {\alpha _{lk,hk}}$

où ${\overline \alpha _{ij}} = \sqrt {\alpha _{ij}^B\alpha _{ij}^V} = \sqrt {\frac{{{P_{s,i}}({T^B}){P_{s,i}}({T^V})}}{{{P_{s,j}}({T^B}){P_{s,j}}({T^V})}}}$.

On évalue les $P_{s,i}$ avec une autre équation non-linéaire qui est:

${P_s} = {10^{ - (A - \frac{B}{{C + T}})}}$

où la température apparaît!

Étape 2: on calcule le $\xi$ correspondant avec l’équation non-linéaire suivante:

${\xi _j} = \alpha _{j,hk}^{{N_{\min }}}{\xi _{hk}}/\left[ {1 + (\alpha _{j,hk}^{{N_{\min }}} - 1){\xi _{hk}}} \right]$

où $\xi_{hk}$ est une constante (écrit dans le code)

C’est la que j’ai besoin des fractions molaires (donc des mu). Chaque nouvelle température sera calculée à partir de là avec une relation (que j’ai pas écrite car j’en suis pas encore sûr.. je dois demander vérification).

J’espère que c’est (finalement) clair!