A mathematical model for the computation of the phase equilibrium related to atmospheric organic aerosols is proposed. The equilibrium is given by the minimum of the Gibbs free energy and is characterized using the notion of phase simplex of its convex hull. A primal-dual interior-point method solving the Karush-Kuhn-Tucker conditions is detailed. Numerical results show the efficiency of our algorithm

Alexandre Caboussat

Jiwen He

John H. Seinfeld

Kee-youn Yoo

Neal R. Amundson

Kee-Youn Yoo

Comptes Rendus Mathématique

1/riced. Theand solidthedergie delides. Nousraissant à),versity ofC. R. Acad. Sci. Paris, Ser. I 340 (2005) 683–686http://france.elsevier.com/direct/CRASSNumerical AnalysisAn optimization problem related to the modeling of atmospheinorganic aerosolsNeal R. Amundsona, Alexandre Caboussata,1, Jiwen Hea,2, John H. Seinfeldb,Kee-Youn Yooca Department of Mathematics, University of Houston, Houston, TX 77204, USAb Department of Chemical Engineering, California Institute of Technology, Pasadena, CA 91125, USAc Department of Chemical Engineering, Seoul National University of Technology, Seoul 139-743, South KoreaReceived 31 December 2004; accepted 25 January 2005Presented by Roland GlowinskiAbstractA mathematical model for the computation of chemical equilibrium of atmospheric inorganic aerosols is proposequilibrium is given by the minimum of the Gibbs free energy for a system involving an aqueous phase, a gas phasesalts. A primal-dual method solving theKarush–Kuhn–Tuckerconditions is detailed. An active set/Newton method permitscomputation of the minimum and track solid salts at the equilibrium.To cite this article: N.R. Amundson et al., C. R. Acad.Sci. Paris, Ser. I 340 (2005). 2005 Académie des sciences. Published by Elsevier SAS. All rights reserved.RésuméUn problème d’optimisation lié à la modélisation d’aérosols inorganiques. Nous proposons un modèle pour l’étudel’équilibre chimique d’un aérosol inorganique dilué dans l’air. L’état d’équilibre est caracterisé par le minimum d’éneGibbs pour un système chimique faisant intervenir une phase liquide, une phase gazeuse et plusieurs phases soprésentons une méthode pour la résolution des conditions deKarush–Kuhn–Tucker. Nous utilisons un algorithme de typeactiveset, couplé avec une méthode de Newton pour déterminer le minimum d’énergie ainsi que les phases solides appal’équilibre. Pour citer cet article : N.R. Amundson et al., C. R. Acad. Sci. Paris, Ser. I 340 (2005). 2005 Académie des sciences. Published by Elsevier SAS. All rights reserved.E-mail addresses:amundson@uh.edu (N.R. Amundson), caboussat@math.uh.edu (A. Caboussat), jiwenhe@math.uh.edu (J. Heseinfeld@caltech.edu (J.H. Seinfeld), kyyoo@snut.ac.kr (K.-Y. Yoo).1 Supported by the Swiss National Science Foundation, Grant PBEL2-103152.2 Supported by the United States Environmental Protection Agency through Cooperative Agreement CR-829068-01 to the UniHouston.1631-073X/$ – see front matter 2005 Académie des sciences. Published by Elsevier SAS. All rights reserved.doi:10.1016/j.crma.2005.01.025684 N.R. Amundson et al. / C. R. Acad. Sci. Paris, Ser. I 340 (2005) 683–686eloped toriori andey oftenr qualityf gas–tainedthat isactivem set ofpressured thestemeouscenhe solidTuckerds tole withspeciesde) isOver the last two decades, a series of thermodynamic modules, see for instance [1,5], has been devpredict the phase transition of atmospheric inorganic aerosols. However, these modules either rely on a pincomplete knowledge of the presence of solid phases or are computationally intensive, implying that thfail to accurately predict the existing solid phases at the equilibrium or are useless for three dimensional aimodels. In this Note, a numerical method for the minimization of the Gibbs free energy in the framework oliquid and solid–liquid equilibrium reactions [6] is presented. The Karush–Kuhn–Tucker (KKT) points are obwith a primal–dual algorithm, where a Newton iteration is applied to the reduced KKT system of equationsprojected on an ‘active’ set of solid phases to find the next primal–dual approximation of the solution. Theset strategy permits to add/delete constraints to/from a working set of active constraints until the equilibriusolid phases is obtained.The chemical equilibrium problem for a closed inorganic aerosol system at constant temperature andand for a specified feed vectorb is the solution of the minimization problemminimize G(nl ,ns ,ng) = nTl µl + nTgµg + nTs µs ,subject to nl > 0, ng > 0, ns  0,Alnl + Agng + Asns = b,(1)wherenα ∈ Rmα , µα ∈ Rmα andAα ∈ Rmc×mα are the concentration vector, the chemical potential vector ancomponent-based formula matrix for species setα, for α = l, g, s respectively. The subscriptsl, g, s denote theliquid, gas and solid phases, respectively. Heremc denotes the number of independent components in the syandmα denotes the number of species in species setα.In (1), the chemical potential vectorµl , is defined as the gradient of the Gibbs free energy of the aquphaseGl :Rml+ → R+, so that,∀nl > 0, µl = ∇Gl(nl). The Gibbs free energyGl is assumed to be at least twicontinuously differentiable in the interior ofRml+ . The first-order homogeneity ofGl is the basis for the relatioGl(nl ) = nTl µl , and theGibbs–Duhem relations∇2Gl(nl)nl = 0. The chemical potentialsµs are constant. TheHessian matrices for the liquid and gas species are respectivelyHl = ∇2Gl(nl ) andHg = diag(1/ng). The aqueousand gas phases are assumed to be present at the equilibrium. The difficulty in solving (1) is to identify tphases existing at the equilibrium.Taking into account the Gibbs–Duhem relations, the solution of (1) is characterized by the Karush–Kuhn–(KKT) system of the first order necessary optimality conditions:µl + ATl λ = 0,µg + ATgλ = 0,ns  0, µs + ATs λ  0, nTs (µs + ATs λ) = 0,Alnl + Agng + Asns = b.(2)A primal–dual solution of the KKT system (2), generally non unique, is called a KKT point. One neeperform phase stability analysis to determine whether a postulated KKT point is thermodynamically stabrespect to any perturbation innl , ng andns . Let (n†l ,n†g,n†s ,λ†) denotes a stable KKT point of (2) and letI †s :={i ∈ {1, . . . , ns}: n†s,i > 0} denotes the corresponding active set of solid phases, which is the set of the solidactually formed in the system at equilibrium. The cardinal ofI †s is denoted byms . Let Ās be the matrix composeby the columns ofAs that are inI †s andĀ = [Al ,Ag, Ās]. The matrixĀ is assumed to be of full row rank. Thstability analysis is then based on thesecond order sufficient condition, which states that,pTHp is positive definite for all nonzero vectorp such thatĀp = 0, (3)whereH = ∇2nl ,ng,ns G(n†l ,n†g ,n†s ) is the Hessian matrix of the Gibbs free energy of the system.Recall that the inertia of a symmetric matrix is an ordered triple(i+, i−, io), wherei+ is the number of positiveeigenvalues,i− the number of negative eigenvalues, andi0 the number of zero eigenvalues. Relationship (3equivalent to requiring the so-called KKT matrixK = ( H ĀT), to have a certain inertia.Ā 0N.R. Amundson et al. / C. R. Acad. Sci. Paris, Ser. I 340 (2005) 683–686 685esncurrentyxsnhed(2)hatise,t. TheLemma 1. If the matrixĀ is of full rank, the condition(3) is equivalent toinertia(ZTĀHZĀ) = (ml + mg + ms,0,0), (4)whereZĀ is a null-space matrix for̄A.Based on an inertia result of Gould [3], the inertia ofK is given by inertia(K) = inertia(ZTĀHZĀ)+ (mc,mc,0).Hence (3) implies thatK is invertible.In the active set algorithm, a sequence of guesses of the active set of solid phasesI †s i constructed. Thissequence of the so-called working sets, denoted byIs , is defined in thedualsense byIs := {j ∈ {1, . . . ,ms}: (µs +ATs λ)j = 0}. The dual variableλ, together with the primal variable(nl ,ng,ns) consists of a sequence of iteratthat converges to the optimal primal–dual solution(n†l ,n†g ,n†s ,λ†), as the setIs , starting from an approximatioof I †s , is expected to converge toI †s [4].Along the process of applying the active-set strategy, the KKT system (2) is first projected onto theworking setIs to form a reduced KKT system. The concentration vectorns is decomposed intons = (n̄s , n̂s)wheren̂s = 0 andn̄s > 0. Let Ās be the matrix composed by the columns ofAs with indices inIs and Âs thematrix composed by the columns ofAs which are not inĀs . Similarly, defineµ̄s andµ̂s as the vector composed bthe components ofµs with indices that are inIs and that are not inIs respectively. LetZĀs be a null-space matriof Ās , Azl = ZĀs Al andAzg = ZĀs As . Let us also defineη by λ = λ∗ + ZĀs η, whereλ∗ is a particular solutionsatisfyingµ̄s + Āsλ∗ = 0 andµ̂s + Âsλ∗  0. The system (2) projected onIs consists of the equality constraintµl + ATzlη = 0, µg + ATzgη = 0, Azlnl + Azgng = ZĀs b, (5)together with inequalities constraintsµ̂s + ÂTs λ  0. (6)Let us ignore for the moment the fact thatn̄s must remain non-negative and apply oneN wton iterationto thereduced system (5) to find the nextprimal–dualapproximation(nl ,ng,η) of the solution. The resulting Newtosystem for the displacementspl , pg , andpη is Hl 0 ATzl0 Hg ATzgAzl Azg 0( plpgpη)=( blbgbη), (7)wherebl = −(µl +ATzlη), bg = −(µg +ATzgη) andbη = ZĀs b−Azlnl −Azgng . The displacement inλ is obtainedas a displacement in the null-space, defined bypλ = ZĀs pη. The displacementpλ is restricted to a certain lengtso that the nextλ+ stays feasible with respect to (6). Finally, the next working setI+s is obtained by addingconstraints that are encountered by the newλ+, i.e. I+s = Is ∪ Ias , with Ias = {j /∈ Is : (µ̂s + ÂTs λ+)j = 0}.Constraints are added in the working set until the sequence of(nl ,ng,λ) converges to a solution of the reducKKT system (5). The concentrations¯ s in the working setIs are then computed based on the last equality ofvia n̄s = (ĀTs )−1(b − Alnl − Agng) where(ĀTs )−1 is the left pseudo-inverse of̄ATs . The non-negativeness ofn̄sis enforced byremoving a saturated saltfrom the working setIs when its concentration becomes negative, tis I+s = Is\Ids , with Ids = {j ∈ Is : ns,j < 0}. If Ids = ∅, a new loop of Newton iterations is restarted; otherwthe algorithm stops at a stable KKT point and the active set of solid phasesI †s at the equilibrium is obtained. Iis assumed that̄As remains of full rank by addition or deletion of constraints in the working set, see also [2]solvability of the linear system (7) and the condition (3) are ensured by the following lemma.Lemma 2. The condition(3) is equivalent toZT(H + S )Z > 0, (8)l l l l686 N.R. Amundson et al. / C. R. Acad. Sci. Paris, Ser. I 340 (2005) 683–686ht:uilibrium.evolution) (1999)33 (1)problem,5–52.gional air. Chem.Fig. 1. Modeling of a sulfate aerosol. Left: reconstruction of the phase diagram at 25◦C with tracking of the presence of solid phases. Rigevolution of the particle mass in function of therelative humidityRH for three different feed vectorsb.whereSl = ATzlQgR−Tg HgR−1g QTgAzl andZl = (Q̃l ,QlR−Tl Qg) with (Qg Q̃g), Rg and (Ql Q̃l), Rl being the QRfactorization ofAzg andATzl respectively. Under condition(8), the system(7) is solvable.To ensure that the primal–dual algorithm converges to aminimum or stable equilibrium, the condition (3) isenforced via (8), by replacingHl with a modified matrix̃Hl so thatZTl (H̃l + Sl)Zl > 0.To illustrate the efficiency of the algorithm, a numerical example of atmospheric sulfate aerosol (NH4)2SO4–H2SO4–H2O is considered here. An aqueous phase is present and three solid phases may appear at eqFig. 1 shows, on the left side, the reconstructed sulfate aerosol phase diagram, and, on the right side, theof the particle mass in function of ther lative humidityRH for three different feed vectorsb. It illustrates that themethod allows us to track phase changes very accurately.References[1] A.S. Ansari, S.N. Pandis, Prediction of multicomponent inorganic atmospheric aerosol behavior, Atmospheric Environment 33 (5745–757.[2] P.E. Gill, W. Murray, M.A. Saunders, M.H. Wright, Inertia-controlling methods for general quadratic programming, SIAM Rev.(1991) 1–36.[3] N.I.M. Gould, On practical conditions for the existence and uniqueness of solutions to the general equality quadratic programmingMath. Programming 32 (1985) 90–99.[4] K. Kunisch, F. Rendl, An infeasible active set method for quadratic problems with simple bounds, SIAM J. Optim. 14 (1) (2003) 3[5] A. Nenes, S.N. Pandis, C. Pilinis, Continued development and testing of a new thermodynamic aerosol module for urban and requality models, Atmospheric Environment 33 (10) (1999) 1553–1560.[6] K.-Y. Yoo, J.-W. He, N.R. Amundson, Canonical form and mathematical interpretation of electrolyte solution systems, Korean JEng. 21 (2) (2004) 303–307.

An optimization problem related to the modeling of atmospheric inorganic aerosols

Amundson, Neal R.

Caboussat, Alexandre

He, Jiwen

Seinfeld, John H.

Yoo, Kee-Youn

Numérisation de Documents Anciens Mathématiques

http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.119.5935

d’Aérosols Organiques.

d’Aérosols Organiques.

Abstract

Similar works

Full text

Available Versions

Comptes Rendus Mathématique

Numérisation de Documents Anciens Mathématiques