Modelování 3: Soustava nelineárních algebraických rovnic

V této části pokračujeme s algebraickými rovnicemi a vysvětlíme jak MATLABem řešit soustavu nelineárních rovnic. Postup řešení ilustrujeme na příkladu rovnováhy kapalina/pára ve třísložkové směsi.

Soustava algebraických rovnic v obecném tvaru

Navážeme na minulý díl o řešení jedné nelineární rovnice a rozšíříme naše dovednosti na soustavu více algebraických rovnic. Každý model, jehož popis vede na soustavu N algebraických rovnic můžeme zapsat v obecném tvaru takto

$\begin{aligned} F_1(x_1,x_2,\dots,x_N) &= 0\\ F_2(x_1,x_2,\dots,x_N) &= 0\\ \dots & \\ F_N(x_1,x_2,\dots,x_N) &= 0 \end{aligned}$

Pokud ovládáte postup při řešení jedné algebraické rovnice, nemělo by vám pochopení postupu řešení soustavy rovnic dělat potíže. Tam, kde u jedné rovnice byla jedna neznámá x, máme nyní vektor N neznámých x₁ až x_N. Místo jedné funkce f(x) teď máme N funkcí F₁ až F_N, jejichž hodnota může záviset na všech neznámých x₁ až x_N. Soustava rovnic je splněna, pokud hodnota všech funkcí F (tj. funkcí residuálu) se blíží nule. Soustava algebraických rovnic zapsaná pomocí vektorů bude

$\mathbf{F}(\mathbf{x}) = \mathbf{0}$ ,

kde $\mathbf{x}$ je vektor neznámých x₁ až x_N a $\mathbf{F}(\mathbf{x})$ je vektorová funkce, jejímž výsledkem je vektor residuálů F₁ až F_N. Nakonec $\mathbf{0}$ je vektor, jehož všechny prvky jsou nuly.

Postup řešení a struktura programu zůstane podobná s tím rozdílem, skaláry používané při řešení jedné rovnice nahradíme u soustavy rovnic vektory. To znamená, že

místo proměnné x (neznámá veličina) budeme používat proměnnou xvec(vektor neznámých)
“modelová” funkce nebude ve výstupní proměnné vracet jako residuál jedno číslo, ale sloupcový vektor N reziduálů,
vstupní parametr, kterým se numerické funkci zadává počáteční odhad řešení nebude skalár, ale vektor.

Také budeme volat jinou numerickou funkci než fzero.

Modelování rovnováhy kapalina-pára vícesložkové směsi

Stavové proměnné

Při odvozování matematického modelu libovolného systému budeme často používat pojem stavové proměnné a proto je důležité si teď na začátku ujasnit, co přesně si budeme pod tímto pojmem představovat. Stavové proměnné u modelu, ve smyslu s jakým s nimi budeme v předmětu “Modelování” pracovat, definujeme tak, že se jedná o vlastnosti (parametry) systému, pomocí kterých systém přesně (kompletně) určíme. Pokud známe stav systému (tj. všechny jeho stavové proměnné), můžeme z těchto stavových proměných vypočítat všechny ostatní požadované vlastnosti. Například u N-složkové směsi konkrétních látek, které mohou být přítomné v kapalné i parní fázi potřebujeme k úplnému určení stavu tyto stavové proměnné

teplotu, t,
tlak, p,
N-1 molárních zlomků jednotlivých složek v kapalné fázi, x_i, (tj. složení kapaliny) a
N-1 molárních zlomků jednotlivých složek v plynné fázi, y_i (tj. složení páry).

Celkový počet stavových proměnných je tedy 2N.

Všimněme si, že za stavové proměnné v tuto chvíli nepovažujeme molární zlomky N-té složky, protože k určení složení fáze je nepotřebujeme a v případě potřeby je můžeme dopočítat z molárních zlomků ostatních složek dopočítáním součtu zlomků do jedničky.

Mezi stavovými proměnnými mohou být vzájemné závislosti. Jestliže existuje rovnováha mezi kapalnou a parní fází, znamená to že chemický potenciál každé složky musí být v obou fázích stejný. To poskytuje N rovnic (vazných podmínek) typu $\mu_i^\mathrm{(l)}=\mu_i^\mathrm{(g)}$ , kde chemické potenciály závisí na teplotě, tlaku a složení příslušné fáze. V případě předpokladu ideálního chování systému se z podmínky rovnosti chemických potenciálů odvodí Raoultův zákon

$py_i=p^\circ_i(t)x_i$ .

Pro 2N stavových proměných platí N podmínek (rovnic), takže systém má N stupňů volnosti. To znamená, že prvních N stavových proměnných si můžeme zvolit (parametry modelu) a touto volbou zároveň určíme i zbývající stavové proměnné (neznámé stavové proměnné). K získání hodnot neznámých proměnných ale potřebujeme vyřešit soustavu rovnic.

Kromě stavových proměnných se v modelu budou vyskytovat také tzv. odvozené proměnné, což jsou vlastnosti systému, které můžeme explicitně určit ze znalostí stavových proměnných (parametrů modelu a neznámých stavových proměnných). Odvozenou proměnnou v našem případě bude například tenze par $p_i^\circ$ nebo molární zlomek poslední složky.

Proměnné můžeme tedy rozdělit do těchto kategorií:

Parametry modelu: Část stavových proměnných, které jsou zadány a známé. V programu jsou to konstanty, které zadáváme ve vektoru parametrů pars nebo jsou jejich hodnoty přímo zapsány v kódu funkce.
Neznámé stavové proměnné: Zbývající část stavových proměnných, jejichž hodnotu chceme nalézt řešením soustavy rovnic. Počet těchto proměnných musí být stejný jako počet rovnic v soustavě. V programu budeme tyto proměnné přenášet ve vektoru neznámých xvec. (Někdy budeme neznámé stavové proměnné označovat zkráceně jako stavové proměnné.)
Odvozené proměnné: Ostatní proměnné, jejichž hodnoty v programu explicitně vypočítáme z ostatních proměnných. Jestliže je výpočet složitější, napíšeme pro výpočet zvláštní funkce.

Zadání příkladu

Vypočítejte zbývající vlastnosti (stavové proměnné) u rovnováhy kapalina/pára směsi hexan/heptan/oktan pro tlak 1 atm, teplotu 90°C a molární zlomek hexanu v kapalině x₁=0,20. Předpokládejte platnost Raoultova zákona pro výpočet rovnováhy.

Začneme tím, že si vypíšeme seznam parametrů a neznámých

Parametry modelu: tlak, teplota a molární zlomek hexanu v kapalině (x₁)
Neznámé stavové proměnné: molární zlomek heptanu v kapalině (x₂), molární zlomek hexanu v páře (y₁) a molární zlomek heptanu v páře (y₂)
Seznam rovnic: Pro každou složku můžeme napsat podmínku rovnováhy kapalina/pára ve tvaru Raoultova zákona. Budeme potřebovat funkci pro výpočet tenze par složky v závislosti na teplotě.

Funkce pro výpočet vektoru reziduálu

Budeme potřebovat funkci, která po zadání vektoru stavových proměnných vypočítá vektor reziduálů. Soubor uložíme do aktuálního adresáře pod jménem res_lg.m

function res = res_lg(xvec, pars)
%RES_LG  Residual function for liquid/vapour equilibrium problem
%  Liquid/vapour equilibrium calculation hexane/heptane/octane mixture
%    xvec is vector of uknown state variables
%      xvec(1) - liquid mole fraction heptane (x2)
%      xvec(2) - vapour mole fraction hexane  (y1)
%      xvec(3) - vapour mole fraction heptane (y2)
%    pars is vector of parameters
%      pars(1) - pressure in kPa             (patm)
%      pars(2) - temperature in Celsius      (t)
%      pars(3) - liquid mole fraction hexane (x1)
%    Raoult law for liquid/vapour equilibrium assumed
 
% unpack model parameters
patm = pars(1);    % pressure in kPa
t    = pars(2);    % temperature in Celsius
x1   = pars(3);    % liquid mole fraction hexane
 
% unpack state variables (unknowns)
x2   = xvec(1);    % liquid mole fraction heptane
y1   = xvec(2);    % vapour mole fraction hexane
y2   = xvec(3);    % vapour mole fraction heptane
 
% calculate auxilary variables - last component mole fractions
x3 = 1-x1-x2;      % liquid mole fraction octane
y3 = 1-y1-y2;      % vapour mole fraction octane
 
% calculate auxilary variables - vapour pressures
[A,B,C] = ptension_coef(); % get Antoine equation coefficients
p1 = ptension(A(1), B(1), C(1), t); % hexane vapour pressure in kPa
p2 = ptension(A(2), B(2), C(2), t); % heptane vapour pressure in kPa
p3 = ptension(A(3), B(3), C(3), t); % octane vapour pressure in kPa
 
% calculate the vector of residuals - implement Raoult's law eqauations
res(1) = patm*y1 - p1*x1;
res(2) = patm*y2 - p2*x2;
res(3) = patm*y3 - p3*x3;
 
% fsolve function expects output to be a column vector
res = res';
 
end

K výpisu je potřeba doplnit několik poznámek.

Obě vstupní proměnné naší funkce jsou vektory. Vektor neznámých xvec pří volání z numerické funkce je sloupcový vektor. Ve vlastním programu nejdříve pro přehlednost překopírujeme jednotlivé prvky vektorů do vhodně pojmenovaných proměnných (rozbalení vektoru parametrů a vektoru neznámých).
Zvolili jsme si, že molární zlomky oktanu nebudou stavové proměnné, ale tzv. odvozené proměnné. Můžeme je snadno vypočítat ze zbývajích molárních zlomků. Podmínka, že součet molárních zlomků je roven jedné tak bude vždy splněna.
Tenzi par vypočítáme jako obvykle z Antoineho rovnice. Opět k tomu použijeme dříve vytvořené funkce. Už víte, že soubory ptension.m a ptension_coef.m musíme nakopírovat do pracovního adresáře.
Z každé rovnice vypočítáme residuál, který posupně ukládáme do vektoru res. Počet prvků vektoru musí být stejný jako počet prvků vektoru neznámých. Na pořadí rovnic nezáleží.
Protože numerická funkce očekává ve výstupní proměnné (vektor reziduálů) sloupcový vektor musíme nakonec transponovat řádkový vektor na sloupcový pomocí apostrofu ‘

Než budeme pokračovat dále, tak si nově vytvořenou funcki vyzkoušíme. Na příkazové řádce uložte parametry modelu do proměnné pars a vyzkoušejte zavolat funkci reziduálů pro několik různých hodnot neznámých stavových proměnných, například

>> pars = [101.325, 90, 0.2] >> res_lg([0.3, 0.3, 0.3], pars) >> res_lg([0.3, 0.1, 0.5], pars)

Ověřte, že funkce nehlásí žádné chyby a vrácený reziduál je sloupcový vektor. Hodnoty vektoru reziduálu se vám nepodaří vynulovat, dokud “netrefíte” takovou kombinaci hodnot neznámých, která bude zároveň řešením soustavy. Samozřejmě, že na hledání řešení použijeme funkci z numerické knihovny. Můžeme si představit, že příslušná numerická funkce nedělá nic jiného, než že podle nějakého algoritmu zkouší zadávat různé hodnoty neznámých tak, aby se vracené hodnoty ve vektoru reziduálů zmenšovaly.

Volání numerické funkce

Zbývá ukázat, jak volat numerickou funkci, která vyřeší soustavu algebraických rovnic. Bohužel taková funkce není v základní verzi MATLABu, ale je potřeba mít rozšíření Optimization Toolbox. Seznam všech rozšíření ve vaší instalaci MATLABu zobrazíte příkazem ver

>> ver

Pokud se Optimization Toolbox nachází ve vaší instalaci MATLABu, je možné použít funkci fsolve z této knihovny.

Jestliže Optimization Toolbox na vašem počítači není instalován, tak je možné použít volně šiřitelnou funkci mmfsolve. Tato funkce není součástí MATLABu, takže je nutné ji nejdříve stáhnout a soubor mmfsolve.m nakopírovat do pracovního adresáře.

Základní používání obou funkcí fsolve a mmfsolve, které si ukážeme je stejné. Pokud je potřeba měnit výchozí nastavení, tak postup se u obou funkcí liší a obě mají různé možnosti nastavení. V takovém případě je potřeba načíst příslušnou dokumentaci k funkci. Z vlastní zkušenosti mohu uvést, že fsolve je trochu více robustnější a s jeho použitím je větší šance na dosažení konvergence řešení u komplikovaných systémů. Nicméně pro velkou většinu “rozumně” naformulovaných modelů, které si ukážeme v předmětu “Modelování”, by mmfsolve mělo fungovat stejně dobře jako fsolve a budeme je proto používat i v našich ukázkách.

K hledání řešení úlohy potřebujeme opět počáteční odhad, kterým je vektor x0. Řešení uložíme do vektoru xsol.

>> pars = [101.325, 90, 0.2] >> x0 = [0.5 0.3 0.3] >> xsol = mmfsolve(@(x) res_lg(x,pars), x0)

U soustavy rovnic je opravdu důležité ověřit, zda metoda zkonvergovala a nabízené řešení vyhovuje zadaným rovnicím

>> res_lg(xsol, pars)

Jestliže úloha nezkonvegovala, pak hledání příčiny může být mimořádně obtížné a frustrující. Je možné, že počáteční odhad je příliš daleko od řešení, potom se vyplatí vyzkoušet jiný odhad. Nebo může být chyba ve funkci residuálů a v důsledku této chyby nebude mít soustava rovnic řešení.

Tento díl ukončíme ukázkou hlavního programu (skriptu) s úplným řešením.

% Main script for liquid/vapour equilibrium calculation problem
% vector of parameters             [p,  t,  x1]
% vector of uknown state variables [x2, y1, y2]
 
% set model parameters
p  = 101.325; % kPa
t  = 90;      % Celsius
x1 = 0.20;    % hexane liquid mole fraction
pars = [p, t, x1] % pack vector of parameters
 
% guess unknowns
x2 = 0.5;
y1 = 0.3;
y2 = 0.3;
xvec0 = [x2, y1, y2] % pack vector of initial guess
 
% solve set of nonlinear equations by "fsolve" (or "mmfsolve")
xsol = mmfsolve(@(x) res_lg(x,pars), xvec0)
 
% check that solution has converged ("res" should be close to zero)
res = res_lg(xsol, pars)
 
% unpack and show the solution
x2 = xsol(1);
y1 = xsol(2);
y2 = xsol(3);
x3 = 1-x1-x2;
y3 = 1-y1-y2;
x_fractions = [x1,x2,x3]
y_fractions = [y1,y2,y3]

A to jsou pro začátek všechny potřebné informace k řešení soustavy algebraických rovnic a v dalších dílech je můžeme začít využívat pří vytváření modelů různých zařízení.

Print This Post